spirix 0.1.1

use crate::core::integer::*;
use crate::core::undefined::*;
use crate::{Integer, Scalar, ScalarConstants};
use core::ops::*;
use i256::I256;
use num_traits::{AsPrimitive, WrappingAdd, WrappingMul, WrappingNeg, WrappingSub};

#[allow(private_bounds)]
impl<
        F: Integer
            + FullInt
            + Shl<isize, Output = F>
            + Shr<isize, Output = F>
            + Shl<F, Output = F>
            + Shr<F, Output = F>
            + Shl<E, Output = F>
            + Shr<E, Output = F>
            + WrappingNeg
            + WrappingAdd
            + WrappingMul
            + WrappingSub,
        E: Integer
            + FullInt
            + Shl<isize, Output = E>
            + Shr<isize, Output = E>
            + Shl<E, Output = E>
            + Shr<E, Output = E>
            + Shl<F, Output = E>
            + Shr<F, Output = E>
            + WrappingNeg
            + WrappingAdd
            + WrappingMul
            + WrappingSub,
    > Scalar<F, E>
where
    Scalar<F, E>: ScalarConstants,
    u8: AsPrimitive<F>,
    u16: AsPrimitive<F>,
    u32: AsPrimitive<F>,
    u64: AsPrimitive<F>,
    u128: AsPrimitive<F>,
    usize: AsPrimitive<F>,
    i8: AsPrimitive<F>,
    i16: AsPrimitive<F>,
    i32: AsPrimitive<F>,
    i64: AsPrimitive<F>,
    i128: AsPrimitive<F>,
    isize: AsPrimitive<F>,
    I256: From<F>,
    u8: AsPrimitive<E>,
    u16: AsPrimitive<E>,
    u32: AsPrimitive<E>,
    u64: AsPrimitive<E>,
    u128: AsPrimitive<E>,
    usize: AsPrimitive<E>,
    i8: AsPrimitive<E>,
    i16: AsPrimitive<E>,
    i32: AsPrimitive<E>,
    i64: AsPrimitive<E>,
    i128: AsPrimitive<E>,
    isize: AsPrimitive<E>,
    I256: From<E>,
{
    /// Calculates the mathematical modulus of this Scalar
    ///
    /// # Description
    ///
    /// Returns the unique value `r` such that `numerator ≡ r (mod period)`, where `r` is the canonical representative of the numerator's congruence class with sign following the period.
    ///
    /// `⬆` denotes *transfinite* (either exploded `[↑]` or infinite `[∞]`). Undefined prefixes collapse both into one tag since the distinction isn't preserved in the stored undefined class.
    ///
    /// # Special Cases (first matching rule wins, in order):
    ///
    /// 1. `[℘?]` numerator or period → the first undefined encountered
    /// 2. `[0]` numerator or period → `[0]`
    /// 3. Transfinite numerator (`[↑]` or `[∞]`):
    ///    - transfinite period → `[℘⬆%⬆]`
    ///    - otherwise → `[℘⬆%]`
    /// 4. Vanished period `[↓]`:
    ///    - vanished numerator → `[℘↓%↓]`
    ///    - otherwise → `[℘%↓]`
    /// 5. Infinite period `[∞]` → `[℘%⬆]` (signless period, no sign to floor against)
    /// 6. Exploded period `[↑]`:
    ///    - signs match → numerator (preserved)
    ///    - signs differ → `[℘%⬆]`
    /// 7. Both finite (normal or vanished with normal period) → integer floored modulus
    ///
    /// # Returns
    ///
    /// - `[#] % [#]` ➔ `[0]`, `[↓]`, or `[#]`
    /// - `[#] % [↓]` ➔ `[℘%↓]`
    /// - `[↓] % [↓]` ➔ `[℘↓%↓]`
    /// - `[#] % [↑]` signs match ➔ `[#]`
    /// - `[#] % [↑]` signs differ ➔ `[℘%⬆]`
    /// - `[↓] % [↑]` signs match ➔ `[↓]`
    /// - `[↓] % [↑]` signs differ ➔ `[℘%⬆]`
    /// - `[?] % [∞]` ➔ `[℘%⬆]` (signless period)
    /// - `[↑] % [#]` or `[↑] % [↓]` or `[∞] % [#]` or `[∞] % [↓]` ➔ `[℘⬆%]`
    /// - `[↑] % [↑]` or `[↑] % [∞]` or `[∞] % [↑]` or `[∞] % [∞]` ➔ `[℘⬆%⬆]`
    /// - `[0] % [?]` or `[?] % [0]` ➔ `[0]`
    ///
    /// # Examples
    ///
    /// ```rust
    /// use spirix::{Scalar, ScalarF5E3};
    ///
    /// // Basic modulus let a = ScalarF5E3::from(7_i32); let b = ScalarF5E3::from(3_i32); assert!(a % b == 1_i32);  // 7 % 3 = 1
    ///
    /// // Result takes sign of period let neg_a = ScalarF5E3::from(-7_i32); assert!(neg_a % b == 2_i32);   // -7 % 3 = 2 let neg_b = ScalarF5E3::from(-3_i32); assert!(a % neg_b == -2_i32);  // 7 % -3 = -2
    ///
    /// // Exploded period, matching signs: numerator preserved let exploded: ScalarF5E3 = ScalarF5E3::MAX * 2_i32; assert!((ScalarF5E3::PI % exploded) == ScalarF5E3::PI);
    ///
    /// // Exploded period, differing signs: undefined assert!((ScalarF5E3::PI % -exploded).is_undefined());
    ///
    /// // Vanished period: undefined let vanished: ScalarF5E3 = ScalarF5E3::MIN_POS / 19_i32; assert!((ScalarF5E3::from(42_i32) % vanished).is_undefined());
    ///
    /// // Zero cases assert!((ScalarF5E3::ZERO % ScalarF5E3::PI).is_zero()); assert!((ScalarF5E3::PI % ScalarF5E3::ZERO).is_zero());
    /// ```
    pub(crate) fn scalar_modulus_scalar(&self, modulus: &Scalar<F, E>) -> Scalar<F, E> {
        if !self.is_normal() || !modulus.is_normal() {
            if self.is_undefined() {
                return *self;
            }
            if modulus.is_undefined() {
                return *modulus;
            }
            if self.is_zero() || modulus.is_zero() {
                return Self::ZERO;
            }
            // Rule 3: Transfinite numerator ([↑] or [∞])
            if self.is_transfinite() {
                if modulus.is_transfinite() {
                    return Self {
                        fraction: TRANSFINITE_MODULUS_TRANSFINITE.prefix.sa(),
                        exponent: Self::ambiguous_exponent(),
                    };
                }
                return Self {
                    fraction: TRANSFINITE_MODULUS.prefix.sa(),
                    exponent: Self::ambiguous_exponent(),
                };
            }
            // Rule 4: Vanished period
            if modulus.vanished() {
                if self.vanished() {
                    return Self {
                        fraction: VANISHED_MODULUS_VANISHED.prefix.sa(),
                        exponent: Self::ambiguous_exponent(),
                    };
                }
                return Self {
                    fraction: MODULUS_VANISHED.prefix.sa(),
                    exponent: Self::ambiguous_exponent(),
                };
            }
            // Rule 5: Infinite period (signless — no sign to floor against)
            if modulus.is_infinite() {
                return Self {
                    fraction: MODULUS_TRANSFINITE.prefix.sa(),
                    exponent: Self::ambiguous_exponent(),
                };
            }
            // Rule 6: Exploded period
            if modulus.exploded() {
                if self.is_negative() == modulus.is_negative() {
                    return *self;
                }
                if self.vanished() {
                    return *modulus;
                }
                return Self {
                    fraction: MODULUS_TRANSFINITE.prefix.sa(),
                    exponent: Self::ambiguous_exponent(),
                };
            }
            // Vanished numerator with normal period: |↓| < |#| always
            if self.vanished() {
                if self.is_negative() == modulus.is_negative() {
                    return *self;
                }
                return *modulus;
            }
            return *self;
        }

        // Both operands are normal. Compute floored modulus using proper restoring-divider style remainder: align fractions by exponent, do integer modulo on inflated wide values, apply floored sign rule.
        //
        // Floored mod: result has the sign of the divisor. same signs:    result = a_mag mod b_mag, signed like a/b diff signs:    result = b_mag - (a_mag mod b_mag), signed like b
        //
        // |a| < |b| short-circuit:
        //   same signs:    result = a (already in [0, b) magnitude) diff signs:    result = a + b (one b-step over to land on b's side)

        let a_neg = self.is_negative();
        let b_neg = modulus.is_negative();
        let signs_differ = a_neg != b_neg;

        // |a| < |b|: shortcut with sign correction. Unsigned compare on AMBIG=0 cycle position gives the correct cyclic-magnitude ordering.
        if self.exponent.into_unsigned() < modulus.exponent.into_unsigned() {
            if !signs_differ {
                // Same sign: a is already the remainder.
                return *self;
            }
            // Diff signs: result = a + b. Inlined add (skips abnormal checks — we already know both are normal — and the result can't underflow to zero: |a+b| = |b|-|a| > 0 since |a|<|b|).
            let exp_diff_ba = modulus.exponent.wrapping_sub(&self.exponent);
            // Overflow in E (true diff > E::MAX): a utterly negligible vs b → result ≈ b.
            if exp_diff_ba.is_negative() {
                return *modulus;
            }
            let shift_ba: isize = exp_diff_ba.saturate();
            if shift_ba >= Self::fraction_bits() {
                // a is negligible at b's precision: result = b.
                return *modulus;
            }
            let mut big_f = modulus.fraction.inflate(true);
            big_f.w_shl_assign(shift_ba);
            let sum = big_f.w_add(self.fraction.inflate(true));
            let leading = sum.leading_same();
            let delta_i: isize = Self::fraction_bits().wrapping_sub(leading);
            let delta_e: E = delta_i.as_();
            let offset = self.exponent.wrapping_add(&delta_e);
            // AMBIG=0 cycle-position underflow check (mirrors scalar bitwise_normal): result wraps past the AMBIG sentinel when delta is negative and `offset - 1` lands ≥ self.exp in cycle-position space. The earlier `mod.exp.is_negative()` test predated AMBIG=0 and now means "|mod| ≥ 1" which is unrelated to underflow.
            let one_e: E = 1u8.as_();
            let underflowed = (delta_i.is_negative()
                && offset.wrapping_sub(&one_e).into_unsigned() >= self.exponent.into_unsigned())
                || offset == Self::ambiguous_exponent();
            // Folded net shift (shl(L).shr(fb) as one signed shift) to sidestep Rust's shift-overflow semantics when L == wide_bits.
            let fb = Self::fraction_bits();
            let shl_amount = if underflowed {
                leading.wrapping_sub(2).wrapping_sub(fb)
            } else {
                leading.wrapping_sub(fb)
            };
            let canonical = if shl_amount >= 0 {
                sum.w_shl(shl_amount)
            } else {
                sum.w_shr(shl_amount.wrapping_neg())
            };
            if underflowed {
                return Self {
                    fraction: canonical.deflate(),
                    exponent: Self::ambiguous_exponent(),
                };
            }
            return Self {
                fraction: canonical.deflate(),
                exponent: offset,
            };
        }

        let exp_diff_e = self.exponent.wrapping_sub(&modulus.exponent);
        let exp_diff: isize = exp_diff_e.saturate();

        // Inflate both fractions to wide effective values, take magnitudes.
        let a_wide = self.fraction.inflate(true);
        let b_wide = modulus.fraction.inflate(true);
        let a_mag = if a_neg { a_wide.w_neg() } else { a_wide };
        let b_mag = if b_neg { b_wide.w_neg() } else { b_wide };

        // Compute (a_mag << exp_diff) mod b_mag. exp_diff can be astronomically large for wide exponent types (up to ~2^127 at E7), so linear chunk-shifting is unusable — the old loop needed exp_diff/FRAC iterations and effectively hung. Instead use (a · (2^exp_diff mod b)) mod b with square-and-multiply over exp_diff's BITS: O(E_BITS) modular steps.
        // modmul is a shift-and-conditional-subtract peasant multiply: every intermediate stays < 2·b_mag, well inside the wide width, and reads the operands unsigned-positive (all values here are magnitudes < 2^(FRAC+1)).
        let _ = exp_diff; // superseded by bit-iteration over exp_diff_e (the saturate() could also flip sign for large diffs)
        let fb = Self::fraction_bits();
        let modmul = |x: F::Wide, y: F::Wide, m: F::Wide| -> F::Wide {
            let mut acc = F::zero().inflate(false);
            // y < 2^(FRAC+1): walk its bits MSB-first.
            let bits = fb.wrapping_add(2);
            let mut k = bits;
            while k > 0 {
                k -= 1;
                acc = acc.w_shl(1);
                if acc >= m {
                    acc = acc.w_sub(m);
                }
                let one_w = F::one().inflate(false);
                if y.w_shr_logical(k).w_and(one_w) == one_w {
                    acc = acc.w_add(x);
                    if acc >= m {
                        acc = acc.w_sub(m);
                    }
                }
            }
            acc
        };
        let mut rem = a_mag.w_rem_unsigned(b_mag);
        // pow2 = 2^exp_diff mod b_mag via square-and-multiply on the UNSIGNED cycle-position difference (bit k of the signed E value is bit k of the unsigned one, so raw-bit iteration is exact even where saturate() would have flipped sign).
        let mut pow2 = F::one().inflate(false).w_shl(1); // 2
        if pow2 >= b_mag {
            pow2 = pow2.w_sub(b_mag);
        }
        let mut acc_pow = F::one().inflate(false); // 1 (b_mag ≥ 2 for any normalized divisor)
        let e_bits = Self::exponent_bits();
        let mut k: isize = 0;
        while k < e_bits {
            let one_e: E = 1u8.as_();
            if (exp_diff_e >> k) & one_e == one_e {
                acc_pow = modmul(acc_pow, pow2, b_mag);
            }
            k += 1;
            if k < e_bits {
                pow2 = modmul(pow2, pow2, b_mag);
            }
        }
        rem = modmul(rem, acc_pow, b_mag);
        let r_mag = rem;

        if r_mag.w_is_zero() {
            return Self::ZERO;
        }

        // Floored sign correction.
        let result_wide = if signs_differ {
            // result magnitude = |b| - r_mag, sign of b
            let mag = b_mag.w_sub(r_mag);
            if b_neg {
                mag.w_neg()
            } else {
                mag
            }
        } else {
            // result magnitude = r_mag, sign of b (= sign of a)
            if b_neg {
                r_mag.w_neg()
            } else {
                r_mag
            }
        };

        if result_wide.w_is_zero() {
            return Self::ZERO;
        }

        // Normalize result at b's exponent. Same pattern as scalar_add_scalar.
        let leading = result_wide.leading_same();
        let delta_i: isize = fb.wrapping_sub(leading);
        let delta_e: E = delta_i.as_();
        let offset = modulus.exponent.wrapping_add(&delta_e);

        // AMBIG=0 cycle-position underflow check (mirrors scalar bitwise_normal): result wraps past the AMBIG sentinel when delta is negative and `offset - 1` lands ≥ modulus.exp in cycle-position space. The earlier `mod.exp.is_negative()` test predated AMBIG=0.
        let one_e: E = 1u8.as_();
        let underflowed = (delta_i.is_negative()
            && offset.wrapping_sub(&one_e).into_unsigned() >= modulus.exponent.into_unsigned())
            || offset == Self::ambiguous_exponent();
        // Folded net shift (shl(L).shr(fb) as one signed shift) to sidestep Rust's shift-overflow semantics when L == wide_bits.
        let shl_amount = if underflowed {
            leading.wrapping_sub(2).wrapping_sub(fb)
        } else {
            leading.wrapping_sub(fb)
        };
        let canonical = if shl_amount >= 0 {
            result_wide.w_shl(shl_amount)
        } else {
            result_wide.w_shr(shl_amount.wrapping_neg())
        };
        if underflowed {
            return Self {
                fraction: canonical.deflate(),
                exponent: Self::ambiguous_exponent(),
            };
        }
        Self {
            fraction: canonical.deflate(),
            exponent: offset,
        }
    }
}