ogdoad 1.0.4

//! Commutative coefficient rings and their optional algebraic capabilities.
//!
//! [`Scalar`] is the borrow-based ring interface used by the generic Clifford
//! engine. Concrete backends are grouped by mathematical role and re-exported
//! from this module:
//!
//! - [`exact`]: fixed-width `Integer` and `Rational` arithmetic;
//! - [`finite_field`]: prime and extension fields, nimbers, and truncated Witt
//!   vectors;
//! - [`big`]: finite-support surreal, omnific, and checked ordinal-nimber models;
//! - [`small`]: capped-relative `p`-adic and unramified local models;
//! - [`global`]: rational-function and represented adele models;
//! - [`functor`]: complex, Laurent, ramified, and Gauss adjunctions.
//!
//! Exactness, valuations, residue fields, integral subrings, extensions, and
//! root-taking are opt-in traits; they are not assumptions of [`Scalar`].
//! [`Tropical`] is deliberately separate because an idempotent semiring has no
//! additive inverses. Arbitrary partizan games likewise do not implement
//! `Scalar`: they form an additive group, not a commutative ring.

pub mod analytic;
pub mod big;
pub mod exact;
pub mod exactness;
pub mod extension;
pub mod finite_field;
pub mod functor;
pub mod global;
pub mod integrality;
pub mod newton;
pub mod poly;
pub(crate) mod poly_factor;
pub mod residue;
pub mod small;
pub mod tropical;
pub mod valued;

pub use analytic::*;
pub use big::*;
pub use exact::*;
pub use exactness::*;
pub use extension::*;
pub use finite_field::*;
pub use functor::*;
pub use global::*;
pub use integrality::*;
pub use newton::*;
pub use poly::*;
pub use residue::*;
pub use small::*;
pub use tropical::*;
pub use valued::*;

use std::fmt::{Debug, Display};
use std::ops::{Add, BitXor, Mul, Neg, Sub};

pub(crate) fn mod_inverse_u128(a: u128, modulus: u128) -> Option<u128> {
    if modulus <= 1 {
        return None;
    }
    let (mut t, mut new_t) = (0u128, 1u128 % modulus);
    let (mut r, mut new_r) = (modulus, a % modulus);
    while new_r != 0 {
        let quotient = r / new_r;
        let q_new_t = mul_mod_u128(quotient % modulus, new_t, modulus);
        let next_t = sub_mod_u128(t, q_new_t, modulus);
        std::mem::swap(&mut t, &mut new_t);
        new_t = next_t;
        r -= quotient * new_r;
        std::mem::swap(&mut r, &mut new_r);
    }
    (r == 1).then_some(t)
}

pub(crate) fn add_mod_u128(a: u128, b: u128, modulus: u128) -> u128 {
    debug_assert!(modulus > 0);
    let a = a % modulus;
    let b = b % modulus;
    if a >= modulus - b {
        a - (modulus - b)
    } else {
        a + b
    }
}

pub(crate) fn sub_mod_u128(a: u128, b: u128, modulus: u128) -> u128 {
    debug_assert!(modulus > 0);
    let a = a % modulus;
    let b = b % modulus;
    if a >= b {
        a - b
    } else {
        modulus - (b - a)
    }
}

pub(crate) fn mul_mod_u128(mut a: u128, mut b: u128, modulus: u128) -> u128 {
    debug_assert!(modulus > 0);
    if modulus == 1 {
        return 0;
    }
    a %= modulus;
    let mut acc = 0u128;
    while b > 0 {
        if b & 1 == 1 {
            acc = add_mod_u128(acc, a, modulus);
        }
        b >>= 1;
        if b > 0 {
            a = add_mod_u128(a, a, modulus);
        }
    }
    acc
}

pub(crate) fn reduce_i128_mod_u128(n: i128, modulus: u128) -> u128 {
    debug_assert!(modulus > 0);
    if n >= 0 {
        (n as u128) % modulus
    } else {
        let r = n.unsigned_abs() % modulus;
        if r == 0 {
            0
        } else {
            modulus - r
        }
    }
}

pub(crate) fn is_prime_u128(p: u128) -> bool {
    if p < 2 {
        return false;
    }
    if p.is_multiple_of(2) {
        return p == 2;
    }
    let mut d = 3u128;
    while d <= p / d {
        if p.is_multiple_of(d) {
            return false;
        }
        d += 2;
    }
    true
}

/// Checked factorial in the exact `i128` carrier.
///
/// Negative inputs are outside the factorial domain. `33!` is the largest
/// factorial represented by `i128`; `34!` returns `None`.
pub fn checked_factorial_i128(n: i128) -> Option<i128> {
    if n < 0 {
        return None;
    }
    let mut acc = 1i128;
    for k in 2..=n {
        acc = acc.checked_mul(k)?;
    }
    Some(acc)
}

/// Factorial computed inside a scalar world via the `Z -> S` ring map.
///
/// This is the finite-field-friendly path for grundy `!n`: in positive
/// characteristic the product is immediately zero once a factor equal to the
/// characteristic appears, so no host integer overflow is involved.
pub fn factorial_in_scalar<S: Scalar>(n: i128) -> Option<S> {
    if n < 0 {
        return None;
    }
    let characteristic = S::characteristic();
    if characteristic > 0 && n.unsigned_abs() >= characteristic {
        return Some(S::zero());
    }
    let mut acc = S::one();
    for k in 2..=n {
        acc = acc.mul(&S::from_int(k));
    }
    Some(acc)
}

/// Generate the owned-value operators `+`, `-` (binary and unary), `*`, and
/// `^ u128` (power) for a [`Scalar`] backend by forwarding to its trait
/// methods, so downstream code can write `a + b`, `a * b`, `-a`, `a ^ 3`
/// instead of `a.add(&b)`, `a.mul(&b)`, `a.neg()`, and an explicit loop.
/// Only use this for backends whose multiplication is total on the represented
/// domain; `Ordinal` gets hand-written additive operators and keeps its partial
/// product behind `nim_mul` / `nim_pow`.
///
/// **`^` is power (grundy `↑`), not XOR.** The RHS is deliberately `u128`
/// so that `x ^ y` never compiles when `y` has the same element type as `x` —
/// the type system enforces the "no element-element XOR" rule (on `Nimber`,
/// `x ^ x` would silently mean nim-*addition*). The exponent is an unsigned
/// meta-integer: `x ^ 0 == one()`.
///
/// **Precedence caveat (§5 `grundy/docs/spec.md`):** Rust's `^` binds looser than
/// `*`. `a * b ^ 3` is `a * (b ^ 3)` in grundy but `(a * b) ^ 3` in Rust.
/// Parenthesize when mixing product and power operators.
///
/// Deliberately *not* a [`Scalar`] supertrait bound: these are concrete-type
/// conveniences for callers (`Surreal + Surreal`, `-nimber`), so generic engine
/// code over `S: Scalar` keeps resolving `.add(&x)` / `.mul(&x)` to the `&self`
/// trait methods — operators-on-`S` would shadow them at owned-receiver sites and
/// force clones the borrow-based engine avoids. Division stays a method
/// ([`Scalar::inv`] is partial — `Div` would have to panic), and the
/// by-reference operator forms are omitted for the same reason.
///
/// The generic-backend form takes its generic clause in brackets, e.g.
/// `impl_scalar_ops!([const P: u128] Fp<P>)` or `impl_scalar_ops!([S: Scalar]
/// Surcomplex<S>)`; the bare form `impl_scalar_ops!(Rational)` is for the
/// concrete ones. (Brackets, not `<…>`, so the matcher stays unambiguous.)
macro_rules! impl_scalar_ops {
    ([$($gen:tt)*] $ty:ty) => {
        impl<$($gen)*> Add for $ty {
            type Output = $ty;
            #[inline]
            fn add(self, rhs: $ty) -> $ty { <$ty as $crate::scalar::Scalar>::add(&self, &rhs) }
        }
        impl<$($gen)*> Sub for $ty {
            type Output = $ty;
            #[inline]
            fn sub(self, rhs: $ty) -> $ty { <$ty as $crate::scalar::Scalar>::sub(&self, &rhs) }
        }
        impl<$($gen)*> Mul for $ty {
            type Output = $ty;
            #[inline]
            fn mul(self, rhs: $ty) -> $ty { <$ty as $crate::scalar::Scalar>::mul(&self, &rhs) }
        }
        impl<$($gen)*> Neg for $ty {
            type Output = $ty;
            #[inline]
            fn neg(self) -> $ty { <$ty as $crate::scalar::Scalar>::neg(&self) }
        }
        impl<$($gen)*> BitXor<u128> for $ty {
            type Output = $ty;
            /// Square-and-multiply power: `x ^ 0 == one()`, `x ^ k` via [`Scalar::pow`].
            ///
            /// `^` is power (grundy `↑`). The RHS is `u128` so element-element `^`
            /// does not compile — no [`BitXor<Self>`] impl exists on any backend.
            /// **Precedence caveat:** Rust's `^` binds looser than `*`; parenthesize
            /// when mixing with product.
            #[inline]
            fn bitxor(self, k: u128) -> $ty {
                <$ty as $crate::scalar::Scalar>::pow(&self, k)
            }
        }
    };
    ($ty:ty) => {
        impl Add for $ty {
            type Output = $ty;
            #[inline]
            fn add(self, rhs: $ty) -> $ty { <$ty as $crate::scalar::Scalar>::add(&self, &rhs) }
        }
        impl Sub for $ty {
            type Output = $ty;
            #[inline]
            fn sub(self, rhs: $ty) -> $ty { <$ty as $crate::scalar::Scalar>::sub(&self, &rhs) }
        }
        impl Mul for $ty {
            type Output = $ty;
            #[inline]
            fn mul(self, rhs: $ty) -> $ty { <$ty as $crate::scalar::Scalar>::mul(&self, &rhs) }
        }
        impl Neg for $ty {
            type Output = $ty;
            #[inline]
            fn neg(self) -> $ty { <$ty as $crate::scalar::Scalar>::neg(&self) }
        }
        impl BitXor<u128> for $ty {
            type Output = $ty;
            /// Square-and-multiply power: `x ^ 0 == one()`, `x ^ k` via [`Scalar::pow`].
            ///
            /// `^` is power (grundy `↑`). The RHS is `u128` so element-element `^`
            /// does not compile — no [`BitXor<Self>`] impl exists on any backend.
            /// **Precedence caveat:** Rust's `^` binds looser than `*`; parenthesize
            /// when mixing with product.
            #[inline]
            fn bitxor(self, k: u128) -> $ty {
                <$ty as $crate::scalar::Scalar>::pow(&self, k)
            }
        }
    };
}

/// Borrow-based interface for a commutative ring with a partial inverse operation.
///
/// Generic algebra code uses these methods instead of concrete operator impls.
/// [`inv`](Self::inv) returns `None` for nonunits and for inverses outside a
/// backend's represented domain.
pub trait Scalar: Clone + PartialEq + Debug + Display {
    /// The additive identity.
    fn zero() -> Self;
    /// The multiplicative identity.
    fn one() -> Self;
    /// Ring addition.
    fn add(&self, rhs: &Self) -> Self;
    /// Additive inverse.
    fn neg(&self) -> Self;
    /// Ring multiplication.
    fn mul(&self, rhs: &Self) -> Self;

    /// Ring characteristic: 0 for characteristic-0 domains, a positive additive
    /// order of `1` for finite fields and finite quotient rings (`Z/p^k`,
    /// truncated Witt vectors, etc.). The engine itself gets signs from
    /// [`Scalar::neg`]; callers that care about characteristic must distinguish
    /// fields from local rings separately.
    fn characteristic() -> u128;

    /// Multiplicative inverse, or `None` if not invertible (zero) or not
    /// finitely representable in this backend (e.g. a non-monomial surreal,
    /// whose inverse is an infinite Hahn series).
    fn inv(&self) -> Option<Self>;

    /// Whether this element is the additive identity.
    fn is_zero(&self) -> bool {
        *self == Self::zero()
    }

    /// Ring subtraction.
    fn sub(&self, rhs: &Self) -> Self {
        self.add(&rhs.neg())
    }

    /// The unital ring homomorphism ℤ → R (the unique ring homomorphism from
    /// the initial ring ℤ into any unital ring).
    ///
    /// The default implementation uses double-and-add over [`Scalar::one`] and
    /// [`Scalar::neg`], so **for characteristic-2 worlds `from_int(n) = n mod 2`**
    /// automatically — do NOT override for `Nimber`/`Ordinal` with a bit-cast of
    /// `n as u128`, which would produce a REPRESENTATION constructor (which nimber)
    /// rather than the ℤ-embedding. Override only where a direct construction is
    /// faster AND semantically identical (e.g. `Rational::from_int(n)`, `Integer(n)`).
    fn from_int(n: i128) -> Self {
        if n == 0 {
            return Self::zero();
        }
        let neg = n < 0;
        let abs = n.unsigned_abs();
        // double-and-add
        let mut base = Self::one();
        let mut acc = Self::zero();
        let mut remaining = abs;
        while remaining > 0 {
            if remaining & 1 == 1 {
                acc = acc.add(&base);
            }
            remaining >>= 1;
            if remaining > 0 {
                base = base.add(&base);
            }
        }
        if neg {
            acc.neg()
        } else {
            acc
        }
    }

    /// `self^exp` by square-and-multiply over [`Scalar::mul`]/[`Scalar::one`];
    /// `x.pow(0) == one()`. The same default-method precedent as
    /// [`Scalar::from_int`]: one correct implementation for every backend whose
    /// `mul` is total, with per-backend overrides only where a genuinely sharper
    /// algorithm exists (e.g. [`Nimber`]'s Fermat-tower `nim_pow`, reached through
    /// [`FiniteField::pow`]).
    ///
    /// `Ordinal` implements `Scalar` with a panic-on-escape `mul` (the represented
    /// Kummer tower's honest boundary), so it inherits this default `pow` as a
    /// checked/panicking path — consistent with its `Scalar` impl, but the
    /// concrete-type `^` operator (the `impl_scalar_ops!` macro) stays
    /// deliberately absent on `Ordinal`. This method is the generic entry point
    /// every other backend's `^` forwards to.
    fn pow(&self, exp: u128) -> Self {
        if exp == 0 {
            return Self::one();
        }
        let mut acc = Self::one();
        let mut base = self.clone();
        let mut e = exp;
        while e > 0 {
            if e & 1 == 1 {
                acc = acc.mul(&base);
            }
            e >>= 1;
            if e > 0 {
                base = base.mul(&base);
            }
        }
        acc
    }
}

// The operator manifest: every backend gets `+ - *` and unary `-` forwarded to
// its `Scalar` methods, so the whole table reads uniformly (`a + b`, `-a`). The
// const-/type-generic backends carry their generic clause; the concrete ones
// don't. (See [`impl_scalar_ops`].)
impl_scalar_ops!(Rational);
impl_scalar_ops!(Integer);
impl_scalar_ops!(Surreal);
impl_scalar_ops!(Omnific);
impl_scalar_ops!(Nimber);
impl Add for Ordinal {
    type Output = Ordinal;
    #[inline]
    fn add(self, rhs: Ordinal) -> Ordinal {
        <Ordinal as Scalar>::add(&self, &rhs)
    }
}
impl Sub for Ordinal {
    type Output = Ordinal;
    #[inline]
    fn sub(self, rhs: Ordinal) -> Ordinal {
        <Ordinal as Scalar>::sub(&self, &rhs)
    }
}
impl Neg for Ordinal {
    type Output = Ordinal;
    #[inline]
    fn neg(self) -> Ordinal {
        <Ordinal as Scalar>::neg(&self)
    }
}
impl_scalar_ops!([const P: u128] Fp<P>);
impl_scalar_ops!([const P: u128, const N: usize] Fpn<P, N>);
impl_scalar_ops!([const P: u128, const N: usize, const F: usize] WittVec<P, N, F>);
impl_scalar_ops!([const P: u128, const K: u128] Qp<P, K>);
impl_scalar_ops!([const P: u128, const K: u128] Zp<P, K>);
impl_scalar_ops!([const P: u128, const N: usize, const F: usize] Qq<P, N, F>);
impl_scalar_ops!([S: Scalar] Surcomplex<S>);
impl_scalar_ops!([S: Scalar, const K: usize] Laurent<S, K>);
impl_scalar_ops!([S: Valued, const E: usize] Ramified<S, E>);
impl_scalar_ops!([S: Valued] Gauss<S>);
impl_scalar_ops!(Adele);
impl_scalar_ops!([S: crate::scalar::ExactFieldScalar] RationalFunction<S>);
impl_scalar_ops!([S: Scalar] Poly<S>);

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

    /// Operators must agree with the `Scalar` trait methods they forward to —
    /// over a char-0 field and a char-2 one (where `-a = a`).
    #[test]
    fn operators_match_trait_methods() {
        let (a, b) = (Rational::new(2, 3), Rational::new(1, 6));
        assert_eq!(a.clone() + b.clone(), Scalar::add(&a, &b));
        assert_eq!(a.clone() - b.clone(), Scalar::sub(&a, &b));
        assert_eq!(a.clone() * b.clone(), Scalar::mul(&a, &b));
        assert_eq!(-a.clone(), Scalar::neg(&a));
        assert_eq!(a.clone() - a.clone(), Rational::zero());

        // char 2: `+` is XOR, `*` the nim product, and unary `-` is identity.
        let (x, y) = (Nimber(6), Nimber(3));
        assert_eq!(x + y, Scalar::add(&x, &y));
        assert_eq!(x * y, Scalar::mul(&x, &y));
        assert_eq!(-x, x);
    }

    #[test]
    fn scalar_power_operator_basic_cases() {
        // Nimber: *2 ^ 2 = nim_mul(2,2) = 3  (since 2*2=3 in nim arithmetic)
        // i.e. Nimber(2) ^ 2 == Nimber(3)
        assert_eq!(Nimber(2) ^ 2u128, Nimber(3));
        // x ^ 0 == one() for all total-product backends
        assert_eq!(Nimber(5) ^ 0u128, Nimber::one());
        assert_eq!(Rational::from_int(7) ^ 0u128, Rational::one());
        // Fp case: 3 ^ 2 == 9 mod p = 2 in F_5
        use crate::scalar::Fp;
        let three: Fp<5> = Fp::from_int(3);
        assert_eq!(three ^ 2u128, Fp::from_int(4)); // 3^2 = 9 ≡ 4 mod 5
                                                    // consistency with repeated mul
        let r2 = Rational::from_int(2);
        let r8 = Rational::from_int(8);
        assert_eq!(r2 ^ 3u128, r8);
    }

    /// Reaches [`Scalar::pow`] through the trait bound alone — no `impl_scalar_ops!`
    /// `^` operator is reachable inside a function generic only over `S: Scalar`,
    /// so this pins the default method itself, independent of any concrete-type
    /// operator convenience.
    fn generic_pow_via_trait<S: Scalar>(x: &S, n: u128) -> S {
        x.pow(n)
    }

    #[test]
    fn scalar_pow_default_matches_repeated_mul_generically() {
        let r = Rational::new(2, 3);
        let mut expected = Rational::one();
        for _ in 0..4 {
            expected = Scalar::mul(&expected, &r);
        }
        assert_eq!(generic_pow_via_trait(&r, 4), expected);
        assert_eq!(generic_pow_via_trait(&r, 0), Rational::one());

        let x = Nimber(5);
        let mut expected_n = Nimber::one();
        for _ in 0..3 {
            expected_n = Scalar::mul(&expected_n, &x);
        }
        assert_eq!(generic_pow_via_trait(&x, 3), expected_n);
    }

    #[test]
    fn ordinal_pow_default_method_has_no_operator_to_fall_back_on() {
        // `Ordinal` deliberately carries no `^` operator (multiplication is
        // checked/partial at the Kummer boundary — see `impl_scalar_ops!`'s
        // doc), so `.pow` here can only be the `Scalar` trait default, not an
        // operator-forwarding convenience. Stick to 0/1 so the checked `mul`
        // never has a chance to escape the verified boundary.
        assert_eq!(Ordinal::zero().pow(0), Ordinal::one());
        assert_eq!(Ordinal::zero().pow(3), Ordinal::zero());
        assert_eq!(Ordinal::one().pow(5), Ordinal::one());
    }

    #[test]
    fn checked_factorial_i128_has_the_grundy_roof() {
        assert_eq!(checked_factorial_i128(-1), None);
        assert_eq!(checked_factorial_i128(0), Some(1));
        assert_eq!(checked_factorial_i128(5), Some(120));
        assert!(checked_factorial_i128(33).is_some());
        assert_eq!(checked_factorial_i128(34), None);
    }

    #[test]
    fn factorial_in_scalar_uses_the_world_ring_map() {
        assert_eq!(factorial_in_scalar::<Integer>(5), Some(Integer(120)));
        assert_eq!(factorial_in_scalar::<Fp<7>>(6), Some(Fp::<7>::from_int(-1)));
        assert_eq!(factorial_in_scalar::<Fp<7>>(7), Some(Fp::<7>::zero()));
        assert_eq!(factorial_in_scalar::<Nimber>(4), Some(Nimber::zero()));
    }

    #[test]
    fn modular_helpers_cover_full_u128_range() {
        let m = (5u128).pow(55);
        assert!(m > i128::MAX as u128);
        assert_eq!(add_mod_u128(m - 1, m - 1, m), m - 2);
        assert_eq!(sub_mod_u128(1, 2, m), m - 1);
        assert_eq!(mul_mod_u128(m - 1, m - 1, m), 1);
        let inv = mod_inverse_u128(2, m).expect("2 is a unit modulo 5^55");
        assert_eq!(mul_mod_u128(2, inv, m), 1);
        assert_eq!(reduce_i128_mod_u128(-2, m), m - 2);
    }
}