ffnt 0.11.0

// TODO: code duplication with z64
#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};

use std::{
    fmt::{self, Display},
    iter::{Product, Sum},
    num::{IntErrorKind, TryFromIntError},
    ops::{
        Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign,
    },
    str::FromStr,
};

use crate::{ParseIntError, Z64, z64::TryDiv};

/// Element of a finite field with a 32 bit characteristic `P`
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
#[derive(Copy, Clone, Debug, Default, Eq, PartialEq, Ord, PartialOrd, Hash)]
#[repr(transparent)]
pub struct Z32<const P: u32>(u32);

impl<const P: u32> Z32<P> {
    const INFO: Z32Info = Z32Info::new(P);

    /// Minimum field element, i.e. 0
    pub const MIN: Z32<P> = {
        assert!(P > 0);
        Self::new_unchecked(0)
    };
    /// Maximum field element, i.e. `P - 1`
    pub const MAX: Z32<P> = {
        assert!(P > 1);
        Self::new_unchecked(P - 1)
    };

    /// Create a new field element corresponding to some integer
    ///
    /// The integer is reduced modulo the field characteristic `P`
    pub const fn new(z: i32) -> Self {
        let res = remi(z, P, Self::info().red_struct);
        debug_assert!(res >= 0);
        let res = res as u32;
        debug_assert!(res < P);
        Self::new_unchecked(res)
    }

    /// Create a new field element corresponding to some integer
    /// without modular reduction
    ///
    /// # Safety
    ///
    /// The argument should be less than `P`
    pub const fn new_unchecked(z: u32) -> Self {
        assert!(P > 0);
        debug_assert!(z < P);
        Self(z)
    }

    /// The multiplicative inverse `1/z` of a field element `z`.
    ///
    /// # Panics
    ///
    /// Panics if `z` is not invertible. If the characteristic `P` is
    /// a prime power this happens only if `z` is zero.
    pub const fn inv(&self) -> Self {
        self.try_inv()
            .expect("Number has no multiplicative inverse")
    }

    /// The multiplicative inverse `1/z` of a field element `z` or
    /// `None` if the inverse does not exist
    pub const fn try_inv(&self) -> Option<Self> {
        let res = extended_gcd(self.0, Self::modulus());
        if res.gcd != 1 {
            return None;
        }
        let s = res.bezout[0];
        let inv = if s < 0 {
            debug_assert!(s + Self::modulus() as i32 >= 0);
            s + Self::modulus() as i32
        } else {
            s
        } as u32;
        let inv = Self::new_unchecked(inv);
        Some(inv)
    }

    /// Check if a field element `z` has a multiplicative inverse `1/z`
    ///
    /// If you know that the characteristic is a prime it is usually
    /// better to check if `z` is zero.
    pub const fn has_inv(&self) -> bool {
        gcd(self.0, Self::modulus()) == 1
    }

    const fn info() -> &'static Z32Info {
        &Self::INFO
    }

    /// The field characteristic `P`
    pub const fn modulus() -> u32 {
        P
    }

    #[allow(missing_docs)]
    pub const fn modulus_inv() -> SpInverse32 {
        Self::info().p_inv
    }

    /// `z` to some integer power `exp`
    pub fn powi(self, exp: i64) -> Self {
        if exp < 0 {
            self.powu((-exp) as u64).inv()
        } else {
            self.powu(exp as u64)
        }
    }

    /// `z` to some integer power `exp`
    pub fn powu(mut self, mut exp: u64) -> Self {
        let mut res = Self::new_unchecked(1);
        while exp > 0 {
            if exp & 1 != 0 {
                res *= self
            };
            self *= self;
            exp /= 2;
        }
        res
    }

    #[cfg(any(feature = "rand", feature = "num-traits"))]
    pub(crate) const fn repr(self) -> u32 {
        self.0
    }
}

impl<const P: u32, const Q: u64> From<Z64<Q>> for Z32<P> {
    fn from(z: Z64<Q>) -> Self {
        u64::from(z).into()
    }
}

impl<const P: u32> From<Z32<P>> for u128 {
    fn from(i: Z32<P>) -> Self {
        i.0 as _
    }
}

impl<const P: u32> From<Z32<P>> for i128 {
    fn from(i: Z32<P>) -> Self {
        i.0 as _
    }
}

impl<const P: u32> From<Z32<P>> for u64 {
    fn from(i: Z32<P>) -> Self {
        i.0 as _
    }
}

impl<const P: u32> From<Z32<P>> for i64 {
    fn from(i: Z32<P>) -> Self {
        i.0 as _
    }
}

impl<const P: u32> From<Z32<P>> for u32 {
    fn from(i: Z32<P>) -> Self {
        i.0
    }
}

impl<const P: u32> From<Z32<P>> for i32 {
    fn from(i: Z32<P>) -> Self {
        i.0 as i32
    }
}

impl<const P: u32> TryFrom<Z32<P>> for u16 {
    type Error = TryFromIntError;

    fn try_from(i: Z32<P>) -> Result<Self, Self::Error> {
        i.0.try_into()
    }
}

impl<const P: u32> TryFrom<Z32<P>> for i16 {
    type Error = TryFromIntError;

    fn try_from(i: Z32<P>) -> Result<Self, Self::Error> {
        i.0.try_into()
    }
}

impl<const P: u32> TryFrom<Z32<P>> for u8 {
    type Error = TryFromIntError;

    fn try_from(i: Z32<P>) -> Result<Self, Self::Error> {
        i.0.try_into()
    }
}

impl<const P: u32> TryFrom<Z32<P>> for i8 {
    type Error = TryFromIntError;

    fn try_from(i: Z32<P>) -> Result<Self, Self::Error> {
        i.0.try_into()
    }
}

impl<const P: u32> From<u128> for Z32<P> {
    fn from(u: u128) -> Self {
        (u.rem_euclid(P as u128) as u32).into()
    }
}

impl<const P: u32> From<i128> for Z32<P> {
    fn from(i: i128) -> Self {
        (i.rem_euclid(P as i128) as u32).into()
    }
}

impl<const P: u32> From<u64> for Z32<P> {
    fn from(u: u64) -> Self {
        (u.rem_euclid(P as u64) as u32).into()
    }
}

impl<const P: u32> From<i64> for Z32<P> {
    fn from(i: i64) -> Self {
        (i.rem_euclid(P as i64) as u32).into()
    }
}

impl<const P: u32> From<u32> for Z32<P> {
    fn from(u: u32) -> Self {
        let num = remu(u, Self::modulus(), Self::info().red_struct) as u32;
        Self::new_unchecked(num)
    }
}

impl<const P: u32> From<i32> for Z32<P> {
    fn from(i: i32) -> Self {
        Self::new(i)
    }
}

impl<const P: u32> From<i16> for Z32<P> {
    fn from(i: i16) -> Self {
        Self::from(i as i32)
    }
}

impl<const P: u32> From<u16> for Z32<P> {
    fn from(u: u16) -> Self {
        Self::from(u as u32)
    }
}

impl<const P: u32> From<i8> for Z32<P> {
    fn from(i: i8) -> Self {
        Self::from(i as i32)
    }
}

impl<const P: u32> From<u8> for Z32<P> {
    fn from(u: u8) -> Self {
        Self::from(u as u32)
    }
}

impl<'a, const P: u32> TryFrom<&'a str> for Z32<P> {
    type Error = ParseIntError;

    fn try_from(s: &'a str) -> Result<Self, Self::Error> {
        s.parse()
    }
}

impl<const P: u32> FromStr for Z32<P> {
    type Err = ParseIntError;

    fn from_str(s: &str) -> Result<Self, Self::Err> {
        let z = s.parse()?;
        if z >= P {
            return Err(IntErrorKind::PosOverflow.into());
        }
        // # Safety
        // we just checked that z < P
        Ok(Self::new_unchecked(z))
    }
}

impl<const P: u32> Display for Z32<P> {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        write!(f, "{}", self.0)
    }
}

impl<const P: u32> AddAssign for Z32<P> {
    fn add_assign(&mut self, rhs: Self) {
        *self = *self + rhs;
    }
}

impl<const P: u32> AddAssign<&Z32<P>> for Z32<P> {
    fn add_assign(&mut self, rhs: &Self) {
        *self = *self + *rhs;
    }
}

impl<const P: u32> SubAssign for Z32<P> {
    fn sub_assign(&mut self, rhs: Self) {
        *self = *self - rhs;
    }
}

impl<const P: u32> SubAssign<&Z32<P>> for Z32<P> {
    fn sub_assign(&mut self, rhs: &Self) {
        *self -= *rhs;
    }
}

impl<const P: u32> MulAssign for Z32<P> {
    fn mul_assign(&mut self, rhs: Self) {
        *self = *self * rhs;
    }
}

impl<const P: u32> MulAssign<&Z32<P>> for Z32<P> {
    fn mul_assign(&mut self, rhs: &Self) {
        *self = *self * *rhs;
    }
}

impl<const P: u32> DivAssign for Z32<P> {
    fn div_assign(&mut self, rhs: Self) {
        *self = *self / rhs;
    }
}

impl<const P: u32> DivAssign<&Z32<P>> for Z32<P> {
    fn div_assign(&mut self, rhs: &Self) {
        *self = *self / *rhs;
    }
}

impl<const P: u32> Add for Z32<P> {
    type Output = Self;

    fn add(self, rhs: Self) -> Self::Output {
        let res = correct_excess((self.0 + rhs.0) as i32, Self::modulus());
        debug_assert!(res >= 0);
        let res = res as u32;
        Self::new_unchecked(res)
    }
}

impl<const P: u32> Add for &Z32<P> {
    type Output = Z32<P>;

    fn add(self, rhs: Self) -> Self::Output {
        *self + *rhs
    }
}

impl<const P: u32> Add<Z32<P>> for &Z32<P> {
    type Output = Z32<P>;

    fn add(self, rhs: Z32<P>) -> Self::Output {
        *self + rhs
    }
}

impl<const P: u32> Add<&Z32<P>> for Z32<P> {
    type Output = Z32<P>;

    fn add(self, rhs: &Z32<P>) -> Self::Output {
        self + *rhs
    }
}

impl<const P: u32> Sub for Z32<P> {
    type Output = Self;

    fn sub(self, rhs: Self) -> Self::Output {
        let res =
            correct_deficit(self.0 as i32 - rhs.0 as i32, Self::modulus());
        debug_assert!(res >= 0);
        let res = res as u32;
        Self::new_unchecked(res)
    }
}

impl<const P: u32> Sub for &Z32<P> {
    type Output = Z32<P>;

    fn sub(self, rhs: Self) -> Self::Output {
        *self - *rhs
    }
}

impl<const P: u32> Sub<Z32<P>> for &Z32<P> {
    type Output = Z32<P>;

    fn sub(self, rhs: Z32<P>) -> Self::Output {
        *self - rhs
    }
}

impl<const P: u32> Sub<&Z32<P>> for Z32<P> {
    type Output = Z32<P>;

    fn sub(self, rhs: &Z32<P>) -> Self::Output {
        self - *rhs
    }
}

impl<const P: u32> Neg for Z32<P> {
    type Output = Self;

    fn neg(self) -> Self::Output {
        Self::default() - self
    }
}

impl<const P: u32> Mul for Z32<P> {
    type Output = Self;

    fn mul(self, rhs: Self) -> Self::Output {
        let num = mul_mod(self.0, rhs.0, Self::modulus(), Self::modulus_inv());
        Self::new_unchecked(num)
    }
}

impl<const P: u32> Mul for &Z32<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: Self) -> Self::Output {
        *self * *rhs
    }
}

impl<const P: u32> Mul<Z32<P>> for &Z32<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: Z32<P>) -> Self::Output {
        *self * rhs
    }
}

impl<const P: u32> Mul<&Z32<P>> for Z32<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: &Z32<P>) -> Self::Output {
        self * *rhs
    }
}

impl<const P: u32> Div for Z32<P> {
    type Output = Self;

    #[allow(clippy::suspicious_arithmetic_impl)]
    fn div(self, rhs: Self) -> Self::Output {
        self * rhs.inv()
    }
}

const fn mul_mod(a: u32, b: u32, n: u32, ninv: SpInverse32) -> u32 {
    let res = normalised_mul_mod(
        a,
        (b as i32) << ninv.shamt,
        ((n as i32) << ninv.shamt) as u32,
        ninv.inv,
    ) >> ninv.shamt;
    res as u32
}

impl<const P: u32> Div for &Z32<P> {
    type Output = Z32<P>;

    fn div(self, rhs: Self) -> Self::Output {
        *self / *rhs
    }
}

impl<const P: u32> Div<Z32<P>> for &Z32<P> {
    type Output = Z32<P>;

    fn div(self, rhs: Z32<P>) -> Self::Output {
        *self / rhs
    }
}

impl<const P: u32> Div<&Z32<P>> for Z32<P> {
    type Output = Z32<P>;

    fn div(self, rhs: &Z32<P>) -> Self::Output {
        self / *rhs
    }
}

impl<const P: u32> TryDiv for Z32<P> {
    type Output = Self;

    fn try_div(self, rhs: Self) -> Option<Self::Output> {
        rhs.try_inv().map(|i| self * i)
    }
}

impl<const P: u32> TryDiv for &Z32<P> {
    type Output = Z32<P>;

    fn try_div(self, rhs: Self) -> Option<Self::Output> {
        (*self).try_div(*rhs)
    }
}

impl<const P: u32> TryDiv<Z32<P>> for &Z32<P> {
    type Output = Z32<P>;

    fn try_div(self, rhs: Z32<P>) -> Option<Self::Output> {
        (*self).try_div(rhs)
    }
}

impl<const P: u32> TryDiv<&Z32<P>> for Z32<P> {
    type Output = Z32<P>;

    fn try_div(self, rhs: &Z32<P>) -> Option<Self::Output> {
        self.try_div(*rhs)
    }
}

impl<const P: u32> Sum for Z32<P> {
    fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
        iter.fold(Self::new_unchecked(0), |a, b| a + b)
    }
}

impl<const P: u32> Product for Z32<P> {
    fn product<I: Iterator<Item = Self>>(iter: I) -> Self {
        iter.fold(Self::new_unchecked(1), |a, b| a * b)
    }
}

const fn normalised_mul_mod(a: u32, b: i32, n: u32, ninv: u32) -> i32 {
    let u = a as u64 * b as u64;
    let h = (u >> (SP_NBITS - 2)) as u32;
    let q = u64_mul_high(h, ninv) >> POST_SHIFT;
    let l = u as u32;
    let r = l.wrapping_sub(q.wrapping_mul(n));
    debug_assert!(r < 2 * n);
    correct_excess(r as i32, n)
}

const fn remu(z: u32, p: u32, red: ReduceStruct) -> i32 {
    let q = u64_mul_high(z, red.ninv);
    let r = (z - q.wrapping_mul(p)) as i32;
    correct_excess(r, p)
}

const fn remi(z: i32, p: u32, red: ReduceStruct) -> i32 {
    let zu = (z as u32) & ((1u32 << (u32::BITS - 1)) - 1);
    let r = remu(zu, p, red);
    let s = i32_sign_mask(z) & (red.sgn as i32);
    correct_deficit(r - s, p)
}

const fn u64_mul_high(a: u32, b: u32) -> u32 {
    u64_get_high(a as u64 * b as u64)
}

const fn u64_get_high(u: u64) -> u32 {
    (u >> u32::BITS) as u32
}

const fn correct_excess(a: i32, p: u32) -> i32 {
    let n = p as i32;
    (a - n) + (i32_sign_mask(a - n) & n)
}

const fn correct_deficit(a: i32, p: u32) -> i32 {
    a + (i32_sign_mask(a) & (p as i32))
}

#[derive(Copy, Clone, Debug, Default, Eq, PartialEq, Ord, PartialOrd, Hash)]
struct ExtendedGCDResult {
    gcd: u32,
    bezout: [i32; 2],
}

const fn extended_gcd(a: u32, b: u32) -> ExtendedGCDResult {
    let mut old_r = a;
    let mut r = b;
    let mut old_s = 1;
    let mut s = 0;
    let mut old_t = 0;
    let mut t = 1;

    while r != 0 {
        let quotient = old_r / r;
        (old_r, r) = (r, old_r - quotient * r);
        (old_s, s) = (s, old_s - quotient as i32 * s);
        (old_t, t) = (t, old_t - quotient as i32 * t);
    }
    ExtendedGCDResult {
        gcd: old_r,
        bezout: [old_s, old_t],
    }
}

const fn gcd(mut a: u32, mut b: u32) -> u32 {
    while b != 0 {
        (a, b) = (b, a % b)
    }
    a
}

const SP_NBITS: u32 = u32::BITS - 2;
const PRE_SHIFT2: u32 = 2 * SP_NBITS + 1;
const POST_SHIFT: u32 = 1;

const fn used_bits(z: u32) -> u32 {
    u32::BITS - z.leading_zeros()
}

#[derive(Copy, Clone, Debug, Default, Eq, PartialEq, Ord, PartialOrd, Hash)]
struct Z32Info {
    p: u32,
    p_inv: SpInverse32,
    red_struct: ReduceStruct,
}

impl Z32Info {
    const fn new(p: u32) -> Self {
        assert!(p > 1);
        assert!(used_bits(p) <= SP_NBITS);

        let p_inv = prep_mul_mod(p);
        let red_struct = prep_rem(p);
        Self {
            p,
            p_inv,
            red_struct,
        }
    }
}

const fn prep_mul_mod(p: u32) -> SpInverse32 {
    let shamt = p.leading_zeros() - (u32::BITS - SP_NBITS);
    let inv = normalised_prep_mul_mod(p << shamt);
    SpInverse32 { inv, shamt }
}

#[derive(Copy, Clone, Debug, Default, Eq, PartialEq, Ord, PartialOrd, Hash)]
struct ReduceStruct {
    ninv: u32,
    sgn: u32,
}

const fn prep_rem(p: u32) -> ReduceStruct {
    let mut q = (1 << (u32::BITS - 1)) / p;
    // r = 2^31 % p
    let r = (1 << (u32::BITS - 1)) - q * p;

    q *= 2;
    q += correct_excess_quo(2 * r as i32, p as i32).0;

    ReduceStruct { ninv: q, sgn: r }
}

const fn correct_excess_quo(a: i32, n: i32) -> (u32, i32) {
    if a >= n { (1, a - n) } else { (0, a) }
}

const fn i32_sign_mask(i: i32) -> i32 {
    i >> (u32::BITS - 1)
}

const fn u32_sign_mask(i: u32) -> i32 {
    i32_sign_mask(i as i32)
}

#[allow(missing_docs)]
#[derive(Copy, Clone, Debug, Default, Eq, PartialEq, Ord, PartialOrd, Hash)]
pub struct SpInverse32 {
    inv: u32,
    shamt: u32,
}

// Adapted from NTL's sp_NormalizedPrepMulMod
//
// Floating-point arithmetic replaced be u64 / i64 to allow `const`.
// The performance impact is not a huge concern since this function
// is only evaluated at compile time and only once for each prime field order.
// This is unlike NTL, where each change triggers a recalculation?
//
// This only works since this function is `const` and can be therefore
// used to compute individual `const INFO` inside `Z32<P>` for each
// `P`. The alternatives `lazy_static!` or `OnceCell` would not be
// recomputed, but instead incorrectly shared between `Z32<P>` with
// different `P`!
const fn normalised_prep_mul_mod(n: u32) -> u32 {
    // NOTE: this is an initial approximation
    //       the true quotient is <= 2^SP_NBITS
    const MAX: u64 = 1u64 << (2 * SP_NBITS - 1);
    let init_quot_approx = MAX / n as u64;

    let approx_rem = MAX - n as u64 * init_quot_approx;

    let approx_rem = (approx_rem << (PRE_SHIFT2 - 2 * SP_NBITS + 1)) - 1;

    let approx_rem_low = approx_rem as u32;
    let s1 = (approx_rem >> u32::BITS) as u32;
    let s2 = approx_rem_low >> (u32::BITS - 1);
    let approx_rem_high = s1.wrapping_add(s2);

    let approx_rem_low = approx_rem_low as i32;
    let approx_rem_high = approx_rem_high as i32;

    let bpl = 1i64 << u32::BITS;

    let fr = approx_rem_low as i64 + approx_rem_high as i64 * bpl;

    // now convert fr*ninv to a long
    // but we have to be careful: fr may be negative.
    // the result should still give floor(r/n) pm 1,
    // and is computed in a way that avoids branching

    let mut q1 = (fr / n as i64) as i32;
    if q1 < 0 {
        // This counteracts the round-to-zero behavior of conversion
        // to i32.  It should be compiled into branch-free code.
        q1 -= 1
    }

    let mut q1 = q1 as u32;
    let approx_rem_low = approx_rem_low as u32;
    let sub = q1.wrapping_mul(n);

    let approx_rem = approx_rem_low.wrapping_sub(sub);

    q1 += (1
        + u32_sign_mask(approx_rem)
        + u32_sign_mask(approx_rem.wrapping_sub(n))) as u32;

    ((init_quot_approx as u32) << (PRE_SHIFT2 - 2 * SP_NBITS + 1))
        .wrapping_add(q1)

    // NTL_PRE_SHIFT1 is 0, so no further shift required
}

/// [Z32] variant for fast repeated multiplication
///
/// Use this variant to speed up repeated multiplication by the same value:
/// ```
/// # use ffnt::{Z32, z32::Z32FastMul};
/// const P: u32 = 10007;
/// let mut numbers = Vec::from_iter((1..1000).map(Z32::<P>::from));
///
/// let factor: Z32FastMul<P> = 12.into();
/// for number in &mut numbers {
///     // same as `number *= Z32::from(12)`, but faster
///     *number *= &factor;
/// }
/// ```
///
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
#[derive(Copy, Clone, Debug, Default, Eq, PartialEq, Ord, PartialOrd, Hash)]
pub struct Z32FastMul<const P: u32> {
    val: Z32<P>,
    val_over_mod_approx: u32,
}

impl<const P: u32> From<Z32<P>> for Z32FastMul<P> {
    fn from(val: Z32<P>) -> Self {
        let val_over_mod_approx = Self::prep_mul_mod_precon(val.0);
        Self {
            val,
            val_over_mod_approx,
        }
    }
}

impl<const P: u32> Z32FastMul<P> {
    fn prep_mul_mod_precon(val: u32) -> u32 {
        let p_inv = Z32::<P>::INFO.p_inv;
        normalized_prep_mul_mod_precon(
            val << p_inv.shamt,
            P << p_inv.shamt,
            p_inv.inv,
        ) << (u32::BITS - SP_NBITS)
    }
}

fn normalized_prep_mul_mod_precon(val: u32, p: u32, p_inv: u32) -> u32 {
    let h = val << 2;
    let q = u64_mul_high(h, p_inv);
    let q = q >> POST_SHIFT;
    let l = val << SP_NBITS;
    let r = l.wrapping_sub(q.wrapping_mul(p)); // r in [0..2*p)
    debug_assert!(r < 2 * p);

    q.saturating_add_signed(1 + i32_sign_mask(r as i32 - p as i32)) // NOTE: not shifted
}

impl<const P: u32> Mul<Z32FastMul<P>> for Z32<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: Z32FastMul<P>) -> Self::Output {
        let res = mul_mod_precon(self.0, rhs.val.0, P, rhs.val_over_mod_approx);
        Z32::new_unchecked(res as u32)
    }
}

impl<const P: u32> Mul<Z32<P>> for Z32FastMul<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: Z32<P>) -> Self::Output {
        rhs * self
    }
}

impl<const P: u32> Mul<Z32FastMul<P>> for &Z32<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: Z32FastMul<P>) -> Self::Output {
        *self * rhs
    }
}

impl<const P: u32> Mul<Z32<P>> for &Z32FastMul<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: Z32<P>) -> Self::Output {
        *self * rhs
    }
}

impl<const P: u32> Mul<&Z32FastMul<P>> for Z32<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: &Z32FastMul<P>) -> Self::Output {
        self * *rhs
    }
}

impl<const P: u32> Mul<&Z32<P>> for Z32FastMul<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: &Z32<P>) -> Self::Output {
        self * *rhs
    }
}

impl<'a, const P: u32> Mul<&'a Z32FastMul<P>> for &Z32<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: &'a Z32FastMul<P>) -> Self::Output {
        *self * *rhs
    }
}

impl<'a, const P: u32> Mul<&'a Z32<P>> for &Z32FastMul<P> {
    type Output = Z32<P>;

    fn mul(self, rhs: &'a Z32<P>) -> Self::Output {
        *self * *rhs
    }
}

impl<const P: u32> MulAssign<Z32FastMul<P>> for Z32<P> {
    fn mul_assign(&mut self, rhs: Z32FastMul<P>) {
        *self = *self * rhs
    }
}

impl<const P: u32> MulAssign<&Z32FastMul<P>> for Z32<P> {
    fn mul_assign(&mut self, rhs: &Z32FastMul<P>) {
        *self = *self * rhs
    }
}

fn mul_mod_precon(lhs: u32, rhs: u32, p: u32, rhs_over_mod_approx: u32) -> i32 {
    let q = u64_mul_high(lhs, rhs_over_mod_approx);
    let lhs_times_rhs = lhs.wrapping_mul(rhs);
    let q_times_p = q.wrapping_mul(p);
    let r = lhs_times_rhs.wrapping_sub(q_times_p);
    correct_excess(r as i32, p)
}

macro_rules! impl_fastmul_from {
    ( $( $t:ty ),* ) => {
        $(
            impl<const P: u32> From<$t> for Z32FastMul<P> {
                fn from(t: $t) -> Self {
                    Self::from(Z32::from(t))
                }
            }
        )*
    }
}

impl_fastmul_from!(i8, i16, i32, i64, i128, u8, u16, u32, u64, u128);

impl<const P: u32> From<Z32FastMul<P>> for Z32<P> {
    fn from(z: Z32FastMul<P>) -> Self {
        z.val
    }
}

#[cfg(test)]
mod tests {

    use ::rand::{Rng, SeedableRng};
    use once_cell::sync::Lazy;
    use rug::{Integer, ops::Pow};

    use super::*;

    const PRIMES: [u32; 3] = [3, 65521, 1073741789];

    #[test]
    fn z32_has_inv() {
        type Z = Z32<6>;
        assert!(!Z::from(0).has_inv());
        assert!(Z::from(1).has_inv());
        assert!(!Z::from(2).has_inv());
        assert!(!Z::from(3).has_inv());
        assert!(!Z::from(4).has_inv());
        assert!(Z::from(5).has_inv());
        assert_eq!(Z::from(6), Z::from(0));
    }

    #[test]
    #[should_panic]
    fn z32_inv0() {
        type Z = Z32<6>;
        Z::from(0).inv();
    }

    #[test]
    #[should_panic]
    fn z32_inv2() {
        type Z = Z32<6>;
        Z::from(2).inv();
    }

    #[test]
    fn z32_constr() {
        let z: Z32<3> = 2.into();
        assert_eq!(u32::from(z), 2);
        let z: Z32<3> = (-1).into();
        assert_eq!(u32::from(z), 2);
        let z: Z32<3> = 5.into();
        assert_eq!(u32::from(z), 2);

        let z: Z32<3> = 0.into();
        assert_eq!(u32::from(z), 0);
        let z: Z32<3> = 3.into();
        assert_eq!(u32::from(z), 0);

        let z: Z32<3> = 2u32.into();
        assert_eq!(u32::from(z), 2);
        let z: Z32<3> = 5u32.into();
        assert_eq!(u32::from(z), 2);

        let z: Z32<3> = 0u32.into();
        assert_eq!(u32::from(z), 0);
        let z: Z32<3> = 3u32.into();
        assert_eq!(u32::from(z), 0);
    }

    static POINTS: Lazy<[i32; 1000]> = Lazy::new(|| {
        let mut pts = [0; 1000];
        let mut rng = rand_xoshiro::Xoshiro256StarStar::seed_from_u64(0);
        for pt in &mut pts {
            *pt = rng.random();
        }
        pts
    });

    #[test]
    fn tst_conv() {
        for pt in *POINTS {
            let z: Z32<{ PRIMES[0] }> = pt.into();
            let z: i32 = z.into();
            assert_eq!(z, pt.rem_euclid(PRIMES[0] as i32));
        }

        for pt in *POINTS {
            let z: Z32<{ PRIMES[1] }> = pt.into();
            let z: i32 = z.into();
            assert_eq!(z, pt.rem_euclid(PRIMES[1] as i32));
        }

        for pt in *POINTS {
            let z: Z32<{ PRIMES[2] }> = pt.into();
            let z: i32 = z.into();
            assert_eq!(z, pt.rem_euclid(PRIMES[2] as i32));
        }
    }

    #[test]
    fn tst_add() {
        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[0] }> = pt1.into();
            let pt1 = pt1 as i64;
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[0] }> = pt2.into();
                let pt2 = pt2 as i64;
                let sum1: i32 = (z1 + z2).into();
                let sum2 = (pt1 + pt2).rem_euclid(PRIMES[0] as i64) as i32;
                assert_eq!(sum1, sum2);
            }
        }

        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[1] }> = pt1.into();
            let pt1 = pt1 as i64;
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[1] }> = pt2.into();
                let pt2 = pt2 as i64;
                let sum1: i32 = (z1 + z2).into();
                let sum2 = (pt1 + pt2).rem_euclid(PRIMES[1] as i64) as i32;
                assert_eq!(sum1, sum2);
            }
        }

        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[2] }> = pt1.into();
            let pt1 = pt1 as i64;
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[2] }> = pt2.into();
                let pt2 = pt2 as i64;
                let sum1: i32 = (z1 + z2).into();
                let sum2 = (pt1 + pt2).rem_euclid(PRIMES[2] as i64) as i32;
                assert_eq!(sum1, sum2);
            }
        }
    }

    #[test]
    fn tst_sub() {
        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[0] }> = pt1.into();
            let pt1 = pt1 as i64;
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[0] }> = pt2.into();
                let pt2 = pt2 as i64;
                let sum1: i32 = (z1 - z2).into();
                let sum2 = (pt1 - pt2).rem_euclid(PRIMES[0] as i64) as i32;
                assert_eq!(sum1, sum2);
            }
        }

        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[1] }> = pt1.into();
            let pt1 = pt1 as i64;
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[1] }> = pt2.into();
                let pt2 = pt2 as i64;
                let sum1: i32 = (z1 - z2).into();
                let sum2 = (pt1 - pt2).rem_euclid(PRIMES[1] as i64) as i32;
                assert_eq!(sum1, sum2);
            }
        }

        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[2] }> = pt1.into();
            let pt1 = pt1 as i64;
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[2] }> = pt2.into();
                let pt2 = pt2 as i64;
                let sum1: i32 = (z1 - z2).into();
                let sum2 = (pt1 - pt2).rem_euclid(PRIMES[2] as i64) as i32;
                assert_eq!(sum1, sum2);
            }
        }
    }

    #[test]
    fn tst_mul() {
        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[0] }> = pt1.into();
            let pt1 = pt1 as i64;
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[0] }> = pt2.into();
                let pt2 = pt2 as i64;
                let prod1: i32 = (z1 * z2).into();
                let prod2 = (pt1 * pt2).rem_euclid(PRIMES[0] as i64) as i32;
                assert_eq!(prod1, prod2);
            }
        }

        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[1] }> = pt1.into();
            let pt1 = pt1 as i64;
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[1] }> = pt2.into();
                let pt2 = pt2 as i64;
                let prod1: i32 = (z1 * z2).into();
                let prod2 = (pt1 * pt2).rem_euclid(PRIMES[1] as i64) as i32;
                assert_eq!(prod1, prod2);
            }
        }

        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[2] }> = pt1.into();
            let pt1 = pt1 as i64;
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[2] }> = pt2.into();
                let pt2 = pt2 as i64;
                let prod1: i32 = (z1 * z2).into();
                let prod2 = (pt1 * pt2).rem_euclid(PRIMES[2] as i64) as i32;
                assert_eq!(prod1, prod2);
            }
        }
    }

    #[test]
    fn tst_fastmul() {
        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[0] }> = pt1.into();
            let fast_z1 = Z32FastMul::from(z1);
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[0] }> = pt2.into();
                assert_eq!(z1 * z2, fast_z1 * z2);
            }
        }

        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[1] }> = pt1.into();
            let fast_z1 = Z32FastMul::from(z1);
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[1] }> = pt2.into();
                assert_eq!(z1 * z2, fast_z1 * z2);
            }
        }

        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[2] }> = pt1.into();
            let fast_z1 = Z32FastMul::from(z1);
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[2] }> = pt2.into();
                assert_eq!(z1 * z2, fast_z1 * z2);
            }
        }
    }

    #[test]
    fn tst_div() {
        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[0] }> = pt1.into();
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[0] }> = pt2.into();
                if i32::from(z2) == 0 {
                    continue;
                }
                let div = z1 / z2;
                assert_eq!(z1, div * z2);
            }
        }

        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[1] }> = pt1.into();
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[1] }> = pt2.into();
                if i32::from(z2) == 0 {
                    continue;
                }
                let div = z1 / z2;
                assert_eq!(z1, div * z2);
            }
        }

        for pt1 in *POINTS {
            let z1: Z32<{ PRIMES[2] }> = pt1.into();
            for pt2 in *POINTS {
                let z2: Z32<{ PRIMES[2] }> = pt2.into();
                if i32::from(z2) == 0 {
                    continue;
                }
                let div = z1 / z2;
                assert_eq!(z1, div * z2);
            }
        }
    }

    #[test]
    fn tst_pow() {
        let mut rng = rand_xoshiro::Xoshiro256StarStar::seed_from_u64(2849);
        for pt1 in *POINTS {
            let base = Integer::from(pt1);
            for _ in 0..100 {
                let exp: u8 = rng.random();
                let pow = base.clone().pow(exp as u32);
                // ensure remainder is positive and less than the mod
                let ref_pow0 =
                    (pow.clone() % PRIMES[0] + PRIMES[0]) % PRIMES[0];
                let ref_pow0: u32 = ref_pow0.try_into().unwrap();
                let z: Z32<{ PRIMES[0] }> = pt1.into();
                let pow0: u32 = z.powu(exp as u64).into();
                assert_eq!(pow0, ref_pow0);

                let ref_pow0 =
                    (pow.clone() % PRIMES[1] + PRIMES[1]) % PRIMES[1];
                let ref_pow0: u32 = ref_pow0.try_into().unwrap();
                let z: Z32<{ PRIMES[1] }> = pt1.into();
                let pow0: u32 = z.powu(exp as u64).into();
                assert_eq!(pow0, ref_pow0);

                let ref_pow0 = (pow % PRIMES[2] + PRIMES[2]) % PRIMES[2];
                let ref_pow0: u32 = ref_pow0.try_into().unwrap();
                let z: Z32<{ PRIMES[2] }> = pt1.into();
                let pow0: u32 = z.powu(exp as u64).into();
                assert_eq!(pow0, ref_pow0);
            }
        }
    }
}