geoit 0.0.2

//! Arbitrary-precision integer arithmetic.
//!
//! `BigUint` is the unsigned core: `Vec<u64>` limbs, little-endian, no leading zeros.
//! `BigInt` is the signed wrapper: sign + magnitude.
//!
//! Zero dependencies. Every line is ours.
#![allow(clippy::needless_range_loop)]

use std::cmp::Ordering;

// ═══════════════════════════════════════════════════════════
// BIGUINT
// ═══════════════════════════════════════════════════════════

/// Unsigned arbitrary-precision integer.
///
/// Representation: `limbs[0]` is the least significant u64.
/// Empty `limbs` means zero. No leading zero limbs (invariant maintained by `normalize`).
#[derive(Clone, Debug)]
pub struct BigUint {
    limbs: Vec<u64>,
}

impl BigUint {
    // ─── Constructors ───

    pub const ZERO: BigUint = BigUint { limbs: Vec::new() };

    pub fn zero() -> Self {
        BigUint { limbs: Vec::new() }
    }

    pub fn one() -> Self {
        BigUint { limbs: vec![1] }
    }

    pub fn from_u64(v: u64) -> Self {
        if v == 0 {
            Self::zero()
        } else {
            BigUint { limbs: vec![v] }
        }
    }

    pub fn from_u128(v: u128) -> Self {
        if v == 0 {
            Self::zero()
        } else {
            let lo = v as u64;
            let hi = (v >> 64) as u64;
            let mut limbs = vec![lo];
            if hi != 0 {
                limbs.push(hi);
            }
            BigUint { limbs }
        }
    }

    /// Construct from little-endian u64 limbs. Normalizes.
    pub fn from_limbs(limbs: Vec<u64>) -> Self {
        let mut b = BigUint { limbs };
        b.normalize();
        b
    }

    // ─── Accessors ───

    pub fn is_zero(&self) -> bool {
        self.limbs.is_empty()
    }

    pub fn is_one(&self) -> bool {
        self.limbs.len() == 1 && self.limbs[0] == 1
    }

    /// Number of significant limbs.
    pub fn len(&self) -> usize {
        self.limbs.len()
    }

    /// Whether this BigUint has no limbs (should not occur after normalization).
    pub fn is_empty(&self) -> bool {
        self.limbs.is_empty()
    }

    /// The limbs, little-endian.
    pub fn limbs(&self) -> &[u64] {
        &self.limbs
    }

    /// Try to convert to u128. Returns None if value doesn't fit.
    pub fn to_u128(&self) -> Option<u128> {
        match self.limbs.len() {
            0 => Some(0),
            1 => Some(self.limbs[0] as u128),
            2 => Some(self.limbs[0] as u128 | ((self.limbs[1] as u128) << 64)),
            _ => None,
        }
    }

    /// Number of significant bits.
    pub fn bit_len(&self) -> u64 {
        if self.is_zero() {
            return 0;
        }
        let top = self.limbs.len() - 1;
        (top as u64) * 64 + (64 - self.limbs[top].leading_zeros() as u64)
    }

    /// Get bit at position `i` (0-indexed from LSB).
    pub fn bit(&self, i: u64) -> bool {
        let limb_idx = (i / 64) as usize;
        let bit_idx = i % 64;
        if limb_idx >= self.limbs.len() {
            return false;
        }
        (self.limbs[limb_idx] >> bit_idx) & 1 == 1
    }

    /// True if the value is even (or zero).
    pub fn is_even(&self) -> bool {
        self.is_zero() || (self.limbs[0] & 1 == 0)
    }

    // ─── Normalization ───

    fn normalize(&mut self) {
        while self.limbs.last() == Some(&0) {
            self.limbs.pop();
        }
    }

    // ─── Addition ───

    /// self += rhs
    pub fn add_assign(&mut self, rhs: &BigUint) {
        let n = self.limbs.len().max(rhs.limbs.len());
        self.limbs.resize(n, 0);
        let mut carry = 0u64;
        for i in 0..n {
            let a = self.limbs[i];
            let b = if i < rhs.limbs.len() { rhs.limbs[i] } else { 0 };
            let (sum1, c1) = a.overflowing_add(b);
            let (sum2, c2) = sum1.overflowing_add(carry);
            self.limbs[i] = sum2;
            carry = (c1 as u64) + (c2 as u64);
        }
        if carry > 0 {
            self.limbs.push(carry);
        }
    }

    /// Returns self + rhs.
    pub fn add(&self, rhs: &BigUint) -> BigUint {
        let mut result = self.clone();
        result.add_assign(rhs);
        result
    }

    // ─── Subtraction (assumes self >= rhs) ───

    /// self -= rhs. Panics if self < rhs.
    pub fn sub_assign(&mut self, rhs: &BigUint) {
        debug_assert!(*self >= *rhs, "BigUint::sub_assign: underflow");
        let mut borrow = 0u64;
        for i in 0..self.limbs.len() {
            let a = self.limbs[i];
            let b = if i < rhs.limbs.len() { rhs.limbs[i] } else { 0 };
            let (diff1, b1) = a.overflowing_sub(b);
            let (diff2, b2) = diff1.overflowing_sub(borrow);
            self.limbs[i] = diff2;
            borrow = (b1 as u64) + (b2 as u64);
        }
        debug_assert_eq!(borrow, 0, "BigUint::sub_assign: borrow remaining");
        self.normalize();
    }

    /// Returns self - rhs. Panics if self < rhs.
    pub fn sub(&self, rhs: &BigUint) -> BigUint {
        let mut result = self.clone();
        result.sub_assign(rhs);
        result
    }

    /// Checked subtraction: None if self < rhs.
    pub fn checked_sub(&self, rhs: &BigUint) -> Option<BigUint> {
        if self < rhs {
            None
        } else {
            Some(self.sub(rhs))
        }
    }

    // ─── Shift ───

    /// Left shift by `bits` positions.
    pub fn shl(&self, bits: u64) -> BigUint {
        if self.is_zero() || bits == 0 {
            return self.clone();
        }
        let limb_shift = (bits / 64) as usize;
        let bit_shift = (bits % 64) as u32;

        let mut result = vec![0u64; limb_shift + self.limbs.len() + 1];
        if bit_shift == 0 {
            result[limb_shift..(self.limbs.len() + limb_shift)].copy_from_slice(&self.limbs);
        } else {
            let mut carry = 0u64;
            for i in 0..self.limbs.len() {
                let shifted = self.limbs[i] << bit_shift;
                result[i + limb_shift] = shifted | carry;
                carry = self.limbs[i] >> (64 - bit_shift);
            }
            if carry > 0 {
                result[self.limbs.len() + limb_shift] = carry;
            }
        }
        BigUint::from_limbs(result)
    }

    /// Right shift by `bits` positions.
    pub fn shr(&self, bits: u64) -> BigUint {
        if self.is_zero() || bits == 0 {
            return self.clone();
        }
        let limb_shift = (bits / 64) as usize;
        let bit_shift = (bits % 64) as u32;

        if limb_shift >= self.limbs.len() {
            return BigUint::zero();
        }

        let new_len = self.limbs.len() - limb_shift;
        let mut result = vec![0u64; new_len];
        if bit_shift == 0 {
            result[..new_len].copy_from_slice(&self.limbs[limb_shift..(new_len + limb_shift)]);
        } else {
            for i in 0..new_len {
                result[i] = self.limbs[i + limb_shift] >> bit_shift;
                if i + limb_shift + 1 < self.limbs.len() {
                    result[i] |= self.limbs[i + limb_shift + 1] << (64 - bit_shift);
                }
            }
        }
        BigUint::from_limbs(result)
    }

    // ─── Multiplication ───

    /// Schoolbook multiplication: O(n*m).
    fn mul_schoolbook(&self, rhs: &BigUint) -> BigUint {
        if self.is_zero() || rhs.is_zero() {
            return BigUint::zero();
        }
        let n = self.limbs.len();
        let m = rhs.limbs.len();
        let mut result = vec![0u64; n + m];

        for i in 0..n {
            let mut carry = 0u128;
            for j in 0..m {
                let prod = (self.limbs[i] as u128) * (rhs.limbs[j] as u128)
                    + (result[i + j] as u128)
                    + carry;
                result[i + j] = prod as u64;
                carry = prod >> 64;
            }
            result[i + m] = carry as u64;
        }
        BigUint::from_limbs(result)
    }

    /// Karatsuba multiplication for large operands.
    fn mul_karatsuba(&self, rhs: &BigUint) -> BigUint {
        let n = self.limbs.len().max(rhs.limbs.len());

        // Base case: small enough for schoolbook
        if n < 32 {
            return self.mul_schoolbook(rhs);
        }

        let half = n / 2;

        // Split: self = lo0 + hi0 * B^half, rhs = lo1 + hi1 * B^half
        let (lo0, hi0) = split_at(self, half);
        let (lo1, hi1) = split_at(rhs, half);

        // z0 = lo0 * lo1
        let z0 = lo0.mul_karatsuba(&lo1);
        // z2 = hi0 * hi1
        let z2 = hi0.mul_karatsuba(&hi1);
        // z1 = (lo0 + hi0)(lo1 + hi1) - z0 - z2
        let sum0 = lo0.add(&hi0);
        let sum1 = lo1.add(&hi1);
        let z1_full = sum0.mul_karatsuba(&sum1);
        let z1 = z1_full.sub(&z0).sub(&z2);

        // result = z0 + z1 * B^half + z2 * B^(2*half)
        let shift = (half as u64) * 64;
        z0.add(&z1.shl(shift)).add(&z2.shl(shift * 2))
    }

    /// Multiply: picks schoolbook or Karatsuba based on size.
    pub fn mul(&self, rhs: &BigUint) -> BigUint {
        let n = self.limbs.len().max(rhs.limbs.len());
        if n < 32 {
            self.mul_schoolbook(rhs)
        } else {
            self.mul_karatsuba(rhs)
        }
    }

    /// Multiply by a single u64.
    pub fn mul_u64(&self, rhs: u64) -> BigUint {
        if self.is_zero() || rhs == 0 {
            return BigUint::zero();
        }
        if rhs == 1 {
            return self.clone();
        }
        let mut result = vec![0u64; self.limbs.len() + 1];
        let mut carry = 0u128;
        for i in 0..self.limbs.len() {
            let prod = (self.limbs[i] as u128) * (rhs as u128) + carry;
            result[i] = prod as u64;
            carry = prod >> 64;
        }
        result[self.limbs.len()] = carry as u64;
        BigUint::from_limbs(result)
    }

    // ─── Division ───

    /// Division and remainder: (quotient, remainder).
    /// Panics on division by zero.
    pub fn div_rem(&self, rhs: &BigUint) -> (BigUint, BigUint) {
        assert!(!rhs.is_zero(), "BigUint: division by zero");

        if self < rhs {
            return (BigUint::zero(), self.clone());
        }
        if rhs.limbs.len() == 1 {
            return self.div_rem_u64(rhs.limbs[0]);
        }

        // Knuth Algorithm D
        knuth_div(self, rhs)
    }

    /// Divide by a single u64. Returns (quotient, remainder).
    pub fn div_rem_u64(&self, rhs: u64) -> (BigUint, BigUint) {
        assert!(rhs != 0, "BigUint: division by zero");
        if self.is_zero() {
            return (BigUint::zero(), BigUint::zero());
        }

        let mut quotient = vec![0u64; self.limbs.len()];
        let mut rem = 0u128;
        for i in (0..self.limbs.len()).rev() {
            rem = (rem << 64) | (self.limbs[i] as u128);
            quotient[i] = (rem / rhs as u128) as u64;
            rem %= rhs as u128;
        }
        (BigUint::from_limbs(quotient), BigUint::from_u64(rem as u64))
    }

    pub fn div(&self, rhs: &BigUint) -> BigUint {
        self.div_rem(rhs).0
    }

    pub fn rem(&self, rhs: &BigUint) -> BigUint {
        self.div_rem(rhs).1
    }

    // ─── GCD ───

    /// Binary GCD (Stein's algorithm).
    pub fn gcd(&self, other: &BigUint) -> BigUint {
        if self.is_zero() {
            return other.clone();
        }
        if other.is_zero() {
            return self.clone();
        }

        let mut a = self.clone();
        let mut b = other.clone();

        // Factor out common powers of 2
        let a_tz = trailing_zeros(&a);
        let b_tz = trailing_zeros(&b);
        let shift = a_tz.min(b_tz);

        a = a.shr(a_tz);
        b = b.shr(b_tz);

        loop {
            // Both are odd at this point (or b is after the first iteration)
            let b_tz2 = trailing_zeros(&b);
            b = b.shr(b_tz2);

            // Ensure a <= b
            if a > b {
                std::mem::swap(&mut a, &mut b);
            }

            // b -= a (both odd, so result is even; trailing zeros handled next iter)
            b = b.sub(&a);

            if b.is_zero() {
                return a.shl(shift);
            }
        }
    }
}

/// Split a BigUint at limb position `at` into (lo, hi).
fn split_at(x: &BigUint, at: usize) -> (BigUint, BigUint) {
    if at >= x.limbs.len() {
        return (x.clone(), BigUint::zero());
    }
    let lo = BigUint::from_limbs(x.limbs[..at].to_vec());
    let hi = BigUint::from_limbs(x.limbs[at..].to_vec());
    (lo, hi)
}

/// Count trailing zeros in a BigUint.
fn trailing_zeros(x: &BigUint) -> u64 {
    if x.is_zero() {
        return 0;
    }
    for (i, &limb) in x.limbs.iter().enumerate() {
        if limb != 0 {
            return (i as u64) * 64 + limb.trailing_zeros() as u64;
        }
    }
    0
}

/// Knuth Algorithm D: multi-word division.
/// Returns (quotient, remainder).
fn knuth_div(u: &BigUint, v: &BigUint) -> (BigUint, BigUint) {
    let n = v.limbs.len();
    let m = u.limbs.len() - n;

    // D1: Normalize — shift so v's top bit is set.
    let shift = v.limbs[n - 1].leading_zeros();
    let vn = v.shl(shift as u64);
    let un = u.shl(shift as u64);

    // un may have gained a limb from the shift
    let mut un_limbs = un.limbs;
    if un_limbs.len() <= m + n {
        un_limbs.resize(m + n + 1, 0);
    }

    let mut q = vec![0u64; m + 1];

    for j in (0..=m).rev() {
        // D3: Estimate q_hat
        let u_hi = (un_limbs[j + n] as u128) << 64 | (un_limbs[j + n - 1] as u128);
        let mut q_hat = u_hi / (vn.limbs[n - 1] as u128);
        let mut r_hat = u_hi % (vn.limbs[n - 1] as u128);

        // Refine estimate
        loop {
            if q_hat >= (1u128 << 64) {
                q_hat -= 1;
                r_hat += vn.limbs[n - 1] as u128;
            } else if n >= 2 {
                let test = q_hat * (vn.limbs[n - 2] as u128);
                let compare = (r_hat << 64) | (un_limbs[j + n - 2] as u128);
                if test > compare {
                    q_hat -= 1;
                    r_hat += vn.limbs[n - 1] as u128;
                } else {
                    break;
                }
            } else {
                break;
            }
            if r_hat >= (1u128 << 64) {
                break;
            }
        }

        // D4: Multiply and subtract
        let mut borrow = 0i128;
        for i in 0..n {
            let prod = q_hat * (vn.limbs[i] as u128);
            let diff = (un_limbs[j + i] as i128) - (prod as u64 as i128) - borrow;
            un_limbs[j + i] = diff as u64;
            borrow = (prod >> 64) as i128 - (diff >> 64);
        }
        let diff = (un_limbs[j + n] as i128) - borrow;
        un_limbs[j + n] = diff as u64;

        // D5: Set quotient digit
        q[j] = q_hat as u64;

        // D6: If we subtracted too much, add back
        if diff < 0 {
            q[j] -= 1;
            let mut carry = 0u64;
            for i in 0..n {
                let sum = (un_limbs[j + i] as u128) + (vn.limbs[i] as u128) + (carry as u128);
                un_limbs[j + i] = sum as u64;
                carry = (sum >> 64) as u64;
            }
            un_limbs[j + n] = un_limbs[j + n].wrapping_add(carry);
        }
    }

    // D8: Unnormalize remainder
    un_limbs.truncate(n);
    let remainder = BigUint::from_limbs(un_limbs).shr(shift as u64);

    (BigUint::from_limbs(q), remainder)
}

// ─── Comparison ───

impl PartialEq for BigUint {
    fn eq(&self, other: &Self) -> bool {
        self.limbs == other.limbs
    }
}
impl Eq for BigUint {}

impl PartialOrd for BigUint {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for BigUint {
    fn cmp(&self, other: &Self) -> Ordering {
        let len_cmp = self.limbs.len().cmp(&other.limbs.len());
        if len_cmp != Ordering::Equal {
            return len_cmp;
        }
        // Same length: compare from most significant limb
        for i in (0..self.limbs.len()).rev() {
            let c = self.limbs[i].cmp(&other.limbs[i]);
            if c != Ordering::Equal {
                return c;
            }
        }
        Ordering::Equal
    }
}

impl std::hash::Hash for BigUint {
    fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
        self.limbs.hash(state);
    }
}

// ─── Display ───

impl std::fmt::Display for BigUint {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        if self.is_zero() {
            return write!(f, "0");
        }
        // Decimal conversion: repeatedly divide by 10^18 (largest power of 10 in u64)
        let ten18: u64 = 1_000_000_000_000_000_000;
        let mut chunks = Vec::new();
        let mut n = self.clone();
        while !n.is_zero() {
            let (q, r) = n.div_rem_u64(ten18);
            chunks.push(r.to_u128().unwrap_or(0) as u64);
            n = q;
        }
        // First chunk: no leading zeros
        if let Some(&last) = chunks.last() {
            write!(f, "{}", last)?;
            for &chunk in chunks.iter().rev().skip(1) {
                write!(f, "{:018}", chunk)?;
            }
        }
        Ok(())
    }
}

// ─── From/Into ───

impl From<u64> for BigUint {
    fn from(v: u64) -> Self {
        BigUint::from_u64(v)
    }
}

impl From<u128> for BigUint {
    fn from(v: u128) -> Self {
        BigUint::from_u128(v)
    }
}

// ═══════════════════════════════════════════════════════════
// BIGINT (signed)
// ═══════════════════════════════════════════════════════════

/// Signed arbitrary-precision integer: sign + magnitude.
#[derive(Clone, Debug)]
pub struct BigInt {
    sign: Sign,
    mag: BigUint,
}

#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
enum Sign {
    Zero,
    Positive,
    Negative,
}

impl BigInt {
    // ─── Constructors ───

    pub fn zero() -> Self {
        BigInt {
            sign: Sign::Zero,
            mag: BigUint::zero(),
        }
    }

    pub fn one() -> Self {
        BigInt {
            sign: Sign::Positive,
            mag: BigUint::one(),
        }
    }

    pub fn from_i128(v: i128) -> Self {
        if v == 0 {
            Self::zero()
        } else if v > 0 {
            BigInt {
                sign: Sign::Positive,
                mag: BigUint::from_u128(v as u128),
            }
        } else if v == i128::MIN {
            // -i128::MIN overflows; handle specially
            BigInt {
                sign: Sign::Negative,
                mag: BigUint::from_u128((i128::MAX as u128) + 1),
            }
        } else {
            BigInt {
                sign: Sign::Negative,
                mag: BigUint::from_u128((-v) as u128),
            }
        }
    }

    pub fn from_u128(v: u128) -> Self {
        if v == 0 {
            Self::zero()
        } else {
            BigInt {
                sign: Sign::Positive,
                mag: BigUint::from_u128(v),
            }
        }
    }

    /// Construct from sign and magnitude. Normalizes sign for zero magnitude.
    pub fn from_sign_mag(negative: bool, mag: BigUint) -> Self {
        if mag.is_zero() {
            BigInt {
                sign: Sign::Zero,
                mag,
            }
        } else if negative {
            BigInt {
                sign: Sign::Negative,
                mag,
            }
        } else {
            BigInt {
                sign: Sign::Positive,
                mag,
            }
        }
    }

    // ─── Accessors ───

    pub fn is_zero(&self) -> bool {
        self.sign == Sign::Zero
    }
    pub fn is_positive(&self) -> bool {
        self.sign == Sign::Positive
    }
    pub fn is_negative(&self) -> bool {
        self.sign == Sign::Negative
    }

    pub fn magnitude(&self) -> &BigUint {
        &self.mag
    }

    /// Absolute value.
    pub fn abs(&self) -> BigInt {
        BigInt::from_sign_mag(false, self.mag.clone())
    }

    pub fn to_i128(&self) -> Option<i128> {
        let v = self.mag.to_u128()?;
        match self.sign {
            Sign::Zero => Some(0),
            Sign::Positive => {
                if v <= i128::MAX as u128 {
                    Some(v as i128)
                } else {
                    None
                }
            }
            Sign::Negative => {
                if v == (i128::MAX as u128) + 1 {
                    Some(i128::MIN)
                } else if v <= i128::MAX as u128 {
                    Some(-(v as i128))
                } else {
                    None
                }
            }
        }
    }

    pub fn to_u128(&self) -> Option<u128> {
        match self.sign {
            Sign::Negative => None,
            _ => self.mag.to_u128(),
        }
    }

    // ─── Negation ───

    pub fn neg(&self) -> BigInt {
        BigInt {
            sign: match self.sign {
                Sign::Zero => Sign::Zero,
                Sign::Positive => Sign::Negative,
                Sign::Negative => Sign::Positive,
            },
            mag: self.mag.clone(),
        }
    }

    // ─── Arithmetic ───

    pub fn add(&self, rhs: &BigInt) -> BigInt {
        match (self.sign, rhs.sign) {
            (Sign::Zero, _) => rhs.clone(),
            (_, Sign::Zero) => self.clone(),
            (Sign::Positive, Sign::Positive) => {
                BigInt::from_sign_mag(false, self.mag.add(&rhs.mag))
            }
            (Sign::Negative, Sign::Negative) => BigInt::from_sign_mag(true, self.mag.add(&rhs.mag)),
            (Sign::Positive, Sign::Negative) => match self.mag.cmp(&rhs.mag) {
                Ordering::Greater => BigInt::from_sign_mag(false, self.mag.sub(&rhs.mag)),
                Ordering::Less => BigInt::from_sign_mag(true, rhs.mag.sub(&self.mag)),
                Ordering::Equal => BigInt::zero(),
            },
            (Sign::Negative, Sign::Positive) => match self.mag.cmp(&rhs.mag) {
                Ordering::Greater => BigInt::from_sign_mag(true, self.mag.sub(&rhs.mag)),
                Ordering::Less => BigInt::from_sign_mag(false, rhs.mag.sub(&self.mag)),
                Ordering::Equal => BigInt::zero(),
            },
        }
    }

    pub fn sub(&self, rhs: &BigInt) -> BigInt {
        self.add(&rhs.neg())
    }

    pub fn mul(&self, rhs: &BigInt) -> BigInt {
        if self.is_zero() || rhs.is_zero() {
            return BigInt::zero();
        }
        let neg = self.sign != rhs.sign;
        BigInt::from_sign_mag(neg, self.mag.mul(&rhs.mag))
    }

    /// Division truncated toward zero. Panics on division by zero.
    pub fn div(&self, rhs: &BigInt) -> BigInt {
        assert!(!rhs.is_zero(), "BigInt: division by zero");
        if self.is_zero() {
            return BigInt::zero();
        }
        let neg = self.sign != rhs.sign;
        let (q, _) = self.mag.div_rem(&rhs.mag);
        BigInt::from_sign_mag(neg, q)
    }

    /// Remainder (sign follows dividend). Panics on division by zero.
    pub fn rem(&self, rhs: &BigInt) -> BigInt {
        assert!(!rhs.is_zero(), "BigInt: division by zero");
        if self.is_zero() {
            return BigInt::zero();
        }
        let (_, r) = self.mag.div_rem(&rhs.mag);
        BigInt::from_sign_mag(self.is_negative(), r)
    }

    /// Division and remainder. Panics on division by zero.
    pub fn div_rem(&self, rhs: &BigInt) -> (BigInt, BigInt) {
        assert!(!rhs.is_zero(), "BigInt: division by zero");
        if self.is_zero() {
            return (BigInt::zero(), BigInt::zero());
        }
        let neg = self.sign != rhs.sign;
        let (q, r) = self.mag.div_rem(&rhs.mag);
        (
            BigInt::from_sign_mag(neg, q),
            BigInt::from_sign_mag(self.is_negative(), r),
        )
    }

    /// GCD (always non-negative).
    pub fn gcd(&self, other: &BigInt) -> BigInt {
        BigInt::from_sign_mag(false, self.mag.gcd(&other.mag))
    }
}

// ─── Comparison ───

impl PartialEq for BigInt {
    fn eq(&self, other: &Self) -> bool {
        self.sign == other.sign && self.mag == other.mag
    }
}
impl Eq for BigInt {}

impl PartialOrd for BigInt {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(other))
    }
}

impl Ord for BigInt {
    fn cmp(&self, other: &Self) -> Ordering {
        match (self.sign, other.sign) {
            (Sign::Zero, Sign::Zero) => Ordering::Equal,
            (Sign::Positive, Sign::Positive) => self.mag.cmp(&other.mag),
            (Sign::Negative, Sign::Negative) => other.mag.cmp(&self.mag),
            (Sign::Positive, _) => Ordering::Greater,
            (_, Sign::Positive) => Ordering::Less,
            (Sign::Zero, Sign::Negative) => Ordering::Greater,
            (Sign::Negative, Sign::Zero) => Ordering::Less,
        }
    }
}

impl std::hash::Hash for BigInt {
    fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
        self.sign.hash(state);
        self.mag.hash(state);
    }
}

// ─── Display ───

impl std::fmt::Display for BigInt {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self.sign {
            Sign::Zero => write!(f, "0"),
            Sign::Positive => write!(f, "{}", self.mag),
            Sign::Negative => write!(f, "-{}", self.mag),
        }
    }
}

// ─── Operator impls ───

impl std::ops::Neg for BigInt {
    type Output = BigInt;
    fn neg(self) -> BigInt {
        BigInt::neg(&self)
    }
}

impl std::ops::Neg for &BigInt {
    type Output = BigInt;
    fn neg(self) -> BigInt {
        BigInt::neg(self)
    }
}

impl std::ops::Add for &BigInt {
    type Output = BigInt;
    fn add(self, rhs: Self) -> BigInt {
        BigInt::add(self, rhs)
    }
}

impl std::ops::Sub for &BigInt {
    type Output = BigInt;
    fn sub(self, rhs: Self) -> BigInt {
        BigInt::sub(self, rhs)
    }
}

impl std::ops::Mul for &BigInt {
    type Output = BigInt;
    fn mul(self, rhs: Self) -> BigInt {
        BigInt::mul(self, rhs)
    }
}

// ─── From/Into ───

impl From<i128> for BigInt {
    fn from(v: i128) -> Self {
        BigInt::from_i128(v)
    }
}

impl From<u128> for BigInt {
    fn from(v: u128) -> Self {
        BigInt::from_u128(v)
    }
}

impl From<i64> for BigInt {
    fn from(v: i64) -> Self {
        BigInt::from_i128(v as i128)
    }
}

impl From<u64> for BigInt {
    fn from(v: u64) -> Self {
        BigInt::from_u128(v as u128)
    }
}

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

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

    // ─── BigUint basics ───

    #[test]
    fn biguint_zero() {
        let z = BigUint::zero();
        assert!(z.is_zero());
        assert_eq!(z.bit_len(), 0);
        assert_eq!(format!("{}", z), "0");
    }

    #[test]
    fn biguint_one() {
        let one = BigUint::one();
        assert!(!one.is_zero());
        assert!(one.is_one());
        assert_eq!(one.bit_len(), 1);
        assert_eq!(format!("{}", one), "1");
    }

    #[test]
    fn biguint_from_u128() {
        let v = BigUint::from_u128(u128::MAX);
        assert_eq!(v.limbs.len(), 2);
        assert_eq!(v.to_u128(), Some(u128::MAX));
    }

    #[test]
    fn biguint_from_u64_max() {
        let v = BigUint::from_u64(u64::MAX);
        assert_eq!(v.to_u128(), Some(u64::MAX as u128));
    }

    // ─── BigUint addition ───

    #[test]
    fn biguint_add_basic() {
        let a = BigUint::from_u64(100);
        let b = BigUint::from_u64(200);
        assert_eq!(a.add(&b).to_u128(), Some(300));
    }

    #[test]
    fn biguint_add_carry() {
        let a = BigUint::from_u64(u64::MAX);
        let b = BigUint::from_u64(1);
        let c = a.add(&b);
        assert_eq!(c.to_u128(), Some(u64::MAX as u128 + 1));
        assert_eq!(c.limbs.len(), 2);
    }

    #[test]
    fn biguint_add_large() {
        let a = BigUint::from_u128(u128::MAX);
        let b = BigUint::from_u128(1);
        let c = a.add(&b);
        assert_eq!(c.limbs.len(), 3);
        assert!(c.to_u128().is_none()); // exceeds u128
    }

    #[test]
    fn biguint_add_zero() {
        let a = BigUint::from_u64(42);
        assert_eq!(a.add(&BigUint::zero()), a);
        assert_eq!(BigUint::zero().add(&a), a);
    }

    // ─── BigUint subtraction ───

    #[test]
    fn biguint_sub_basic() {
        let a = BigUint::from_u64(300);
        let b = BigUint::from_u64(100);
        assert_eq!(a.sub(&b).to_u128(), Some(200));
    }

    #[test]
    fn biguint_sub_borrow() {
        let a = BigUint::from_u128((1u128 << 64) + 1);
        let b = BigUint::from_u64(2);
        let c = a.sub(&b);
        assert_eq!(c.to_u128(), Some(u64::MAX as u128));
    }

    #[test]
    fn biguint_sub_to_zero() {
        let a = BigUint::from_u64(42);
        assert!(a.sub(&a).is_zero());
    }

    #[test]
    fn biguint_checked_sub_underflow() {
        let a = BigUint::from_u64(1);
        let b = BigUint::from_u64(2);
        assert!(a.checked_sub(&b).is_none());
    }

    // ─── BigUint multiplication ───

    #[test]
    fn biguint_mul_basic() {
        let a = BigUint::from_u64(1000);
        let b = BigUint::from_u64(2000);
        assert_eq!(a.mul(&b).to_u128(), Some(2_000_000));
    }

    #[test]
    fn biguint_mul_large() {
        let a = BigUint::from_u64(u64::MAX);
        let b = BigUint::from_u64(u64::MAX);
        let c = a.mul(&b);
        let expected = (u64::MAX as u128) * (u64::MAX as u128);
        assert_eq!(c.to_u128(), Some(expected));
    }

    #[test]
    fn biguint_mul_zero() {
        let a = BigUint::from_u64(42);
        assert!(a.mul(&BigUint::zero()).is_zero());
        assert!(BigUint::zero().mul(&a).is_zero());
    }

    #[test]
    fn biguint_mul_one() {
        let a = BigUint::from_u64(42);
        assert_eq!(a.mul(&BigUint::one()), a);
    }

    #[test]
    fn biguint_mul_u64() {
        let a = BigUint::from_u128(u128::MAX);
        let b = a.mul_u64(2);
        // u128::MAX * 2 overflows u128
        assert!(b.to_u128().is_none());
        assert_eq!(b.limbs.len(), 3);
    }

    #[test]
    fn biguint_mul_karatsuba_matches_schoolbook() {
        // Build operands large enough to trigger Karatsuba (32+ limbs)
        let a_limbs: Vec<u64> = (1..=40).collect();
        let b_limbs: Vec<u64> = (41..=80).collect();
        let a = BigUint::from_limbs(a_limbs);
        let b = BigUint::from_limbs(b_limbs);

        let school = a.mul_schoolbook(&b);
        let karat = a.mul_karatsuba(&b);
        assert_eq!(school, karat);
    }

    // ─── BigUint division ───

    #[test]
    fn biguint_div_basic() {
        let a = BigUint::from_u64(100);
        let b = BigUint::from_u64(7);
        let (q, r) = a.div_rem(&b);
        assert_eq!(q.to_u128(), Some(14));
        assert_eq!(r.to_u128(), Some(2));
    }

    #[test]
    fn biguint_div_exact() {
        let a = BigUint::from_u64(100);
        let b = BigUint::from_u64(10);
        let (q, r) = a.div_rem(&b);
        assert_eq!(q.to_u128(), Some(10));
        assert!(r.is_zero());
    }

    #[test]
    fn biguint_div_larger_divisor() {
        let a = BigUint::from_u64(5);
        let b = BigUint::from_u64(10);
        let (q, r) = a.div_rem(&b);
        assert!(q.is_zero());
        assert_eq!(r.to_u128(), Some(5));
    }

    #[test]
    fn biguint_div_multiword() {
        let a = BigUint::from_u128(u128::MAX);
        let b = BigUint::from_u64(3);
        let (q, r) = a.div_rem(&b);
        let qa = q.mul(&b).add(&r);
        assert_eq!(qa, a); // q*b + r == a
    }

    #[test]
    fn biguint_div_multiword_by_multiword() {
        let a = BigUint::from_u128(u128::MAX);
        let b = BigUint::from_u128(u64::MAX as u128 + 2); // two-limb divisor
        let (q, r) = a.div_rem(&b);
        let qa = q.mul(&b).add(&r);
        assert_eq!(qa, a);
    }

    #[test]
    #[should_panic(expected = "division by zero")]
    fn biguint_div_by_zero() {
        let _ = BigUint::from_u64(1).div_rem(&BigUint::zero());
    }

    // ─── BigUint GCD ───

    #[test]
    fn biguint_gcd_basic() {
        let a = BigUint::from_u64(48);
        let b = BigUint::from_u64(18);
        assert_eq!(a.gcd(&b).to_u128(), Some(6));
    }

    #[test]
    fn biguint_gcd_coprime() {
        let a = BigUint::from_u64(17);
        let b = BigUint::from_u64(13);
        assert!(a.gcd(&b).is_one());
    }

    #[test]
    fn biguint_gcd_with_zero() {
        let a = BigUint::from_u64(42);
        assert_eq!(BigUint::zero().gcd(&a), a);
        assert_eq!(a.gcd(&BigUint::zero()), a);
    }

    #[test]
    fn biguint_gcd_power_of_two() {
        let a = BigUint::from_u64(64);
        let b = BigUint::from_u64(48);
        assert_eq!(a.gcd(&b).to_u128(), Some(16));
    }

    #[test]
    fn biguint_gcd_large() {
        let a = BigUint::from_u128(u128::MAX);
        let b = BigUint::from_u128(u128::MAX - 1);
        let g = a.gcd(&b);
        assert!(g.is_one());
    }

    // ─── BigUint shift ───

    #[test]
    fn biguint_shl_basic() {
        let a = BigUint::from_u64(1);
        assert_eq!(a.shl(64).to_u128(), Some(1u128 << 64));
    }

    #[test]
    fn biguint_shr_basic() {
        let a = BigUint::from_u128(1u128 << 64);
        assert_eq!(a.shr(64).to_u128(), Some(1));
    }

    #[test]
    fn biguint_shl_shr_roundtrip() {
        let a = BigUint::from_u128(12345678901234567890);
        assert_eq!(a.shl(37).shr(37), a);
    }

    #[test]
    fn biguint_shr_to_zero() {
        let a = BigUint::from_u64(255);
        assert!(a.shr(100).is_zero());
    }

    // ─── BigUint comparison ───

    #[test]
    fn biguint_ord() {
        let a = BigUint::from_u64(100);
        let b = BigUint::from_u64(200);
        assert!(a < b);
        assert!(b > a);
        assert_eq!(a, a);
    }

    #[test]
    fn biguint_ord_different_lengths() {
        let a = BigUint::from_u64(u64::MAX);
        let b = BigUint::from_u128((u64::MAX as u128) + 1);
        assert!(a < b);
    }

    // ─── BigUint display ───

    #[test]
    fn biguint_display_small() {
        assert_eq!(format!("{}", BigUint::from_u64(12345)), "12345");
    }

    #[test]
    fn biguint_display_large() {
        let a = BigUint::from_u128(u128::MAX);
        assert_eq!(format!("{}", a), format!("{}", u128::MAX));
    }

    // ─── BigUint bit operations ───

    #[test]
    fn biguint_bit() {
        let a = BigUint::from_u64(0b1010);
        assert!(!a.bit(0));
        assert!(a.bit(1));
        assert!(!a.bit(2));
        assert!(a.bit(3));
        assert!(!a.bit(100)); // out of range
    }

    #[test]
    fn biguint_is_even() {
        assert!(BigUint::zero().is_even());
        assert!(!BigUint::one().is_even());
        assert!(BigUint::from_u64(42).is_even());
        assert!(!BigUint::from_u64(43).is_even());
    }

    // ─── BigInt basics ───

    #[test]
    fn bigint_zero() {
        let z = BigInt::zero();
        assert!(z.is_zero());
        assert!(!z.is_positive());
        assert!(!z.is_negative());
        assert_eq!(z.to_i128(), Some(0));
    }

    #[test]
    fn bigint_from_positive() {
        let a = BigInt::from_i128(42);
        assert!(a.is_positive());
        assert_eq!(a.to_i128(), Some(42));
    }

    #[test]
    fn bigint_from_negative() {
        let a = BigInt::from_i128(-42);
        assert!(a.is_negative());
        assert_eq!(a.to_i128(), Some(-42));
    }

    #[test]
    fn bigint_from_i128_min() {
        let a = BigInt::from_i128(i128::MIN);
        assert!(a.is_negative());
        assert_eq!(a.to_i128(), Some(i128::MIN));
    }

    #[test]
    fn bigint_from_i128_max() {
        let a = BigInt::from_i128(i128::MAX);
        assert!(a.is_positive());
        assert_eq!(a.to_i128(), Some(i128::MAX));
    }

    // ─── BigInt addition ───

    #[test]
    fn bigint_add_positive() {
        let a = BigInt::from_i128(100);
        let b = BigInt::from_i128(200);
        assert_eq!(a.add(&b).to_i128(), Some(300));
    }

    #[test]
    fn bigint_add_negative() {
        let a = BigInt::from_i128(-100);
        let b = BigInt::from_i128(-200);
        assert_eq!(a.add(&b).to_i128(), Some(-300));
    }

    #[test]
    fn bigint_add_mixed_positive_result() {
        let a = BigInt::from_i128(100);
        let b = BigInt::from_i128(-30);
        assert_eq!(a.add(&b).to_i128(), Some(70));
    }

    #[test]
    fn bigint_add_mixed_negative_result() {
        let a = BigInt::from_i128(30);
        let b = BigInt::from_i128(-100);
        assert_eq!(a.add(&b).to_i128(), Some(-70));
    }

    #[test]
    fn bigint_add_cancellation() {
        let a = BigInt::from_i128(42);
        let b = BigInt::from_i128(-42);
        assert!(a.add(&b).is_zero());
    }

    #[test]
    fn bigint_add_overflow_i128() {
        let a = BigInt::from_i128(i128::MAX);
        let b = BigInt::from_i128(1);
        let c = a.add(&b);
        // Can't fit in i128 anymore
        assert!(c.to_i128().is_none());
        // But it's correct as BigInt
        assert!(c.is_positive());
        // And subtracting 1 gets us back
        assert_eq!(c.sub(&b).to_i128(), Some(i128::MAX));
    }

    // ─── BigInt subtraction ───

    #[test]
    fn bigint_sub() {
        let a = BigInt::from_i128(100);
        let b = BigInt::from_i128(30);
        assert_eq!(a.sub(&b).to_i128(), Some(70));
    }

    // ─── BigInt multiplication ───

    #[test]
    fn bigint_mul_positive() {
        let a = BigInt::from_i128(100);
        let b = BigInt::from_i128(200);
        assert_eq!(a.mul(&b).to_i128(), Some(20000));
    }

    #[test]
    fn bigint_mul_sign_rules() {
        let pos = BigInt::from_i128(3);
        let neg = BigInt::from_i128(-5);
        assert_eq!(pos.mul(&neg).to_i128(), Some(-15));
        assert_eq!(neg.mul(&pos).to_i128(), Some(-15));
        assert_eq!(neg.mul(&neg).to_i128(), Some(25)); // (-5)*(-5) = 25
    }

    #[test]
    fn bigint_mul_zero() {
        let a = BigInt::from_i128(42);
        assert!(a.mul(&BigInt::zero()).is_zero());
    }

    // ─── BigInt division ───

    #[test]
    fn bigint_div_basic() {
        let a = BigInt::from_i128(100);
        let b = BigInt::from_i128(7);
        let (q, r) = a.div_rem(&b);
        assert_eq!(q.to_i128(), Some(14));
        assert_eq!(r.to_i128(), Some(2));
    }

    #[test]
    fn bigint_div_negative() {
        let a = BigInt::from_i128(-100);
        let b = BigInt::from_i128(7);
        let (q, r) = a.div_rem(&b);
        assert_eq!(q.to_i128(), Some(-14));
        assert_eq!(r.to_i128(), Some(-2));
    }

    // ─── BigInt GCD ───

    #[test]
    fn bigint_gcd() {
        let a = BigInt::from_i128(-48);
        let b = BigInt::from_i128(18);
        let g = a.gcd(&b);
        assert_eq!(g.to_i128(), Some(6));
        assert!(g.is_positive());
    }

    // ─── BigInt ordering ───

    #[test]
    fn bigint_ord() {
        let neg = BigInt::from_i128(-100);
        let zero = BigInt::zero();
        let pos = BigInt::from_i128(100);
        assert!(neg < zero);
        assert!(zero < pos);
        assert!(neg < pos);
    }

    // ─── BigInt display ───

    #[test]
    fn bigint_display() {
        assert_eq!(format!("{}", BigInt::from_i128(42)), "42");
        assert_eq!(format!("{}", BigInt::from_i128(-42)), "-42");
        assert_eq!(format!("{}", BigInt::zero()), "0");
    }

    // ─── Division stress: Fibonacci / factorials ───

    #[test]
    fn biguint_fibonacci_stress() {
        // Compute fib(100) and verify against known value
        let mut a = BigUint::from_u64(0);
        let mut b = BigUint::from_u64(1);
        for _ in 0..100 {
            let c = a.add(&b);
            a = b;
            b = c;
        }
        // After 100 iterations: a = fib(100), b = fib(101)
        // fib(100) = 354224848179261915075
        assert_eq!(format!("{}", a), "354224848179261915075");
    }

    #[test]
    fn biguint_factorial_stress() {
        // 20! = 2432902008176640000
        let mut f = BigUint::one();
        for i in 1..=20u64 {
            f = f.mul_u64(i);
        }
        assert_eq!(f.to_u128(), Some(2432902008176640000u128));
    }

    #[test]
    fn biguint_factorial_50() {
        // 50! is ~65 digits, well beyond u128
        let mut f = BigUint::one();
        for i in 1..=50u64 {
            f = f.mul_u64(i);
        }
        assert!(f.to_u128().is_none()); // doesn't fit
                                        // Verify by dividing back down
        for i in (1..=50u64).rev() {
            let (q, r) = f.div_rem_u64(i);
            assert!(r.is_zero(), "50! should be divisible by {}", i);
            f = q;
        }
        assert!(f.is_one());
    }

    #[test]
    fn biguint_div_roundtrip_stress() {
        // For various a, b: verify a == (a/b)*b + (a%b)
        let values: Vec<u128> = vec![
            0,
            1,
            2,
            u64::MAX as u128,
            u64::MAX as u128 + 1,
            u128::MAX,
            u128::MAX / 3,
            999999999999999999,
        ];
        for &a_val in &values {
            for &b_val in &values {
                if b_val == 0 {
                    continue;
                }
                let a = BigUint::from_u128(a_val);
                let b = BigUint::from_u128(b_val);
                let (q, r) = a.div_rem(&b);
                let reconstruct = q.mul(&b).add(&r);
                assert_eq!(
                    reconstruct, a,
                    "div_rem roundtrip failed: {} / {}",
                    a_val, b_val
                );
                assert!(r < b, "remainder >= divisor: {} % {}", a_val, b_val);
            }
        }
    }
}