geoit 0.0.2

//! Radical extension tower: exact arithmetic for elements of ℚ(α₁, ..., αₖ)
//! where each αᵢ = ⁿ√(expression in prior αs).
//!
//! Replaces `QuadraticSurd` and `AlgebraicNumber` with a single unified type
//! that handles √2 + √3, ∛2, and arbitrary nested radicals without panics.
//!
//! # Design
//!
//! A `RadicalTower` specifies a chain of field extensions:
//!   ℚ → ℚ(α₁) → ℚ(α₁, α₂) → ... → ℚ(α₁, ..., αₖ)
//!
//! Each extension adjoins αᵢ as a root of a polynomial over the previous field.
//! The simplest case: αᵢ = ⁿ√d for rational d, giving minimal polynomial xⁿ - d.
//!
//! A `RadicalElement` is a vector of coefficients in the power basis:
//!   c₀ + c₁α + c₂α² + ... + c_{d-1}α^{d-1}
//! where d = total_degree of the tower.
//!
//! Arithmetic is polynomial multiplication modulo the tower relations.

use crate::scalar::bigrat::BigRat;
use crate::scalar::coeff::Coeff;
use crate::scalar::rat::Rat;
use std::sync::Arc;

// ═══════════════════════════════════════════════════════════
// TOWER SPECIFICATION
// ═══════════════════════════════════════════════════════════

/// A single layer of the radical extension tower.
///
/// Adjoins a root of `min_poly` to the field built by prior layers.
/// For a simple radical ⁿ√d, min_poly = [−d, 0, ..., 0, 1] (degree n).
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct RadicalLayer {
    /// Degree of this extension.
    pub degree: u32,
    /// Minimal polynomial of the adjoined element over the prior field.
    /// Coefficients are rational (in terms of the base field ℚ).
    /// Length = degree + 1. Leading coefficient is 1 (monic).
    pub min_poly: Vec<Rat>,
}

impl RadicalLayer {
    /// Create a layer for ⁿ√d: minimal polynomial xⁿ - d.
    pub fn nth_root(n: u32, d: Rat) -> Self {
        assert!(n >= 2, "RadicalLayer: degree must be ≥ 2, got {}", n);
        let mut min_poly = vec![Rat::ZERO; n as usize + 1];
        min_poly[0] = -d;
        min_poly[n as usize] = Rat::ONE;
        RadicalLayer {
            degree: n,
            min_poly,
        }
    }

    /// Create a layer for √d: minimal polynomial x² - d.
    pub fn sqrt(d: Rat) -> Self {
        Self::nth_root(2, d)
    }

    /// Create a layer for ∛d: minimal polynomial x³ - d.
    pub fn cbrt(d: Rat) -> Self {
        Self::nth_root(3, d)
    }
}

/// A tower of radical extensions: ℚ → ℚ(α₁) → ℚ(α₁, α₂) → ...
///
/// Shared via `Arc` — multiple RadicalElements can reference the same tower.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct RadicalTower {
    /// Each layer in order. Empty = the rationals ℚ.
    pub layers: Vec<RadicalLayer>,
    /// Product of all layer degrees. Cached for performance.
    pub(crate) total_degree: usize,
}

impl RadicalTower {
    /// The trivial tower: just ℚ (no extensions).
    pub fn rational() -> Self {
        RadicalTower {
            layers: Vec::new(),
            total_degree: 1,
        }
    }

    /// Single-layer tower: ℚ(ⁿ√d).
    pub fn single(layer: RadicalLayer) -> Self {
        let deg = layer.degree as usize;
        RadicalTower {
            layers: vec![layer],
            total_degree: deg,
        }
    }

    /// Tower with two layers: ℚ(α₁, α₂).
    pub fn two(layer1: RadicalLayer, layer2: RadicalLayer) -> Self {
        let deg = layer1.degree as usize * layer2.degree as usize;
        RadicalTower {
            layers: vec![layer1, layer2],
            total_degree: deg,
        }
    }

    /// Construct from an arbitrary list of layers.
    pub fn from_layers(layers: Vec<RadicalLayer>) -> Self {
        let total_degree = layers
            .iter()
            .map(|l| l.degree as usize)
            .product::<usize>()
            .max(1);
        RadicalTower {
            layers,
            total_degree,
        }
    }

    /// Total degree = product of all layer degrees.
    /// This is the dimension of the extension field over ℚ.
    pub fn total_degree(&self) -> usize {
        self.total_degree
    }

    /// Number of extension layers.
    pub fn depth(&self) -> usize {
        self.layers.len()
    }

    /// Is this the trivial tower (just ℚ)?
    pub fn is_rational(&self) -> bool {
        self.layers.is_empty()
    }
}

// ═══════════════════════════════════════════════════════════
// RADICAL ELEMENT
// ═══════════════════════════════════════════════════════════

/// An element of a radical extension field.
///
/// Represented as a vector of rational coefficients in the power basis
/// of the tower. For a single-layer tower ℚ(√d), this is [a, b]
/// representing a + b√d. For a two-layer tower ℚ(√2, √3), this is
/// [a, b, c, d] representing a + b√2 + c√3 + d√6.
#[derive(Clone, Debug)]
pub struct RadicalElement {
    /// Coefficients in the power basis. Length = tower.total_degree().
    /// coeffs[0] is the rational part.
    /// v0.0.2: uses Coeff to support BigRat overflow without panicking.
    pub(crate) coeffs: Vec<Coeff>,
    /// The tower specification (shared).
    pub(crate) tower: Arc<RadicalTower>,
}

impl RadicalElement {
    /// Create from a rational number (trivial tower).
    pub fn from_rat(r: Rat) -> Self {
        RadicalElement {
            coeffs: vec![Coeff::Rat(r)],
            tower: Arc::new(RadicalTower::rational()),
        }
    }

    /// Create from a BigRat (trivial tower, no panic on overflow).
    /// v0.0.2: this is the path that previously panicked.
    pub fn from_bigrat(b: BigRat) -> Self {
        RadicalElement {
            coeffs: vec![Coeff::from_bigrat(b)],
            tower: Arc::new(RadicalTower::rational()),
        }
    }

    /// Create from coefficients and a tower.
    pub fn new(coeffs: Vec<Coeff>, tower: Arc<RadicalTower>) -> Self {
        assert_eq!(
            coeffs.len(),
            tower.total_degree(),
            "RadicalElement: coeffs length {} != tower degree {}",
            coeffs.len(),
            tower.total_degree()
        );
        RadicalElement { coeffs, tower }
    }

    /// Create √d as an element of ℚ(√d).
    pub fn sqrt(d: Rat) -> Self {
        let tower = Arc::new(RadicalTower::single(RadicalLayer::sqrt(d)));
        RadicalElement {
            coeffs: vec![Coeff::ZERO, Coeff::ONE], // 0 + 1·√d
            tower,
        }
    }

    /// Create ∛d as an element of ℚ(∛d).
    pub fn cbrt(d: Rat) -> Self {
        let tower = Arc::new(RadicalTower::single(RadicalLayer::cbrt(d)));
        RadicalElement {
            coeffs: vec![Coeff::ZERO, Coeff::ONE, Coeff::ZERO], // 0 + 1·∛d + 0·(∛d)²
            tower,
        }
    }

    /// Is this a purely rational element?
    pub fn is_rational(&self) -> bool {
        if self.tower.is_rational() {
            return true;
        }
        self.coeffs[1..].iter().all(|c| c.is_zero())
    }

    /// Extract as Rat if purely rational.
    pub fn to_rat(&self) -> Option<Rat> {
        if self.is_rational() {
            self.coeffs[0].to_rat()
        } else {
            None
        }
    }

    /// Is this element zero?
    pub fn is_zero(&self) -> bool {
        self.coeffs.iter().all(|c| c.is_zero())
    }

    /// Is this element strictly positive?
    ///
    /// For single sqrt tower ℚ(√d): a + b√d > 0 is computed exactly.
    /// For rational elements: trivial.
    /// For multi-layer: uses conservative rational approximation.
    pub fn is_positive(&self) -> bool {
        if self.is_zero() {
            return false;
        }
        if self.tower.is_rational() {
            return self.coeffs[0].is_positive();
        }
        if self.tower.depth() == 1 && self.tower.layers[0].degree == 2 {
            return self.sign_single_sqrt() > 0;
        }
        // Multi-layer: approximate by evaluating with rational lower bounds for each radical
        // Conservative: check if rational part dominates
        self.approximate_sign() > 0
    }

    /// Is this element strictly negative?
    pub fn is_negative(&self) -> bool {
        if self.is_zero() {
            return false;
        }
        if self.tower.is_rational() {
            return self.coeffs[0].is_negative();
        }
        if self.tower.depth() == 1 && self.tower.layers[0].degree == 2 {
            return self.sign_single_sqrt() < 0;
        }
        self.approximate_sign() < 0
    }

    /// Exact sign for a + b√d: returns -1, 0, or 1.
    fn sign_single_sqrt(&self) -> i8 {
        let a = self.coeffs[0].clone();
        let b = self.coeffs[1].clone();
        let d = Coeff::Rat(-self.tower.layers[0].min_poly[0]); // x² - d, so d = -min_poly[0]

        if b.is_zero() {
            return if a.is_positive() {
                1
            } else if a.is_negative() {
                -1
            } else {
                0
            };
        }
        if a.is_zero() {
            return if (b.is_positive() && d.is_positive()) || (b.is_negative() && d.is_negative()) {
                1
            } else {
                -1
            };
        }

        if d.is_positive() {
            if a.is_positive() && b.is_positive() {
                return 1;
            }
            if a.is_negative() && b.is_negative() {
                return -1;
            }

            let a2 = a.clone() * a.clone();
            let b2d = b.clone() * b.clone() * d;
            if a.is_positive() {
                if a2 > b2d {
                    1
                } else if a2 < b2d {
                    -1
                } else {
                    0
                }
            } else if b2d > a2 {
                1
            } else if b2d < a2 {
                -1
            } else {
                0
            }
        } else if a.is_positive() {
            1
        } else {
            -1
        }
    }

    /// Approximate sign for multi-layer towers.
    fn approximate_sign(&self) -> i8 {
        if self.tower.depth() == 2 {
            let d0 = -self.tower.layers[0].min_poly[0];
            let d1 = -self.tower.layers[1].min_poly[0];
            let scale = Rat::from(10000i64);
            let approx_sqrt = |d: Rat| -> Coeff {
                if d.is_zero() || d.is_negative() {
                    return Coeff::ZERO;
                }
                let n = (d * scale * scale).num() / (d * scale * scale).den() as i128;
                if n <= 0 {
                    return Coeff::ZERO;
                }
                let mut s = 1i128;
                while s * s <= n {
                    s += 1;
                }
                s -= 1;
                Coeff::Rat(Rat::new(s, 10000))
            };

            let s0 = approx_sqrt(d0);
            let s1 = approx_sqrt(d1);
            let s01 = approx_sqrt(d0 * d1);

            let val = self.coeffs[0].clone()
                + self.coeffs[1].clone() * s0
                + self.coeffs[2].clone() * s1
                + self.coeffs[3].clone() * s01;
            if val.is_positive() {
                return 1;
            }
            if val.is_negative() {
                return -1;
            }
        }

        if self.coeffs[0].is_positive() {
            1
        } else if self.coeffs[0].is_negative() {
            -1
        } else {
            0
        }
    }

    /// The tower this element belongs to.
    pub fn tower(&self) -> &RadicalTower {
        &self.tower
    }

    /// Negate.
    pub fn neg(&self) -> Self {
        RadicalElement {
            coeffs: self.coeffs.iter().map(|c| -c.clone()).collect(),
            tower: Arc::clone(&self.tower),
        }
    }

    // ─── Same-tower arithmetic ───

    /// Add two elements in the same tower.
    fn add_same_tower(&self, other: &Self) -> Self {
        debug_assert!(Arc::ptr_eq(&self.tower, &other.tower) || *self.tower == *other.tower);
        let coeffs = self
            .coeffs
            .iter()
            .zip(&other.coeffs)
            .map(|(a, b)| a.clone() + b.clone())
            .collect();
        RadicalElement {
            coeffs,
            tower: Arc::clone(&self.tower),
        }
    }

    /// Subtract two elements in the same tower.
    fn sub_same_tower(&self, other: &Self) -> Self {
        debug_assert!(Arc::ptr_eq(&self.tower, &other.tower) || *self.tower == *other.tower);
        let coeffs = self
            .coeffs
            .iter()
            .zip(&other.coeffs)
            .map(|(a, b)| a.clone() - b.clone())
            .collect();
        RadicalElement {
            coeffs,
            tower: Arc::clone(&self.tower),
        }
    }

    /// Multiply two elements in the same tower.
    ///
    /// Polynomial multiplication modulo the tower relations.
    /// For a single-layer tower ℚ(α) with α^n = d:
    ///   (a₀ + a₁α + ... + a_{n-1}α^{n-1}) * (b₀ + b₁α + ... + b_{n-1}α^{n-1})
    /// = polynomial product, then reduce modulo (α^n - d).
    fn mul_same_tower(&self, other: &Self) -> Self {
        let d = self.tower.total_degree();
        if d == 1 {
            // Trivial tower: just coefficient multiplication
            let c = self.coeffs[0].clone() * other.coeffs[0].clone();
            return RadicalElement {
                coeffs: vec![c],
                tower: Arc::clone(&self.tower),
            };
        }

        if self.tower.depth() == 1 {
            // Single-layer tower: polynomial multiplication mod min_poly
            return self.mul_single_layer(other);
        }

        // Multi-layer tower: use the recursive Kronecker product structure.
        // Each element is indexed by a multi-index (i₁, i₂, ..., iₖ) where
        // 0 ≤ iⱼ < degree_j. Multiplication is polynomial multiplication
        // in each variable, reduced modulo the corresponding min_poly.
        self.mul_multi_layer(other)
    }

    /// Multiplication in a single-layer tower: polynomial mod min_poly.
    fn mul_single_layer(&self, other: &Self) -> Self {
        let n = self.tower.layers[0].degree as usize;
        let min_poly = &self.tower.layers[0].min_poly;

        // Step 1: polynomial multiplication (degree up to 2(n-1))
        let mut product = vec![Coeff::ZERO; 2 * n - 1];
        for (i, ai) in self.coeffs.iter().enumerate() {
            if ai.is_zero() {
                continue;
            }
            for (j, bj) in other.coeffs.iter().enumerate() {
                if bj.is_zero() {
                    continue;
                }
                product[i + j] = product[i + j].clone() + ai.clone() * bj.clone();
            }
        }

        // Step 2: reduce modulo min_poly = x^n - d (where d = -min_poly[0])
        if min_poly.len() == n + 1 && min_poly[n] == Rat::ONE {
            for k in (n..product.len()).rev() {
                let overflow = product[k].clone();
                if overflow.is_zero() {
                    continue;
                }
                product[k] = Coeff::ZERO;
                for i in 0..n {
                    product[k - n + i] =
                        product[k - n + i].clone() - overflow.clone() * Coeff::Rat(min_poly[i]);
                }
            }
        }

        let coeffs = product[..n].to_vec();
        RadicalElement {
            coeffs,
            tower: Arc::clone(&self.tower),
        }
    }

    /// Multiplication in a multi-layer tower using tensor product structure.
    ///
    /// Each element is indexed by a multi-index (i₁, ..., iₖ) where 0 ≤ iⱼ < dⱼ.
    /// Multiplication is polynomial multiplication in each variable,
    /// then reduction modulo each min_poly from last dimension to first.
    fn mul_multi_layer(&self, other: &Self) -> Self {
        let dims: Vec<usize> = self
            .tower
            .layers
            .iter()
            .map(|l| l.degree as usize)
            .collect();
        let k = dims.len();

        // Extended dimensions for unreduced product
        let ext_dims: Vec<usize> = dims.iter().map(|&d| 2 * d - 1).collect();
        let ext_total: usize = ext_dims.iter().product();
        let mut product = vec![Coeff::ZERO; ext_total];

        // Step 1: polynomial multiplication via multi-indices
        for (i, ai) in self.coeffs.iter().enumerate() {
            if ai.is_zero() {
                continue;
            }
            let mi = flat_to_multi(i, &dims);
            for (j, bj) in other.coeffs.iter().enumerate() {
                if bj.is_zero() {
                    continue;
                }
                let mj = flat_to_multi(j, &dims);
                let mk: Vec<usize> = mi.iter().zip(&mj).map(|(a, b)| a + b).collect();
                if let Some(flat) = multi_to_flat(&mk, &ext_dims) {
                    product[flat] = product[flat].clone() + ai.clone() * bj.clone();
                }
            }
        }

        // Step 2: reduce each dimension modulo its min_poly, from last to first
        let mut cur_dims = ext_dims;
        for dim_idx in (0..k).rev() {
            let n = dims[dim_idx];
            let cur_n = cur_dims[dim_idx];
            if cur_n <= n {
                continue;
            }

            let min_poly = &self.tower.layers[dim_idx].min_poly;

            // Strides for iterating fibers along this dimension
            let inner_stride: usize = cur_dims[dim_idx + 1..].iter().product();
            let outer_count: usize = cur_dims[..dim_idx].iter().product();
            let dim_stride = cur_n * inner_stride;

            // New array with this dimension reduced
            let new_dim_stride = n * inner_stride;
            let new_total = outer_count * new_dim_stride;
            let mut reduced = vec![Coeff::ZERO; new_total];

            for outer in 0..outer_count {
                for inner in 0..inner_stride {
                    let mut slice = Vec::with_capacity(cur_n);
                    for e in 0..cur_n {
                        slice.push(product[outer * dim_stride + e * inner_stride + inner].clone());
                    }

                    for e in (n..cur_n).rev() {
                        let c = slice[e].clone();
                        if c.is_zero() {
                            continue;
                        }
                        slice[e] = Coeff::ZERO;
                        for i in 0..n {
                            slice[e - n + i] =
                                slice[e - n + i].clone() - c.clone() * Coeff::Rat(min_poly[i]);
                        }
                    }

                    for e in 0..n {
                        reduced[outer * new_dim_stride + e * inner_stride + inner] =
                            slice[e].clone();
                    }
                }
            }

            product = reduced;
            cur_dims[dim_idx] = n;
        }

        product.truncate(self.tower.total_degree());
        RadicalElement {
            coeffs: product,
            tower: Arc::clone(&self.tower),
        }
    }

    /// Division: self / other in the same tower.
    /// Computed as self * other.recip().
    fn div_same_tower(&self, other: &Self) -> Self {
        let inv = other.recip_same_tower();
        self.mul_same_tower(&inv)
    }

    /// Reciprocal in the same tower.
    /// For single-layer: use extended GCD in the polynomial ring.
    /// For multi-layer: use the norm-based approach.
    fn recip_same_tower(&self) -> Self {
        assert!(!self.is_zero(), "RadicalElement: reciprocal of zero");

        if self.tower.is_rational() {
            return RadicalElement {
                coeffs: vec![self.coeffs[0].clone().recip()],
                tower: Arc::clone(&self.tower),
            };
        }

        if self.tower.depth() == 1 {
            // Try Rat-path first (most common). If coeffs are all Rat, use extended GCD.
            let rat_coeffs: Option<Vec<Rat>> = self.coeffs.iter().map(|c| c.to_rat()).collect();
            if let Some(rc) = rat_coeffs {
                return self.recip_single_layer_rat(&rc);
            }
        }

        // General path: matrix inversion with Coeff arithmetic
        self.recip_multi_layer()
    }

    /// Reciprocal in a single-layer tower via extended GCD (Rat fast path).
    fn recip_single_layer_rat(&self, rat_coeffs: &[Rat]) -> Self {
        let n = self.tower.layers[0].degree as usize;
        let min_poly = &self.tower.layers[0].min_poly;

        let (gcd, s, _) = poly_extended_gcd(rat_coeffs, min_poly);

        let gcd_val = gcd
            .iter()
            .find(|c| !c.is_zero())
            .expect("RadicalElement: element has no inverse (GCD is zero)");
        let inv_gcd = gcd_val.recip();

        let mut inv: Vec<Coeff> = s.iter().map(|c| Coeff::Rat(*c * inv_gcd)).collect();
        inv.truncate(n);
        while inv.len() < n {
            inv.push(Coeff::ZERO);
        }

        RadicalElement {
            coeffs: inv,
            tower: Arc::clone(&self.tower),
        }
    }

    /// Reciprocal via matrix inversion with Coeff arithmetic.
    #[allow(clippy::needless_range_loop)]
    fn recip_multi_layer(&self) -> Self {
        let d = self.tower.total_degree();

        let mut matrix = vec![vec![Coeff::ZERO; d + 1]; d];
        for j in 0..d {
            let mut ej = vec![Coeff::ZERO; d];
            ej[j] = Coeff::ONE;
            let ej_elem = RadicalElement::new(ej, Arc::clone(&self.tower));
            let product = self.mul_same_tower(&ej_elem);
            for i in 0..d {
                matrix[i][j] = product.coeffs[i].clone();
            }
        }

        matrix[0][d] = Coeff::ONE;

        for col in 0..d {
            let pivot = (col..d).find(|&r| !matrix[r][col].is_zero());
            let pivot = pivot.expect("RadicalElement: element not invertible in tower");
            matrix.swap(col, pivot);

            let inv_pivot = matrix[col][col].clone().recip();
            for j in 0..=d {
                matrix[col][j] = matrix[col][j].clone() * inv_pivot.clone();
            }

            for row in 0..d {
                if row == col {
                    continue;
                }
                let factor = matrix[row][col].clone();
                if factor.is_zero() {
                    continue;
                }
                for j in 0..=d {
                    let val = matrix[col][j].clone();
                    matrix[row][j] = matrix[row][j].clone() - factor.clone() * val;
                }
            }
        }

        let coeffs: Vec<Coeff> = (0..d).map(|i| matrix[i][d].clone()).collect();
        RadicalElement {
            coeffs,
            tower: Arc::clone(&self.tower),
        }
    }

    // ─── Cross-tower arithmetic (lazy merge) ───

    /// Add two elements, potentially from different towers.
    pub fn add(&self, other: &Self) -> Self {
        if Arc::ptr_eq(&self.tower, &other.tower) || *self.tower == *other.tower {
            return self.add_same_tower(other);
        }
        // Different towers: embed both into the combined tower
        let (a, b) = embed_into_common_tower(self, other);
        a.add_same_tower(&b)
    }

    /// Subtract two elements, potentially from different towers.
    pub fn sub(&self, other: &Self) -> Self {
        if Arc::ptr_eq(&self.tower, &other.tower) || *self.tower == *other.tower {
            return self.sub_same_tower(other);
        }
        let (a, b) = embed_into_common_tower(self, other);
        a.sub_same_tower(&b)
    }

    /// Multiply two elements, potentially from different towers.
    pub fn mul(&self, other: &Self) -> Self {
        if Arc::ptr_eq(&self.tower, &other.tower) || *self.tower == *other.tower {
            return self.mul_same_tower(other);
        }
        let (a, b) = embed_into_common_tower(self, other);
        a.mul_same_tower(&b)
    }

    /// Divide two elements, potentially from different towers.
    pub fn div(&self, other: &Self) -> Self {
        if Arc::ptr_eq(&self.tower, &other.tower) || *self.tower == *other.tower {
            return self.div_same_tower(other);
        }
        let (a, b) = embed_into_common_tower(self, other);
        a.div_same_tower(&b)
    }
}

// ═══════════════════════════════════════════════════════════
// TOWER MERGING
// ═══════════════════════════════════════════════════════════

/// Embed two elements from different towers into a common tower.
///
/// The combined tower contains all layers from both towers.
/// Duplicate layers (same min_poly) are shared, not duplicated.
fn embed_into_common_tower(
    a: &RadicalElement,
    b: &RadicalElement,
) -> (RadicalElement, RadicalElement) {
    // If either is rational, embed into the other's tower
    if a.tower.is_rational() {
        let embedded_a = embed_rational_into(a.coeffs[0].clone(), &b.tower);
        return (embedded_a, b.clone());
    }
    if b.tower.is_rational() {
        let embedded_b = embed_rational_into(b.coeffs[0].clone(), &a.tower);
        return (a.clone(), embedded_b);
    }

    // Build combined tower: layers from a's tower + new layers from b's tower
    let mut combined_layers = a.tower.layers.clone();
    let mut b_layer_mapping: Vec<Option<usize>> = Vec::new();

    for b_layer in &b.tower.layers {
        // Check if this layer already exists in the combined tower
        let existing = combined_layers.iter().position(|l| l == b_layer);
        match existing {
            Some(idx) => b_layer_mapping.push(Some(idx)),
            None => {
                b_layer_mapping.push(Some(combined_layers.len()));
                combined_layers.push(b_layer.clone());
            }
        }
    }

    let combined_degree: usize = combined_layers.iter().map(|l| l.degree as usize).product();
    let combined = Arc::new(RadicalTower {
        layers: combined_layers,
        total_degree: combined_degree,
    });

    // Embed a: its layers are the first layers of the combined tower.
    // Coefficients are embedded by placing them in the sub-lattice
    // corresponding to the first tower's dimensions.
    let a_embedded = embed_element_into_combined(a, &combined, 0..a.tower.depth());
    let b_embedded = embed_element_into_combined_mapped(b, &combined, &b_layer_mapping);

    (a_embedded, b_embedded)
}

/// Embed a rational into a tower as the constant element.
fn embed_rational_into(r: Coeff, tower: &Arc<RadicalTower>) -> RadicalElement {
    let mut coeffs = vec![Coeff::ZERO; tower.total_degree()];
    coeffs[0] = r;
    RadicalElement {
        coeffs,
        tower: Arc::clone(tower),
    }
}

/// Embed an element whose tower layers are a prefix of the combined tower.
fn embed_element_into_combined(
    elem: &RadicalElement,
    combined: &Arc<RadicalTower>,
    _layer_range: std::ops::Range<usize>,
) -> RadicalElement {
    let src_degree = elem.tower.total_degree();
    let dst_degree = combined.total_degree();

    if src_degree == dst_degree {
        return RadicalElement {
            coeffs: elem.coeffs.clone(),
            tower: Arc::clone(combined),
        };
    }

    // The source element lives in the first `src_degree` dimensions.
    // In the combined tower, these correspond to indices 0, 1, ..., src_degree-1
    // if the source layers are a prefix. The remaining dimensions get zero.
    let stride = dst_degree / src_degree;
    let mut coeffs = vec![Coeff::ZERO; dst_degree];
    for (i, c) in elem.coeffs.iter().enumerate() {
        coeffs[i * stride] = c.clone();
    }

    RadicalElement {
        coeffs,
        tower: Arc::clone(combined),
    }
}

/// Embed an element whose tower layers map into the combined tower via mapping.
fn embed_element_into_combined_mapped(
    elem: &RadicalElement,
    combined: &Arc<RadicalTower>,
    _mapping: &[Option<usize>],
) -> RadicalElement {
    let dst_degree = combined.total_degree();

    if elem.tower.is_rational() {
        return embed_rational_into(elem.coeffs[0].clone(), combined);
    }

    let src_dims: Vec<usize> = elem
        .tower
        .layers
        .iter()
        .map(|l| l.degree as usize)
        .collect();
    let dst_dims: Vec<usize> = combined.layers.iter().map(|l| l.degree as usize).collect();

    let mut coeffs = vec![Coeff::ZERO; dst_degree];

    for (src_flat, c) in elem.coeffs.iter().enumerate() {
        if c.is_zero() {
            continue;
        }
        let src_multi = flat_to_multi(src_flat, &src_dims);
        let mut dst_multi = vec![0usize; dst_dims.len()];
        let offset = dst_dims.len() - src_dims.len();
        for (j, &idx) in src_multi.iter().enumerate() {
            dst_multi[offset + j] = idx;
        }
        if let Some(dst_flat) = multi_to_flat(&dst_multi, &dst_dims) {
            coeffs[dst_flat] = c.clone();
        }
    }

    RadicalElement {
        coeffs,
        tower: Arc::clone(combined),
    }
}

// ═══════════════════════════════════════════════════════════
// MULTI-INDEX HELPERS
// ═══════════════════════════════════════════════════════════

/// Convert flat index to multi-index given dimensions.
fn flat_to_multi(mut flat: usize, dims: &[usize]) -> Vec<usize> {
    let mut result = vec![0usize; dims.len()];
    for i in (0..dims.len()).rev() {
        result[i] = flat % dims[i];
        flat /= dims[i];
    }
    result
}

/// Convert multi-index to flat index. Returns None if any index is out of range.
fn multi_to_flat(multi: &[usize], dims: &[usize]) -> Option<usize> {
    let mut flat = 0;
    let mut stride = 1;
    for i in (0..dims.len()).rev() {
        if multi[i] >= dims[i] {
            return None;
        }
        flat += multi[i] * stride;
        stride *= dims[i];
    }
    Some(flat)
}

// ═══════════════════════════════════════════════════════════
// POLYNOMIAL EXTENDED GCD (for reciprocal computation)
// ═══════════════════════════════════════════════════════════

/// Extended GCD for univariate polynomials over ℚ.
/// Returns (gcd, s, t) such that gcd = s*a + t*b.
fn poly_extended_gcd(a: &[Rat], b: &[Rat]) -> (Vec<Rat>, Vec<Rat>, Vec<Rat>) {
    let mut old_r = a.to_vec();
    let mut r = b.to_vec();
    let mut old_s = vec![Rat::ONE];
    let mut s = vec![Rat::ZERO];
    let mut old_t = vec![Rat::ZERO];
    let mut t = vec![Rat::ONE];

    poly_trim_rats(&mut old_r);
    poly_trim_rats(&mut r);

    while !r.iter().all(|c| c.is_zero()) {
        let (q, rem) = poly_div_rats(&old_r, &r);

        old_r = r;
        r = rem;

        let new_s = poly_sub_rats(&old_s, &poly_mul_rats(&q, &s));
        old_s = s;
        s = new_s;

        let new_t = poly_sub_rats(&old_t, &poly_mul_rats(&q, &t));
        old_t = t;
        t = new_t;
    }

    (old_r, old_s, old_t)
}

fn poly_trim_rats(p: &mut Vec<Rat>) {
    while p.len() > 1 && p.last().is_some_and(|c| c.is_zero()) {
        p.pop();
    }
}

fn poly_mul_rats(a: &[Rat], b: &[Rat]) -> Vec<Rat> {
    if a.is_empty() || b.is_empty() {
        return vec![Rat::ZERO];
    }
    let mut r = vec![Rat::ZERO; a.len() + b.len() - 1];
    for (i, &ca) in a.iter().enumerate() {
        for (j, &cb) in b.iter().enumerate() {
            r[i + j] += ca * cb;
        }
    }
    poly_trim_rats(&mut r);
    r
}

fn poly_sub_rats(a: &[Rat], b: &[Rat]) -> Vec<Rat> {
    let n = a.len().max(b.len());
    let mut r = vec![Rat::ZERO; n];
    for (i, v) in r.iter_mut().enumerate() {
        let ai = if i < a.len() { a[i] } else { Rat::ZERO };
        let bi = if i < b.len() { b[i] } else { Rat::ZERO };
        *v = ai - bi;
    }
    poly_trim_rats(&mut r);
    r
}

fn poly_div_rats(a: &[Rat], b: &[Rat]) -> (Vec<Rat>, Vec<Rat>) {
    let da = a.len().saturating_sub(1);
    let db = b.len().saturating_sub(1);
    if da < db || b.iter().all(|c| c.is_zero()) {
        return (vec![Rat::ZERO], a.to_vec());
    }
    let mut rem = a.to_vec();
    let mut quot = vec![Rat::ZERO; da - db + 1];
    let lead_b = *b.last().unwrap();
    for i in (0..=da - db).rev() {
        let coeff = rem[i + db] / lead_b;
        quot[i] = coeff;
        for j in 0..=db {
            rem[i + j] -= coeff * b[j];
        }
    }
    poly_trim_rats(&mut rem);
    poly_trim_rats(&mut quot);
    (quot, rem)
}

// ═══════════════════════════════════════════════════════════
// TRAIT IMPLEMENTATIONS
// ═══════════════════════════════════════════════════════════

impl PartialEq for RadicalElement {
    fn eq(&self, other: &Self) -> bool {
        if Arc::ptr_eq(&self.tower, &other.tower) || *self.tower == *other.tower {
            return self.coeffs == other.coeffs;
        }
        // Cross-tower: subtract and check zero
        let diff = self.sub(other);
        diff.is_zero()
    }
}

impl Eq for RadicalElement {}

impl std::fmt::Display for RadicalElement {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        if self.is_rational() {
            return write!(f, "{}", self.coeffs[0]);
        }
        if self.tower.depth() == 1 && self.tower.layers[0].degree == 2 {
            let a = &self.coeffs[0];
            let b = &self.coeffs[1];
            let d = -self.tower.layers[0].min_poly[0];
            if a.is_zero() {
                write!(f, "{}√{}", b, d)
            } else if b.is_zero() {
                write!(f, "{}", a)
            } else if b.is_positive() {
                write!(f, "{} + {}√{}", a, b, d)
            } else {
                write!(f, "{} - {}√{}", a, b.clone().abs(), d)
            }
        } else {
            write!(
                f,
                "Radical(deg={}, coeffs={:?})",
                self.tower.total_degree(),
                self.coeffs
            )
        }
    }
}

// ═══════════════════════════════════════════════════════════
// TESTS
// ═══════════════════════════════════════════════════════════

#[cfg(test)]
mod tests {
    use super::*;

    // ─── Basic construction ───

    #[test]
    fn from_rat_roundtrip() {
        let r = RadicalElement::from_rat(Rat::new(3, 7));
        assert!(r.is_rational());
        assert_eq!(r.to_rat(), Some(Rat::new(3, 7)));
    }

    #[test]
    fn sqrt2_is_not_rational() {
        let s = RadicalElement::sqrt(Rat::from(2));
        assert!(!s.is_rational());
        assert_eq!(s.to_rat(), None);
    }

    #[test]
    fn cbrt2_construction() {
        let c = RadicalElement::cbrt(Rat::from(2));
        assert!(!c.is_rational());
        assert_eq!(c.tower.total_degree(), 3);
        assert_eq!(c.coeffs, vec![Coeff::ZERO, Coeff::ONE, Coeff::ZERO]);
    }

    // ─── Same-tower arithmetic ───

    #[test]
    fn rat_add() {
        let a = RadicalElement::from_rat(Rat::from(3));
        let b = RadicalElement::from_rat(Rat::from(4));
        let c = a.add(&b);
        assert_eq!(c.to_rat(), Some(Rat::from(7)));
    }

    #[test]
    fn rat_mul() {
        let a = RadicalElement::from_rat(Rat::from(3));
        let b = RadicalElement::from_rat(Rat::from(5));
        let c = a.mul(&b);
        assert_eq!(c.to_rat(), Some(Rat::from(15)));
    }

    #[test]
    fn sqrt2_add_sqrt2() {
        let s = RadicalElement::sqrt(Rat::from(2));
        let two_s = s.add(&s);
        // √2 + √2 = 2√2 → coeffs [0, 2]
        assert_eq!(two_s.coeffs[0], Coeff::ZERO);
        assert_eq!(two_s.coeffs[1], Coeff::Rat(Rat::from(2)));
    }

    #[test]
    fn sqrt2_times_sqrt2() {
        // √2 * √2 = 2
        let s = RadicalElement::sqrt(Rat::from(2));
        let product = s.mul(&s);
        assert!(product.is_rational());
        assert_eq!(product.to_rat(), Some(Rat::from(2)));
    }

    #[test]
    fn one_plus_sqrt2_times_one_minus_sqrt2() {
        // (1 + √2)(1 - √2) = 1 - 2 = -1
        let tower = Arc::new(RadicalTower::single(RadicalLayer::sqrt(Rat::from(2))));
        let a = RadicalElement::new(vec![Coeff::ONE, Coeff::ONE], Arc::clone(&tower));
        let b = RadicalElement::new(
            vec![Coeff::ONE, Coeff::Rat(Rat::NEG_ONE)],
            Arc::clone(&tower),
        );
        let c = a.mul(&b);
        assert!(c.is_rational());
        assert_eq!(c.to_rat(), Some(Rat::from(-1)));
    }

    #[test]
    fn sqrt2_recip() {
        // 1/√2 = √2/2 → coeffs [0, 1/2]
        let s = RadicalElement::sqrt(Rat::from(2));
        let inv = s.recip_same_tower();
        assert_eq!(inv.coeffs[0], Coeff::ZERO);
        assert_eq!(inv.coeffs[1], Coeff::Rat(Rat::new(1, 2)));
        // Verify: √2 * (√2/2) = 1
        let product = s.mul(&inv);
        assert!(product.is_rational());
        assert_eq!(product.to_rat(), Some(Rat::ONE));
    }

    #[test]
    fn one_plus_sqrt2_recip() {
        // 1/(1 + √2) = -1 + √2 (by conjugate: (1-√2)/(1-2) = -(1-√2) = -1+√2)
        let tower = Arc::new(RadicalTower::single(RadicalLayer::sqrt(Rat::from(2))));
        let a = RadicalElement::new(vec![Coeff::ONE, Coeff::ONE], Arc::clone(&tower));
        let inv = a.recip_same_tower();
        assert_eq!(inv.coeffs[0], Coeff::Rat(Rat::from(-1)));
        assert_eq!(inv.coeffs[1], Coeff::ONE);
        // Verify
        let product = a.mul(&inv);
        assert!(product.is_rational());
        assert_eq!(product.to_rat(), Some(Rat::ONE));
    }

    #[test]
    fn cbrt2_cubed() {
        // (∛2)³ = 2
        let c = RadicalElement::cbrt(Rat::from(2));
        let c2 = c.mul(&c);
        let c3 = c2.mul(&c);
        assert!(c3.is_rational());
        assert_eq!(c3.to_rat(), Some(Rat::from(2)));
    }

    // ─── Cross-tower arithmetic (the √2 + √3 test!) ───

    #[test]
    fn sqrt2_plus_sqrt3_no_panic() {
        // This was the original QuadraticSurd panic. Now it works.
        let s2 = RadicalElement::sqrt(Rat::from(2));
        let s3 = RadicalElement::sqrt(Rat::from(3));
        let sum = s2.add(&s3);
        // Should be in ℚ(√2, √3), degree 4
        assert_eq!(sum.tower.total_degree(), 4);
        assert!(!sum.is_rational());
        assert!(!sum.is_zero());
    }

    #[test]
    fn sqrt2_times_sqrt3() {
        let s2 = RadicalElement::sqrt(Rat::from(2));
        let s3 = RadicalElement::sqrt(Rat::from(3));
        let product = s2.mul(&s3);
        // √2 * √3 = √6, should be in ℚ(√2, √3)
        assert_eq!(product.tower.total_degree(), 4);
        assert!(!product.is_rational());
    }

    #[test]
    fn sqrt2_plus_sqrt3_squared() {
        // (√2 + √3)² = 2 + 2√6 + 3 = 5 + 2√6
        let s2 = RadicalElement::sqrt(Rat::from(2));
        let s3 = RadicalElement::sqrt(Rat::from(3));
        let sum = s2.add(&s3);
        let sq = sum.mul(&sum);
        // Should not be rational (has √6 component)
        assert!(!sq.is_rational());
        // But the rational part should be 5
        assert_eq!(sq.coeffs[0], Coeff::Rat(Rat::from(5)));
    }

    #[test]
    fn rat_plus_sqrt2() {
        // 3 + √2
        let three = RadicalElement::from_rat(Rat::from(3));
        let s2 = RadicalElement::sqrt(Rat::from(2));
        let sum = three.add(&s2);
        assert_eq!(sum.tower.total_degree(), 2);
        assert_eq!(sum.coeffs[0], Coeff::Rat(Rat::from(3)));
        assert_eq!(sum.coeffs[1], Coeff::ONE);
    }

    // ─── Equality ───

    #[test]
    fn same_tower_equality() {
        let tower = Arc::new(RadicalTower::single(RadicalLayer::sqrt(Rat::from(2))));
        let a = RadicalElement::new(
            vec![Coeff::Rat(Rat::from(3)), Coeff::Rat(Rat::from(4))],
            Arc::clone(&tower),
        );
        let b = RadicalElement::new(
            vec![Coeff::Rat(Rat::from(3)), Coeff::Rat(Rat::from(4))],
            Arc::clone(&tower),
        );
        assert_eq!(a, b);
    }

    #[test]
    fn cross_type_rat_equality() {
        let a = RadicalElement::from_rat(Rat::from(3));
        let tower = Arc::new(RadicalTower::single(RadicalLayer::sqrt(Rat::from(2))));
        let b = RadicalElement::new(
            vec![Coeff::Rat(Rat::from(3)), Coeff::ZERO],
            Arc::clone(&tower),
        );
        assert_eq!(a, b);
    }

    // ─── Display ───

    #[test]
    fn display_rational() {
        let r = RadicalElement::from_rat(Rat::new(3, 7));
        assert_eq!(format!("{}", r), "3/7");
    }

    #[test]
    fn display_sqrt() {
        let s = RadicalElement::sqrt(Rat::from(2));
        assert_eq!(format!("{}", s), "1√2"); // 0 + 1·√2, a=0 so shows "b√d"
    }
}