nc_polynomial/

polynomial.rs

// MIT License
//
// Copyright (c) 2026 Raja Lehtihet & Wael El Oraiby
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//! Polynomial arithmetic over `Z_q[x]` with an explicit maximum degree bound.
//!
//! The library stores coefficients in dense form up to `max_degree`, always reduced modulo `q`.
//! It supports:
//! - basic arithmetic (`+`, `-`, negation, scalar multiply),
//! - polynomial multiplication (schoolbook and NTT),
//! - quotient-ring operations modulo another polynomial,
//! - evaluation, derivative, and long division with remainder.

use core::fmt;

/// Dense polynomial representation with fixed metadata (`max_degree`, `modulus`).
///
/// Internally, `coeffs` always has length `max_degree + 1`. Coefficients beyond the true degree
/// are stored as zeros.
#[derive(Debug, Clone, PartialEq, Eq)]
pub struct Polynomial {
    /// Dense coefficient vector where `coeffs[i]` is the coefficient of `x^i`.
    coeffs: Vec<u64>,
    /// Coefficient modulus `q` for arithmetic in `Z_q`.
    modulus: u64,
    /// Maximum supported degree `n`.
    max_degree: usize,
}

/// Errors returned by polynomial operations.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum PolynomialError {
    /// The modulus must be at least `2`.
    InvalidModulus(u64),
    /// Input coefficients requested a degree that exceeds `max_degree`.
    DegreeOverflow { requested: usize, max_degree: usize },
    /// Polynomials are incompatible (different degree bounds or moduli).
    IncompatiblePolynomials,
    /// Multiplication produced a degree higher than `max_degree`.
    ProductDegreeOverflow { degree: usize, max_degree: usize },
    /// NTT requires a non-zero power-of-two transform length.
    NttLengthMustBePowerOfTwo(usize),
    /// NTT length is not supported by the chosen modulus.
    NttLengthUnsupported { length: usize, modulus: u64 },
    /// Provided primitive root cannot generate a valid root of unity for the transform length.
    InvalidPrimitiveRoot {
        primitive_root: u64,
        length: usize,
        modulus: u64,
    },
    /// Division by the zero polynomial is undefined.
    DivisionByZeroPolynomial,
    /// A required multiplicative inverse does not exist under the modulus.
    NonInvertibleCoefficient { coefficient: u64, modulus: u64 },
}

impl fmt::Display for PolynomialError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::InvalidModulus(m) => write!(f, "invalid modulus {m}, expected q >= 2"),
            Self::DegreeOverflow {
                requested,
                max_degree,
            } => write!(
                f,
                "requested degree {requested} exceeds max_degree {max_degree}"
            ),
            Self::IncompatiblePolynomials => {
                write!(f, "polynomials are incompatible (different n or q)")
            }
            Self::ProductDegreeOverflow { degree, max_degree } => {
                write!(f, "product degree {degree} exceeds max_degree {max_degree}")
            }
            Self::NttLengthMustBePowerOfTwo(length) => {
                write!(f, "NTT length {length} must be a non-zero power of two")
            }
            Self::NttLengthUnsupported { length, modulus } => write!(
                f,
                "NTT length {length} is not supported by modulus {modulus} (length must divide q-1)"
            ),
            Self::InvalidPrimitiveRoot {
                primitive_root,
                length,
                modulus,
            } => write!(
                f,
                "primitive root {primitive_root} is invalid for NTT length {length} under modulus {modulus}"
            ),
            Self::DivisionByZeroPolynomial => write!(f, "division by zero polynomial"),
            Self::NonInvertibleCoefficient {
                coefficient,
                modulus,
            } => write!(
                f,
                "coefficient {coefficient} has no multiplicative inverse modulo {modulus}"
            ),
        }
    }
}

impl std::error::Error for PolynomialError {}

impl Polynomial {
    /// Creates a polynomial with degree bound `max_degree` over modulus `modulus`.
    ///
    /// Coefficients are normalized modulo `modulus` and then zero-padded to length
    /// `max_degree + 1`.
    ///
    /// # Errors
    /// Returns:
    /// - [`PolynomialError::InvalidModulus`] when `modulus < 2`,
    /// - [`PolynomialError::DegreeOverflow`] when `coeffs.len() > max_degree + 1`.
    pub fn new(max_degree: usize, modulus: u64, coeffs: &[u64]) -> Result<Self, PolynomialError> {
        if modulus < 2 {
            return Err(PolynomialError::InvalidModulus(modulus));
        }
        if coeffs.len() > max_degree + 1 {
            return Err(PolynomialError::DegreeOverflow {
                requested: coeffs.len().saturating_sub(1),
                max_degree,
            });
        }

        let mut normalized = vec![0_u64; max_degree + 1];
        for (i, value) in coeffs.iter().copied().enumerate() {
            normalized[i] = value % modulus;
        }

        Ok(Self {
            coeffs: normalized,
            modulus,
            max_degree,
        })
    }

    /// Returns the additive identity polynomial `0`.
    pub fn zero(max_degree: usize, modulus: u64) -> Result<Self, PolynomialError> {
        Self::new(max_degree, modulus, &[])
    }

    /// Returns the multiplicative identity polynomial `1`.
    pub fn one(max_degree: usize, modulus: u64) -> Result<Self, PolynomialError> {
        Self::new(max_degree, modulus, &[1])
    }

    /// Returns the configured maximum degree `n`.
    pub fn max_degree(&self) -> usize {
        self.max_degree
    }

    /// Returns the coefficient modulus `q`.
    pub fn modulus(&self) -> u64 {
        self.modulus
    }

    /// Returns the internal dense coefficient slice of length `max_degree + 1`.
    pub fn coefficients(&self) -> &[u64] {
        &self.coeffs
    }

    /// Returns coefficients trimmed to the actual degree.
    ///
    /// Zero polynomial is returned as `[0]`.
    pub fn trimmed_coefficients(&self) -> Vec<u64> {
        match self.degree() {
            Some(d) => self.coeffs[..=d].to_vec(),
            None => vec![0],
        }
    }

    /// Returns coefficient of `x^degree`, or `None` when out of bounds.
    pub fn coeff(&self, degree: usize) -> Option<u64> {
        self.coeffs.get(degree).copied()
    }

    /// Returns the true polynomial degree, or `None` for zero polynomial.
    pub fn degree(&self) -> Option<usize> {
        self.coeffs.iter().rposition(|&c| c != 0)
    }

    /// Returns `true` when all coefficients are zero.
    pub fn is_zero(&self) -> bool {
        self.degree().is_none()
    }

    /// Adds two polynomials coefficient-wise in `Z_q[x]`.
    ///
    /// # Errors
    /// Returns [`PolynomialError::IncompatiblePolynomials`] when `n` or `q` differ.
    pub fn add(&self, rhs: &Self) -> Result<Self, PolynomialError> {
        self.ensure_compatible(rhs)?;

        let mut out = vec![0_u64; self.max_degree + 1];
        for (i, slot) in out.iter_mut().enumerate() {
            *slot = mod_add(self.coeffs[i], rhs.coeffs[i], self.modulus);
        }

        Ok(Self {
            coeffs: out,
            modulus: self.modulus,
            max_degree: self.max_degree,
        })
    }

    /// Subtracts `rhs` from `self` coefficient-wise in `Z_q[x]`.
    ///
    /// # Errors
    /// Returns [`PolynomialError::IncompatiblePolynomials`] when `n` or `q` differ.
    pub fn sub(&self, rhs: &Self) -> Result<Self, PolynomialError> {
        self.ensure_compatible(rhs)?;

        let mut out = vec![0_u64; self.max_degree + 1];
        for (i, slot) in out.iter_mut().enumerate() {
            *slot = mod_sub(self.coeffs[i], rhs.coeffs[i], self.modulus);
        }

        Ok(Self {
            coeffs: out,
            modulus: self.modulus,
            max_degree: self.max_degree,
        })
    }

    /// Returns additive inverse `-self` in `Z_q[x]`.
    pub fn neg(&self) -> Self {
        let mut out = vec![0_u64; self.max_degree + 1];
        for (i, slot) in out.iter_mut().enumerate() {
            let c = self.coeffs[i];
            *slot = if c == 0 { 0 } else { self.modulus - c };
        }

        Self {
            coeffs: out,
            modulus: self.modulus,
            max_degree: self.max_degree,
        }
    }

    /// Multiplies all coefficients by `scalar` in `Z_q`.
    pub fn scalar_mul(&self, scalar: u64) -> Self {
        let mut out = vec![0_u64; self.max_degree + 1];
        let reduced_scalar = scalar % self.modulus;

        for (i, slot) in out.iter_mut().enumerate() {
            *slot = mod_mul(self.coeffs[i], reduced_scalar, self.modulus);
        }

        Self {
            coeffs: out,
            modulus: self.modulus,
            max_degree: self.max_degree,
        }
    }

    /// Schoolbook polynomial multiplication with strict degree bound enforcement.
    ///
    /// # Errors
    /// Returns:
    /// - [`PolynomialError::IncompatiblePolynomials`] when `n` or `q` differ,
    /// - [`PolynomialError::ProductDegreeOverflow`] if exact product degree exceeds `max_degree`.
    pub fn mul(&self, rhs: &Self) -> Result<Self, PolynomialError> {
        self.ensure_compatible(rhs)?;

        if let (Some(d1), Some(d2)) = (self.degree(), rhs.degree()) {
            let d = d1 + d2;
            if d > self.max_degree {
                return Err(PolynomialError::ProductDegreeOverflow {
                    degree: d,
                    max_degree: self.max_degree,
                });
            }
        }

        let mut out = vec![0_u64; self.max_degree + 1];
        for (i, &a) in self.coeffs.iter().enumerate() {
            if a == 0 {
                continue;
            }
            for (j, &b) in rhs.coeffs.iter().enumerate() {
                if b == 0 {
                    continue;
                }
                let idx = i + j;
                if idx > self.max_degree {
                    continue;
                }
                // Convolution accumulation: out[idx] += a_i * b_j (mod q).
                let prod = mod_mul(a, b, self.modulus);
                out[idx] = mod_add(out[idx], prod, self.modulus);
            }
        }

        Ok(Self {
            coeffs: out,
            modulus: self.modulus,
            max_degree: self.max_degree,
        })
    }

    /// NTT-based polynomial multiplication.
    ///
    /// This computes exact convolution via NTT and then places the result in the dense output
    /// buffer (length `max_degree + 1`).
    ///
    /// # Errors
    /// Returns:
    /// - [`PolynomialError::IncompatiblePolynomials`] when `n` or `q` differ,
    /// - [`PolynomialError::ProductDegreeOverflow`] if exact product degree exceeds `max_degree`,
    /// - NTT-specific errors when modulus/root parameters are not valid.
    pub fn mul_ntt(&self, rhs: &Self, primitive_root: u64) -> Result<Self, PolynomialError> {
        self.ensure_compatible(rhs)?;

        let (Some(lhs_degree), Some(rhs_degree)) = (self.degree(), rhs.degree()) else {
            // Zero times anything is zero.
            return Self::zero(self.max_degree, self.modulus);
        };

        let degree = lhs_degree + rhs_degree;
        if degree > self.max_degree {
            return Err(PolynomialError::ProductDegreeOverflow {
                degree,
                max_degree: self.max_degree,
            });
        }

        let lhs = &self.coeffs[..=lhs_degree];
        let rhs = &rhs.coeffs[..=rhs_degree];
        let product = convolution_ntt(lhs, rhs, self.modulus, primitive_root)?;

        let mut out = vec![0_u64; self.max_degree + 1];
        for (i, coeff) in product.into_iter().enumerate() {
            out[i] = coeff;
        }

        Ok(Self {
            coeffs: out,
            modulus: self.modulus,
            max_degree: self.max_degree,
        })
    }

    /// Schoolbook multiplication truncated to degree `max_degree`.
    ///
    /// Unlike [`Self::mul`], this method never returns degree-overflow errors; terms above
    /// `max_degree` are discarded.
    ///
    /// # Errors
    /// Returns [`PolynomialError::IncompatiblePolynomials`] when `n` or `q` differ.
    pub fn mul_truncated(&self, rhs: &Self) -> Result<Self, PolynomialError> {
        self.ensure_compatible(rhs)?;

        let mut out = vec![0_u64; self.max_degree + 1];
        for (i, &a) in self.coeffs.iter().enumerate() {
            if a == 0 {
                continue;
            }
            for (j, &b) in rhs.coeffs.iter().enumerate() {
                if b == 0 {
                    continue;
                }
                let idx = i + j;
                if idx > self.max_degree {
                    // Remaining j values only increase idx, so we can stop this inner loop.
                    break;
                }
                let prod = mod_mul(a, b, self.modulus);
                out[idx] = mod_add(out[idx], prod, self.modulus);
            }
        }

        Ok(Self {
            coeffs: out,
            modulus: self.modulus,
            max_degree: self.max_degree,
        })
    }

    /// Reduces `self` modulo `modulus_poly`, returning `self mod modulus_poly` in `Z_q[x]`.
    ///
    /// # Errors
    /// Returns:
    /// - [`PolynomialError::IncompatiblePolynomials`] when `n` or `q` differ,
    /// - [`PolynomialError::DivisionByZeroPolynomial`] if `modulus_poly` is zero,
    /// - [`PolynomialError::NonInvertibleCoefficient`] if leading coefficient of `modulus_poly`
    ///   is not invertible in `Z_q`.
    pub fn rem_mod_poly(&self, modulus_poly: &Self) -> Result<Self, PolynomialError> {
        self.ensure_compatible(modulus_poly)?;

        let mut reduced = self.coeffs.clone();
        reduce_coefficients_mod_poly(&mut reduced, &modulus_poly.coeffs, self.modulus)?;

        Ok(Self {
            coeffs: reduced,
            modulus: self.modulus,
            max_degree: self.max_degree,
        })
    }

    /// Adds two polynomials and reduces modulo `modulus_poly`.
    pub fn add_mod_poly(&self, rhs: &Self, modulus_poly: &Self) -> Result<Self, PolynomialError> {
        self.add(rhs)?.rem_mod_poly(modulus_poly)
    }

    /// Subtracts two polynomials and reduces modulo `modulus_poly`.
    pub fn sub_mod_poly(&self, rhs: &Self, modulus_poly: &Self) -> Result<Self, PolynomialError> {
        self.sub(rhs)?.rem_mod_poly(modulus_poly)
    }

    /// Multiplies two polynomials and reduces modulo `modulus_poly`.
    ///
    /// This computes the exact product first (in a temporary buffer sized to exact degree),
    /// then applies polynomial reduction.
    ///
    /// # Errors
    /// Returns compatibility and reduction-related errors as described by [`Self::rem_mod_poly`].
    pub fn mul_mod_poly(&self, rhs: &Self, modulus_poly: &Self) -> Result<Self, PolynomialError> {
        self.ensure_compatible(rhs)?;
        self.ensure_compatible(modulus_poly)?;

        let (Some(lhs_degree), Some(rhs_degree)) = (self.degree(), rhs.degree()) else {
            return Self::zero(self.max_degree, self.modulus);
        };

        let mut product = vec![0_u64; lhs_degree + rhs_degree + 1];
        for i in 0..=lhs_degree {
            let a = self.coeffs[i];
            if a == 0 {
                continue;
            }
            for j in 0..=rhs_degree {
                let b = rhs.coeffs[j];
                if b == 0 {
                    continue;
                }
                let idx = i + j;
                // Exact product accumulation before quotient-ring reduction.
                let term = mod_mul(a, b, self.modulus);
                product[idx] = mod_add(product[idx], term, self.modulus);
            }
        }

        reduce_coefficients_mod_poly(&mut product, &modulus_poly.coeffs, self.modulus)?;

        let mut out = vec![0_u64; self.max_degree + 1];
        for (i, coeff) in product.into_iter().enumerate().take(self.max_degree + 1) {
            out[i] = coeff;
        }

        Ok(Self {
            coeffs: out,
            modulus: self.modulus,
            max_degree: self.max_degree,
        })
    }

    /// Evaluates polynomial at `x` using Horner's method in `Z_q`.
    pub fn evaluate(&self, x: u64) -> u64 {
        let x_mod = x % self.modulus;
        let mut acc = 0_u64;

        for &coeff in self.coeffs.iter().rev() {
            acc = mod_mul(acc, x_mod, self.modulus);
            acc = mod_add(acc, coeff, self.modulus);
        }

        acc
    }

    /// Returns formal derivative `d/dx(self)` in `Z_q[x]`.
    pub fn derivative(&self) -> Self {
        let mut out = vec![0_u64; self.max_degree + 1];
        for (deg, &coeff) in self.coeffs.iter().enumerate().skip(1) {
            let factor = deg as u64 % self.modulus;
            out[deg - 1] = mod_mul(factor, coeff, self.modulus);
        }

        Self {
            coeffs: out,
            modulus: self.modulus,
            max_degree: self.max_degree,
        }
    }

    /// Polynomial long division in `Z_q[x]`: returns `(quotient, remainder)`.
    ///
    /// # Errors
    /// Returns:
    /// - [`PolynomialError::IncompatiblePolynomials`] when `n` or `q` differ,
    /// - [`PolynomialError::DivisionByZeroPolynomial`] if divisor is zero,
    /// - [`PolynomialError::NonInvertibleCoefficient`] if divisor leading coefficient has no
    ///   inverse in `Z_q`.
    pub fn div_rem(&self, divisor: &Self) -> Result<(Self, Self), PolynomialError> {
        self.ensure_compatible(divisor)?;
        if divisor.is_zero() {
            return Err(PolynomialError::DivisionByZeroPolynomial);
        }
        if self.is_zero() {
            return Ok((
                Self::zero(self.max_degree, self.modulus)?,
                Self::zero(self.max_degree, self.modulus)?,
            ));
        }

        let divisor_degree = divisor.degree().expect("checked non-zero divisor");
        let lead = divisor.coeffs[divisor_degree];
        let lead_inv =
            mod_inverse(lead, self.modulus).ok_or(PolynomialError::NonInvertibleCoefficient {
                coefficient: lead,
                modulus: self.modulus,
            })?;

        let mut remainder = self.clone();
        let mut quotient = Self::zero(self.max_degree, self.modulus)?;

        while let Some(rem_deg) = remainder.degree() {
            if rem_deg < divisor_degree {
                break;
            }

            // Cancel the current highest remainder term.
            let diff = rem_deg - divisor_degree;
            let rem_lead = remainder.coeffs[rem_deg];
            let factor = mod_mul(rem_lead, lead_inv, self.modulus);
            quotient.coeffs[diff] = mod_add(quotient.coeffs[diff], factor, self.modulus);

            // remainder -= factor * x^diff * divisor
            for i in 0..=divisor_degree {
                let idx = i + diff;
                let scaled = mod_mul(factor, divisor.coeffs[i], self.modulus);
                remainder.coeffs[idx] = mod_sub(remainder.coeffs[idx], scaled, self.modulus);
            }
        }

        Ok((quotient, remainder))
    }

    /// Ensures both polynomials share identical `(max_degree, modulus)` metadata.
    fn ensure_compatible(&self, rhs: &Self) -> Result<(), PolynomialError> {
        if self.modulus != rhs.modulus || self.max_degree != rhs.max_degree {
            return Err(PolynomialError::IncompatiblePolynomials);
        }
        Ok(())
    }
}

/// Computes `(a + b) mod modulus` using widened arithmetic to avoid overflow.
fn mod_add(a: u64, b: u64, modulus: u64) -> u64 {
    ((a as u128 + b as u128) % modulus as u128) as u64
}

/// Computes `(a - b) mod modulus` using widened arithmetic to avoid underflow.
fn mod_sub(a: u64, b: u64, modulus: u64) -> u64 {
    ((a as u128 + modulus as u128 - b as u128) % modulus as u128) as u64
}

/// Computes `(a * b) mod modulus` using widened arithmetic to avoid overflow.
fn mod_mul(a: u64, b: u64, modulus: u64) -> u64 {
    ((a as u128 * b as u128) % modulus as u128) as u64
}

/// Computes `base^exponent mod modulus` via binary exponentiation.
fn mod_pow(mut base: u64, mut exponent: u64, modulus: u64) -> u64 {
    let mut acc = 1_u64 % modulus;
    base %= modulus;

    while exponent > 0 {
        if exponent & 1 == 1 {
            acc = mod_mul(acc, base, modulus);
        }
        base = mod_mul(base, base, modulus);
        exponent >>= 1;
    }

    acc
}

/// Reorders values in-place using bit-reversed indices.
///
/// This is the standard preprocessing step for iterative Cooley-Tukey NTT.
fn bit_reverse_permute(values: &mut [u64]) {
    let n = values.len();
    let mut j = 0_usize;

    for i in 1..n {
        let mut bit = n >> 1;
        while j & bit != 0 {
            j ^= bit;
            bit >>= 1;
        }
        j ^= bit;

        if i < j {
            values.swap(i, j);
        }
    }
}

/// In-place iterative radix-2 NTT.
///
/// `root` must be an `n`-th root of unity in `Z_modulus`, where `n = values.len()`.
fn ntt_in_place(values: &mut [u64], root: u64, modulus: u64) -> Result<(), PolynomialError> {
    let n = values.len();
    if n == 0 || !n.is_power_of_two() {
        return Err(PolynomialError::NttLengthMustBePowerOfTwo(n));
    }

    // The butterfly stages assume bit-reversed layout.
    bit_reverse_permute(values);

    let mut len = 2_usize;
    while len <= n {
        // Twiddle factor increment for this stage.
        let wlen = mod_pow(root, (n / len) as u64, modulus);
        for start in (0..n).step_by(len) {
            let mut w = 1_u64;
            for i in 0..(len / 2) {
                // Cooley-Tukey butterfly:
                // (u, v) -> (u + w*v, u - w*v)
                let u = values[start + i];
                let v = mod_mul(values[start + i + len / 2], w, modulus);
                values[start + i] = mod_add(u, v, modulus);
                values[start + i + len / 2] = mod_sub(u, v, modulus);
                w = mod_mul(w, wlen, modulus);
            }
        }
        len <<= 1;
    }

    Ok(())
}

/// Convolution via NTT for coefficient vectors over `Z_modulus`.
///
/// `primitive_root` should be a primitive root modulo `modulus`.
fn convolution_ntt(
    lhs: &[u64],
    rhs: &[u64],
    modulus: u64,
    primitive_root: u64,
) -> Result<Vec<u64>, PolynomialError> {
    if lhs.is_empty() || rhs.is_empty() {
        return Ok(vec![0]);
    }

    // Convolution size and transform size.
    let out_len = lhs.len() + rhs.len() - 1;
    let ntt_len = out_len.next_power_of_two();
    if !ntt_len.is_power_of_two() {
        return Err(PolynomialError::NttLengthMustBePowerOfTwo(ntt_len));
    }

    let ntt_len_u64 = ntt_len as u64;
    if (modulus - 1) % ntt_len_u64 != 0 {
        return Err(PolynomialError::NttLengthUnsupported {
            length: ntt_len,
            modulus,
        });
    }

    // Derive primitive n-th root of unity from primitive root of the field.
    let root = mod_pow(primitive_root, (modulus - 1) / ntt_len_u64, modulus);
    let is_valid_root = mod_pow(root, ntt_len_u64, modulus) == 1
        && (ntt_len == 1 || mod_pow(root, (ntt_len / 2) as u64, modulus) != 1);
    if !is_valid_root {
        return Err(PolynomialError::InvalidPrimitiveRoot {
            primitive_root,
            length: ntt_len,
            modulus,
        });
    }

    // Inverse transform constants.
    let root_inv = mod_inverse(root, modulus).ok_or(PolynomialError::NonInvertibleCoefficient {
        coefficient: root,
        modulus,
    })?;
    let n_inv = mod_inverse(ntt_len_u64 % modulus, modulus).ok_or(
        PolynomialError::NonInvertibleCoefficient {
            coefficient: ntt_len_u64 % modulus,
            modulus,
        },
    )?;

    // Zero-pad inputs to the transform size.
    let mut fa = vec![0_u64; ntt_len];
    let mut fb = vec![0_u64; ntt_len];
    for (i, coeff) in lhs.iter().copied().enumerate() {
        fa[i] = coeff % modulus;
    }
    for (i, coeff) in rhs.iter().copied().enumerate() {
        fb[i] = coeff % modulus;
    }

    // Forward transforms.
    ntt_in_place(&mut fa, root, modulus)?;
    ntt_in_place(&mut fb, root, modulus)?;

    // Point-wise multiplication in evaluation domain.
    for (a, b) in fa.iter_mut().zip(fb.iter()) {
        *a = mod_mul(*a, *b, modulus);
    }

    // Inverse transform and normalization by 1/n.
    ntt_in_place(&mut fa, root_inv, modulus)?;
    for coeff in &mut fa {
        *coeff = mod_mul(*coeff, n_inv, modulus);
    }
    // Trim back to exact convolution length.
    fa.truncate(out_len);

    Ok(fa)
}

/// Computes multiplicative inverse of `a mod modulus` using extended Euclid.
///
/// Returns `None` if `gcd(a, modulus) != 1`.
fn mod_inverse(a: u64, modulus: u64) -> Option<u64> {
    let mut t = 0_i128;
    let mut new_t = 1_i128;
    let mut r = modulus as i128;
    let mut new_r = (a % modulus) as i128;

    while new_r != 0 {
        let quotient = r / new_r;
        (t, new_t) = (new_t, t - quotient * new_t);
        (r, new_r) = (new_r, r - quotient * new_r);
    }

    if r != 1 {
        return None;
    }

    if t < 0 {
        t += modulus as i128;
    }
    Some(t as u64)
}

/// Returns the highest index whose coefficient is non-zero.
fn degree_of(coeffs: &[u64]) -> Option<usize> {
    coeffs.iter().rposition(|&c| c != 0)
}

/// Reduces a mutable coefficient buffer modulo `modulus_poly` in-place.
///
/// This performs polynomial long-division style elimination on the highest degree terms.
fn reduce_coefficients_mod_poly(
    coeffs: &mut [u64],
    modulus_poly: &[u64],
    modulus: u64,
) -> Result<(), PolynomialError> {
    let Some(modulus_degree) = degree_of(modulus_poly) else {
        return Err(PolynomialError::DivisionByZeroPolynomial);
    };

    let leading = modulus_poly[modulus_degree];
    let leading_inverse =
        mod_inverse(leading, modulus).ok_or(PolynomialError::NonInvertibleCoefficient {
            coefficient: leading,
            modulus,
        })?;

    while let Some(current_degree) = degree_of(coeffs) {
        if current_degree < modulus_degree {
            break;
        }

        // Eliminate the current leading term using shifted/scaled modulus polynomial.
        let shift = current_degree - modulus_degree;
        let factor = mod_mul(coeffs[current_degree], leading_inverse, modulus);
        if factor == 0 {
            continue;
        }

        for (i, &modulus_coeff) in modulus_poly.iter().take(modulus_degree + 1).enumerate() {
            let idx = shift + i;
            if idx >= coeffs.len() || modulus_coeff == 0 {
                continue;
            }
            let scaled = mod_mul(factor, modulus_coeff, modulus);
            coeffs[idx] = mod_sub(coeffs[idx], scaled, modulus);
        }
    }

    Ok(())
}

#[cfg(test)]
mod tests {
    use super::{Polynomial, PolynomialError, mod_inverse, mod_mul, mod_pow, ntt_in_place};

    fn p(n: usize, q: u64, coeffs: &[u64]) -> Polynomial {
        Polynomial::new(n, q, coeffs).expect("polynomial should build")
    }

    fn lcg_next(state: &mut u64) -> u64 {
        *state = state
            .wrapping_mul(6_364_136_223_846_793_005)
            .wrapping_add(1_442_695_040_888_963_407);
        *state
    }

    fn enumerate_polynomials(
        max_degree: usize,
        modulus: u64,
        values: &[u64],
        used_len: usize,
    ) -> Vec<Polynomial> {
        assert!(
            used_len <= max_degree + 1,
            "used_len exceeds polynomial capacity"
        );
        assert!(!values.is_empty(), "values cannot be empty");

        let total = values.len().pow(used_len as u32);
        let mut out = Vec::with_capacity(total);
        for mut state in 0..total {
            let mut coeffs = vec![0_u64; used_len];
            for slot in &mut coeffs {
                let digit = state % values.len();
                *slot = values[digit];
                state /= values.len();
            }
            out.push(p(max_degree, modulus, &coeffs));
        }
        out
    }

    #[test]
    fn constructor_normalizes_and_pads_coefficients() {
        let poly = p(5, 7, &[10, 15, 6]);
        assert_eq!(poly.coefficients(), &[3, 1, 6, 0, 0, 0]);
    }

    #[test]
    fn constructor_rejects_invalid_modulus() {
        let err = Polynomial::new(3, 1, &[1, 2]).expect_err("expected error");
        assert_eq!(err, PolynomialError::InvalidModulus(1));
    }

    #[test]
    fn constructor_rejects_degree_overflow() {
        let err = Polynomial::new(2, 11, &[1, 2, 3, 4]).expect_err("expected error");
        assert_eq!(
            err,
            PolynomialError::DegreeOverflow {
                requested: 3,
                max_degree: 2
            }
        );
    }

    #[test]
    fn degree_and_zero_behaviour() {
        let zero = p(4, 13, &[]);
        assert!(zero.is_zero());
        assert_eq!(zero.degree(), None);
        assert_eq!(zero.trimmed_coefficients(), vec![0]);

        let poly = p(4, 13, &[0, 2, 0, 9]);
        assert!(!poly.is_zero());
        assert_eq!(poly.degree(), Some(3));
        assert_eq!(poly.trimmed_coefficients(), vec![0, 2, 0, 9]);
    }

    #[test]
    fn addition_and_subtraction_work_modulo_q() {
        let a = p(5, 17, &[16, 5, 0, 1]);
        let b = p(5, 17, &[3, 13, 6, 10]);

        let sum = a.add(&b).expect("add should work");
        assert_eq!(sum.trimmed_coefficients(), vec![2, 1, 6, 11]);

        let diff = a.sub(&b).expect("sub should work");
        assert_eq!(diff.trimmed_coefficients(), vec![13, 9, 11, 8]);
    }

    #[test]
    fn unary_negation_and_scalar_multiplication() {
        let a = p(4, 19, &[0, 4, 18, 3]);
        let neg = a.neg();
        assert_eq!(neg.trimmed_coefficients(), vec![0, 15, 1, 16]);

        let scaled = a.scalar_mul(7);
        assert_eq!(scaled.trimmed_coefficients(), vec![0, 9, 12, 2]);
    }

    #[test]
    fn incompatible_polynomials_fail_for_binary_ops() {
        let a = p(3, 11, &[1, 2]);
        let b = p(4, 11, &[1, 2]);
        let err = a.add(&b).expect_err("expected incompatibility");
        assert_eq!(err, PolynomialError::IncompatiblePolynomials);
    }

    #[test]
    fn multiplication_checked_and_truncated() {
        let a = p(5, 23, &[1, 2, 3]);
        let b = p(5, 23, &[4, 5]);

        let prod = a.mul(&b).expect("mul should work");
        assert_eq!(prod.trimmed_coefficients(), vec![4, 13, 22, 15]);

        let c = p(3, 29, &[1, 2, 3, 4]);
        let d = p(3, 29, &[1, 1, 1, 1]);
        let err = c.mul(&d).expect_err("degree should overflow");
        assert_eq!(
            err,
            PolynomialError::ProductDegreeOverflow {
                degree: 6,
                max_degree: 3
            }
        );

        let truncated = c.mul_truncated(&d).expect("mul_truncated should work");
        assert_eq!(truncated.coefficients(), &[1, 3, 6, 10]);
    }

    #[test]
    fn schoolbook_convolution_matches_expected_coefficients() {
        let a = p(10, 998_244_353, &[1, 2, 3, 4]);
        let b = p(10, 998_244_353, &[5, 6, 7]);
        let product = a.mul(&b).expect("schoolbook multiplication should work");

        // [1,2,3,4] * [5,6,7] = [5,16,34,52,45,28]
        assert_eq!(product.trimmed_coefficients(), vec![5, 16, 34, 52, 45, 28]);
    }

    #[test]
    fn ntt_round_trip_recovers_original_vector() {
        let modulus = 998_244_353_u64;
        let primitive_root = 3_u64;
        let n = 8_usize;

        let omega = mod_pow(primitive_root, (modulus - 1) / n as u64, modulus);
        let omega_inv = mod_inverse(omega, modulus).expect("omega must be invertible");
        let n_inv = mod_inverse(n as u64, modulus).expect("length must be invertible");

        let mut values = vec![7, 11, 19, 23, 31, 2, 5, 13];
        let original = values.clone();

        ntt_in_place(&mut values, omega, modulus).expect("forward NTT should work");
        ntt_in_place(&mut values, omega_inv, modulus).expect("inverse NTT should work");
        for value in &mut values {
            *value = mod_mul(*value, n_inv, modulus);
        }

        assert_eq!(values, original);
    }

    #[test]
    fn ntt_multiplication_matches_schoolbook_convolution() {
        let modulus = 998_244_353_u64;
        let primitive_root = 3_u64;

        let a = p(31, modulus, &[4, 1, 9, 16, 25, 36, 49, 64, 81, 100]);
        let b = p(31, modulus, &[3, 14, 15, 92, 65, 35, 89, 79]);

        let schoolbook = a.mul(&b).expect("schoolbook multiplication should work");
        let ntt = a
            .mul_ntt(&b, primitive_root)
            .expect("NTT multiplication should work");

        assert_eq!(ntt.coefficients(), schoolbook.coefficients());
    }

    #[test]
    fn ntt_matches_schoolbook_on_many_deterministic_cases() {
        let modulus = 998_244_353_u64;
        let primitive_root = 3_u64;
        let mut seed = 0xC0FFEE_u64;

        for _ in 0..200 {
            let len_a = (lcg_next(&mut seed) % 48 + 1) as usize;
            let len_b = (lcg_next(&mut seed) % 48 + 1) as usize;
            let max_degree = len_a + len_b;

            let mut coeffs_a = vec![0_u64; len_a];
            let mut coeffs_b = vec![0_u64; len_b];
            for coeff in &mut coeffs_a {
                *coeff = lcg_next(&mut seed) % modulus;
            }
            for coeff in &mut coeffs_b {
                *coeff = lcg_next(&mut seed) % modulus;
            }

            let a = p(max_degree, modulus, &coeffs_a);
            let b = p(max_degree, modulus, &coeffs_b);
            let schoolbook = a.mul(&b).expect("schoolbook multiplication should work");
            let ntt = a
                .mul_ntt(&b, primitive_root)
                .expect("NTT multiplication should work");

            assert_eq!(ntt.coefficients(), schoolbook.coefficients());
        }
    }

    #[test]
    fn ntt_errors_are_reported_for_bad_parameters() {
        // out_len = 21, next power of two is 32, but 32 does not divide q-1 = 16.
        let coeffs = vec![1_u64; 11];
        let a = p(20, 17, &coeffs);
        let b = p(20, 17, &coeffs);
        let err = a
            .mul_ntt(&b, 3)
            .expect_err("unsupported NTT length should fail");
        assert_eq!(
            err,
            PolynomialError::NttLengthUnsupported {
                length: 32,
                modulus: 17
            }
        );

        let c = p(16, 998_244_353, &[1, 2, 3, 4]);
        let d = p(16, 998_244_353, &[5, 6, 7, 8]);
        let err = c
            .mul_ntt(&d, 1)
            .expect_err("invalid primitive root should fail");
        assert_eq!(
            err,
            PolynomialError::InvalidPrimitiveRoot {
                primitive_root: 1,
                length: 8,
                modulus: 998_244_353
            }
        );
    }

    #[test]
    fn remainder_mod_polynomial_reduces_degree() {
        // Modulus m(x) = x^4 + 1 over mod 17.
        let m = p(6, 17, &[1, 0, 0, 0, 1]);
        let a = p(6, 17, &[2, 0, 0, 0, 1, 3]); // 2 + x^4 + 3x^5
        let reduced = a.rem_mod_poly(&m).expect("reduction should work");

        // x^4 = -1 and x^5 = -x, so result is 1 - 3x = 1 + 14x mod 17.
        assert_eq!(reduced.trimmed_coefficients(), vec![1, 14]);
    }

    #[test]
    fn quotient_ring_add_sub_and_mul() {
        // Work in Z_17[x] / (x^4 + 1).
        let m = p(4, 17, &[1, 0, 0, 0, 1]);
        let a = p(4, 17, &[0, 0, 0, 16]); // -x^3
        let b = p(4, 17, &[0, 0, 0, 2]); // 2x^3

        let added = a.add_mod_poly(&b, &m).expect("add mod poly should work");
        assert_eq!(added.trimmed_coefficients(), vec![0, 0, 0, 1]);

        let subbed = a.sub_mod_poly(&b, &m).expect("sub mod poly should work");
        assert_eq!(subbed.trimmed_coefficients(), vec![0, 0, 0, 14]);

        let c = p(4, 17, &[1, 0, 0, 1]); // 1 + x^3
        let d = p(4, 17, &[1, 0, 0, 1]); // 1 + x^3
        let err = c.mul(&d).expect_err("plain mul overflows max_degree");
        assert_eq!(
            err,
            PolynomialError::ProductDegreeOverflow {
                degree: 6,
                max_degree: 4
            }
        );

        let ring_product = c.mul_mod_poly(&d, &m).expect("mul mod poly should work");
        // (1 + x^3)^2 = 1 + 2x^3 + x^6, x^6 = -x^2 over x^4 + 1.
        assert_eq!(ring_product.trimmed_coefficients(), vec![1, 0, 16, 2]);
    }

    #[test]
    fn polynomial_modulus_errors_are_reported() {
        let a = p(4, 11, &[1, 2, 3]);
        let b = p(4, 11, &[3, 4]);
        let zero_poly = p(4, 11, &[]);
        let err = a
            .mul_mod_poly(&b, &zero_poly)
            .expect_err("zero modulus polynomial should fail");
        assert_eq!(err, PolynomialError::DivisionByZeroPolynomial);

        // Leading coefficient 2 has no inverse mod 8.
        let x = p(4, 8, &[0, 1]);
        let non_invertible_modulus = p(4, 8, &[1, 2]); // 1 + 2x
        let err = x
            .rem_mod_poly(&non_invertible_modulus)
            .expect_err("non-invertible modulus lead should fail");
        assert_eq!(
            err,
            PolynomialError::NonInvertibleCoefficient {
                coefficient: 2,
                modulus: 8
            }
        );

        let c = p(5, 11, &[1, 2]);
        let err = a
            .add_mod_poly(&b, &c)
            .expect_err("incompatible polynomial settings should fail");
        assert_eq!(err, PolynomialError::IncompatiblePolynomials);
    }

    #[test]
    fn evaluate_uses_horner_rule_modulo_q() {
        let poly = p(4, 31, &[7, 0, 3, 4]); // 7 + 3x^2 + 4x^3
        let value = poly.evaluate(10);
        // 7 + 3*100 + 4*1000 = 4307, 4307 mod 31 = 29
        assert_eq!(value, 29);
    }

    #[test]
    fn derivative_is_computed_modulo_q() {
        let poly = p(6, 11, &[3, 5, 7, 9, 2]); // 3 + 5x + 7x^2 + 9x^3 + 2x^4
        let deriv = poly.derivative(); // 5 + 14x + 27x^2 + 8x^3 mod 11
        assert_eq!(deriv.trimmed_coefficients(), vec![5, 3, 5, 8]);
    }

    #[test]
    fn long_division_exact_case() {
        // (x^3 + 2x^2 + 6x + 5) / (x + 1) over mod 7 = x^2 + x + 5, remainder 0
        let dividend = p(5, 7, &[5, 6, 2, 1]);
        let divisor = p(5, 7, &[1, 1]);
        let (quotient, remainder) = dividend.div_rem(&divisor).expect("division should work");

        assert_eq!(quotient.trimmed_coefficients(), vec![5, 1, 1]);
        assert!(remainder.is_zero());
    }

    #[test]
    fn long_division_with_remainder() {
        // (3 + 0x + x^2) / (2 + x) over mod 5
        let dividend = p(4, 5, &[3, 0, 1]);
        let divisor = p(4, 5, &[2, 1]);
        let (quotient, remainder) = dividend.div_rem(&divisor).expect("division should work");

        assert_eq!(quotient.trimmed_coefficients(), vec![3, 1]);
        assert_eq!(remainder.trimmed_coefficients(), vec![2]);

        // Validate dividend = divisor * quotient + remainder
        let reconstructed = divisor
            .mul(&quotient)
            .expect("reconstruction product should fit")
            .add(&remainder)
            .expect("reconstruction sum should fit");
        assert_eq!(reconstructed.coefficients(), dividend.coefficients());
    }

    #[test]
    fn division_errors_are_reported() {
        let a = p(4, 9, &[1, 2, 3]);
        let zero = p(4, 9, &[]);
        let err = a.div_rem(&zero).expect_err("division by zero should fail");
        assert_eq!(err, PolynomialError::DivisionByZeroPolynomial);

        // Leading coefficient is 2, not invertible mod 8.
        let dividend = p(4, 8, &[1, 0, 1]);
        let divisor = p(4, 8, &[0, 2]);
        let err = dividend
            .div_rem(&divisor)
            .expect_err("division should fail when inverse does not exist");
        assert_eq!(
            err,
            PolynomialError::NonInvertibleCoefficient {
                coefficient: 2,
                modulus: 8
            }
        );
    }

    #[test]
    fn mod_q_algebraic_laws_hold_on_small_exhaustive_domain() {
        let modulus = 5_u64;
        let max_degree = 2_usize;
        let all = enumerate_polynomials(max_degree, modulus, &[0, 1, 2, 3, 4], 3);
        let reps = enumerate_polynomials(max_degree, modulus, &[0, 1, 4], 3);
        let zero = p(max_degree, modulus, &[]);
        let one = p(max_degree, modulus, &[1]);

        for a in &all {
            assert_eq!(a.add(&zero).expect("a + 0 should work"), *a);
            assert_eq!(zero.add(a).expect("0 + a should work"), *a);
            assert_eq!(a.sub(&zero).expect("a - 0 should work"), *a);
            assert!(
                a.add(&a.neg()).expect("a + (-a) should work").is_zero(),
                "additive inverse should cancel"
            );
            assert_eq!(a.scalar_mul(modulus + 2), a.scalar_mul(2));
            assert_eq!(a.scalar_mul(0), zero);
            assert_eq!(
                a.mul_truncated(&one).expect("a * 1 should work"),
                *a,
                "multiplicative identity should hold"
            );
        }

        for a in &all {
            for b in &all {
                assert_eq!(
                    a.add(b).expect("a+b should work"),
                    b.add(a).expect("b+a should work"),
                    "addition should commute"
                );
                assert_eq!(
                    a.mul_truncated(b).expect("a*b should work"),
                    b.mul_truncated(a).expect("b*a should work"),
                    "truncated multiplication should commute"
                );
                assert_eq!(
                    a.sub(b).expect("a-b should work"),
                    a.add(&b.neg()).expect("a+(-b) should work"),
                    "subtraction should match addition with negation"
                );
            }
        }

        for a in &reps {
            for b in &reps {
                for c in &reps {
                    assert_eq!(
                        a.add(&b.add(c).expect("b+c should work"))
                            .expect("a+(b+c) should work"),
                        a.add(b)
                            .expect("a+b should work")
                            .add(c)
                            .expect("(a+b)+c should work"),
                        "addition should associate"
                    );

                    let lhs = a
                        .mul_truncated(&b.add(c).expect("b+c should work"))
                        .expect("a*(b+c) should work");
                    let rhs = a
                        .mul_truncated(b)
                        .expect("a*b should work")
                        .add(&a.mul_truncated(c).expect("a*c should work"))
                        .expect("ab+ac should work");
                    assert_eq!(lhs, rhs, "multiplication should distribute over addition");
                }
            }
        }
    }

    #[test]
    fn evaluate_matches_naive_formula_on_exhaustive_small_domain() {
        let modulus = 11_u64;
        let max_degree = 3_usize;
        let polys = enumerate_polynomials(max_degree, modulus, &[0, 1, 2], 4);

        for poly in &polys {
            for x in 0..modulus {
                let mut acc = 0_u64;
                let mut x_pow = 1_u64;
                for &coeff in poly.coefficients() {
                    acc = (acc + (coeff * x_pow) % modulus) % modulus;
                    x_pow = (x_pow * x) % modulus;
                }
                assert_eq!(poly.evaluate(x), acc);
            }
        }
    }

    #[test]
    fn ntt_matches_schoolbook_exhaustive_small_domains() {
        // Degree <= 1 under q=5 gives out_len <= 3 -> ntt_len=4, supported since 4 | (5-1).
        let small_q_polys = enumerate_polynomials(2, 5, &[0, 1, 2, 3, 4], 2);
        for a in &small_q_polys {
            for b in &small_q_polys {
                let schoolbook = a.mul(b).expect("schoolbook should work");
                let ntt = a.mul_ntt(b, 2).expect("NTT should work");
                assert_eq!(ntt.coefficients(), schoolbook.coefficients());
            }
        }

        // Degree <= 2 under q=17 gives out_len <= 5 -> ntt_len=8, supported since 8 | (17-1).
        let larger_polys = enumerate_polynomials(4, 17, &[0, 1, 2, 3], 3);
        for a in &larger_polys {
            for b in &larger_polys {
                let schoolbook = a.mul(b).expect("schoolbook should work");
                let ntt = a.mul_ntt(b, 3).expect("NTT should work");
                assert_eq!(ntt.coefficients(), schoolbook.coefficients());
            }
        }
    }

    #[test]
    fn mod_polynomial_reduction_is_bounded_and_idempotent() {
        let modulus = 17_u64;
        let max_degree = 4_usize;
        let modulus_poly = p(max_degree, modulus, &[1, 0, 0, 0, 1]); // x^4 + 1
        let all = enumerate_polynomials(max_degree, modulus, &[0, 1, 2], 5);
        let sample: Vec<_> = all.iter().take(64).cloned().collect();

        for poly in &all {
            let reduced = poly
                .rem_mod_poly(&modulus_poly)
                .expect("reduction should succeed");
            if let Some(degree) = reduced.degree() {
                assert!(
                    degree < 4,
                    "canonical representative degree must be < deg(modulus)"
                );
            }
            assert_eq!(
                reduced,
                reduced
                    .rem_mod_poly(&modulus_poly)
                    .expect("reduction idempotence should hold")
            );
        }

        for a in &sample {
            for b in &sample {
                let lhs_add = a
                    .add_mod_poly(b, &modulus_poly)
                    .expect("add mod poly should work");
                let rhs_add = a
                    .rem_mod_poly(&modulus_poly)
                    .expect("reduction should work")
                    .add_mod_poly(
                        &b.rem_mod_poly(&modulus_poly)
                            .expect("reduction should work"),
                        &modulus_poly,
                    )
                    .expect("add mod poly should work");
                assert_eq!(lhs_add, rhs_add);

                let lhs_mul = a
                    .mul_mod_poly(b, &modulus_poly)
                    .expect("mul mod poly should work");
                let rhs_mul = a
                    .rem_mod_poly(&modulus_poly)
                    .expect("reduction should work")
                    .mul_mod_poly(
                        &b.rem_mod_poly(&modulus_poly)
                            .expect("reduction should work"),
                        &modulus_poly,
                    )
                    .expect("mul mod poly should work");
                assert_eq!(lhs_mul, rhs_mul);
            }
        }
    }

    #[test]
    fn long_division_invariants_hold_on_exhaustive_small_domain() {
        let modulus = 5_u64;
        let max_degree = 4_usize;
        let dividends = enumerate_polynomials(max_degree, modulus, &[0, 1, 2], 5);
        let divisors: Vec<_> = enumerate_polynomials(max_degree, modulus, &[0, 1, 2], 3)
            .into_iter()
            .filter(|poly| !poly.is_zero())
            .collect();

        for dividend in &dividends {
            for divisor in &divisors {
                let (quotient, remainder) =
                    dividend.div_rem(divisor).expect("division should succeed");

                let reconstructed = divisor
                    .mul(&quotient)
                    .expect("product should fit in max degree")
                    .add(&remainder)
                    .expect("sum should fit in max degree");
                assert_eq!(reconstructed.coefficients(), dividend.coefficients());

                if let Some(rem_degree) = remainder.degree() {
                    let divisor_degree = divisor.degree().expect("divisors are non-zero");
                    assert!(
                        rem_degree < divisor_degree,
                        "remainder degree must be strictly less than divisor degree"
                    );
                }
            }
        }
    }
}