malachite-base 0.13.0

// Copyright © 2026 Mikhail Hogrefe
//
// This file is part of Malachite.
//
// Malachite is free software: you can redistribute it and/or modify it under the terms of the GNU
// Lesser General Public License (LGPL) as published by the Free Software Foundation; either version
// 3 of the License, or (at your option) any later version. See <https://www.gnu.org/licenses/>.

use crate::num::arithmetic::mod_mul::{
    mod_mul_precompute_shoup, mod_mul_shoup, mod_mul_shoup_lazy,
};
use crate::num::arithmetic::traits::{ModIsReduced, ModPowerOf2IsReduced};
use crate::num::basic::integers::PrimitiveInt;
use crate::num::basic::unsigneds::PrimitiveUnsigned;
use crate::num::conversion::traits::ExactFrom;
use crate::polynomial::{ModEvaluate, ModEvaluateGeometric, ModEvaluateMany, ModPowerOf2Evaluate};
use crate::unsigned_polynomial::UnsignedPolynomial;
use crate::unsigned_polynomial::arithmetic::mod_mul::ModData;
use crate::unsigned_polynomial::arithmetic::mod_mul_middle::mod_mul_middle_karatsuba;
use crate::unsigned_polynomial::arithmetic::mod_mul_truncated::mod_mul_truncated_to_out;
use alloc::vec;
use alloc::vec::Vec;
use core::cmp::max;

// Evaluates a polynomial at x modulo 2^pow, after checking that pow fits in `T` and that the
// coefficients and x are reduced.
fn mod_power_of_2_evaluate<T: PrimitiveUnsigned>(p: &UnsignedPolynomial<T>, x: T, pow: u64) -> T {
    assert!(pow <= T::WIDTH);
    assert!(
        p.mod_power_of_2_is_reduced(pow),
        "self must be reduced mod 2^pow, but {p} has a coefficient >= 2^{pow}"
    );
    assert!(
        x.significant_bits() <= pow,
        "x must be reduced mod 2^pow, but {x} >= 2^{pow}"
    );
    let mut value = T::ZERO;
    for &c in p.coefficients.iter().rev() {
        value = value.wrapping_mul(x).wrapping_add(c);
    }
    value.mod_power_of_2(pow)
}

impl<T: PrimitiveUnsigned> ModPowerOf2Evaluate<T> for &UnsignedPolynomial<T> {
    type Output = T;

    /// Evaluates an [`UnsignedPolynomial`] at a value of its coefficient type, modulo $2^k$, taking
    /// the polynomial by reference. The coefficients and the value must already be reduced modulo
    /// $2^k$, and $k$ may be at most the width of the type.
    ///
    /// $$
    /// f(p, x, k) = \sum_{i=0}^{n-1} c_i x^i \bmod 2^k,
    /// $$
    ///
    /// where $c_i$ is the coefficient of $x^i$ in $p$ and $n$ is its length. The zero polynomial
    /// evaluates to 0 everywhere.
    ///
    /// Reducing modulo $2^k$ commutes with wrapping arithmetic, which works modulo $2^w$ for the
    /// type's width $w \geq k$, so Horner's rule is carried out with wrapping multiplications and
    /// additions and the value is reduced once, at the end.
    ///
    /// # Worst-case complexity
    /// $T(n) = O(n)$
    ///
    /// $M(n) = O(1)$
    ///
    /// where $T$ is time, $M$ is additional memory, and $n$ is `self.len()`.
    ///
    /// # Panics
    /// Panics if `pow` is greater than `T::WIDTH`, or if any coefficient of `self` or `x` is
    /// greater than or equal to $2^k$.
    ///
    /// # Examples
    /// ```
    /// use core::str::FromStr;
    /// use malachite_base::polynomial::ModPowerOf2Evaluate;
    /// use malachite_base::unsigned_polynomial::UnsignedPolynomial;
    ///
    /// let p = UnsignedPolynomial::<u8>::from_str("5*x^2+3*x+7").unwrap();
    /// // 5 * 36 + 3 * 6 + 7 = 205, which is 13 mod 16.
    /// assert_eq!((&p).mod_power_of_2_evaluate(6, 4), 13);
    /// assert_eq!((&p).mod_power_of_2_evaluate(0, 4), 7);
    /// // All 8 bits of a u8: 205 itself.
    /// assert_eq!((&p).mod_power_of_2_evaluate(6, 8), 205);
    /// ```
    ///
    /// This is equivalent to `nmod_poly_evaluate_nmod` from `nmod_poly/evaluate_nmod.c`, FLINT
    /// 3.6.0, with the modulus $2^k$, except that the value must be reduced.
    #[inline]
    fn mod_power_of_2_evaluate(self, x: T, pow: u64) -> T {
        mod_power_of_2_evaluate(self, x, pow)
    }
}

impl<T: PrimitiveUnsigned> ModPowerOf2Evaluate<T> for UnsignedPolynomial<T> {
    type Output = T;

    /// Evaluates an [`UnsignedPolynomial`] at a value of its coefficient type, modulo $2^k$, taking
    /// the polynomial by value. The coefficients and the value must already be reduced modulo
    /// $2^k$, and $k$ may be at most the width of the type.
    ///
    /// $$
    /// f(p, x, k) = \sum_{i=0}^{n-1} c_i x^i \bmod 2^k,
    /// $$
    ///
    /// where $c_i$ is the coefficient of $x^i$ in $p$ and $n$ is its length. The zero polynomial
    /// evaluates to 0 everywhere.
    ///
    /// Reducing modulo $2^k$ commutes with wrapping arithmetic, which works modulo $2^w$ for the
    /// type's width $w \geq k$, so Horner's rule is carried out with wrapping multiplications and
    /// additions and the value is reduced once, at the end.
    ///
    /// # Worst-case complexity
    /// $T(n) = O(n)$
    ///
    /// $M(n) = O(1)$
    ///
    /// where $T$ is time, $M$ is additional memory, and $n$ is `self.len()`.
    ///
    /// # Panics
    /// Panics if `pow` is greater than `T::WIDTH`, or if any coefficient of `self` or `x` is
    /// greater than or equal to $2^k$.
    ///
    /// # Examples
    /// ```
    /// use core::str::FromStr;
    /// use malachite_base::polynomial::ModPowerOf2Evaluate;
    /// use malachite_base::unsigned_polynomial::UnsignedPolynomial;
    ///
    /// let p = UnsignedPolynomial::<u8>::from_str("5*x^2+3*x+7").unwrap();
    /// // 5 * 36 + 3 * 6 + 7 = 205, which is 13 mod 16.
    /// assert_eq!(p.clone().mod_power_of_2_evaluate(6, 4), 13);
    /// assert_eq!(p.clone().mod_power_of_2_evaluate(0, 4), 7);
    /// // All 8 bits of a u8: 205 itself.
    /// assert_eq!(p.mod_power_of_2_evaluate(6, 8), 205);
    /// ```
    ///
    /// This is equivalent to `nmod_poly_evaluate_nmod` from `nmod_poly/evaluate_nmod.c`, FLINT
    /// 3.6.0, with the modulus $2^k$, except that the value must be reduced.
    #[inline]
    fn mod_power_of_2_evaluate(self, x: T, pow: u64) -> T {
        mod_power_of_2_evaluate(&self, x, pow)
    }
}

// The shortest polynomial evaluated with Shoup's method, when `T` is wider than 32 bits: Shoup's
// precomputation is a two-by-one division, which is only paid back from this length on. For
// narrower types every polynomial of length 2 or more uses it. Tuned on Apple M-series for `u64`
// and `u128`; FLINT's `FLINT_MULMOD_SHOUP_THRESHOLD` is 10.
const MOD_EVALUATE_SHOUP_THRESHOLD: usize = 3;

// Evaluates a polynomial at `x` modulo `m` with Horner's rule, reducing after every step.
// `coefficients` must be nonempty.
//
// This is equivalent to `_nmod_poly_evaluate_nmod_horner` from `nmod_poly/evaluate_nmod.c`, FLINT
// 3.6.0.
crate_test_fn! {mod_evaluate_horner<T: PrimitiveUnsigned>(coefficients: &[T], x: T, m: T) -> T {
    let data = T::precompute_mod_mul_data(&m);
    let (&last, rest) = coefficients.split_last().unwrap();
    let mut value = last;
    for &c in rest.iter().rev() {
        value.mod_mul_precomputed_assign(x, m, &data);
        value.mod_add_assign(c, m);
    }
    value
}}

// Evaluates a polynomial at `x` modulo `m` with Horner's rule, multiplying by `x` with Shoup's
// method. `coefficients` must be nonempty, `x_precomp` must be `mod_mul_precompute_shoup(x, m)`,
// and the top bit of `m` must be clear.
//
// This is equivalent to `_nmod_poly_evaluate_nmod_precomp` from `nmod_poly/evaluate_nmod.c`, FLINT
// 3.6.0.
crate_test_fn! {mod_evaluate_shoup<T: PrimitiveUnsigned>(
    coefficients: &[T],
    x: T,
    x_precomp: T,
    m: T,
) -> T {
    let (&last, rest) = coefficients.split_last().unwrap();
    let mut value = last;
    for &c in rest.iter().rev() {
        value = mod_mul_shoup(x, value, x_precomp, m);
        value.mod_add_assign(c, m);
    }
    value
}}

// Evaluates a polynomial at `x` modulo `m` like `mod_evaluate_shoup`, but reduces only partially:
// the result is congruent to the polynomial's value and less than $3m - 1$. `coefficients` must be
// nonempty, `x_precomp` must be `mod_mul_precompute_shoup(x, m)`, and `m` must be at most `T::MAX /
// 3`, so that $3m - 1$ values fit.
//
// This is equivalent to `_nmod_poly_evaluate_nmod_precomp_lazy` from `nmod_poly/evaluate_nmod.c`,
// FLINT 3.6.0.
crate_test_fn! {mod_evaluate_shoup_lazy<T: PrimitiveUnsigned>(
    coefficients: &[T],
    x: T,
    x_precomp: T,
    m: T,
) -> T {
    let (&last, rest) = coefficients.split_last().unwrap();
    let mut value = last;
    for &c in rest.iter().rev() {
        // value is x * value mod m, or that plus m
        value = mod_mul_shoup_lazy(x, value, x_precomp, m);
        // value is now less than 3m - 1, since c < m
        value.wrapping_add_assign(c);
    }
    value
}}

// Evaluates the polynomial with the given coefficients at `x` modulo `m`. The coefficients and `x`
// must be less than `m`; this is not checked. It is public, but hidden, because malachite-nz
// evaluates an `IntegerPolynomial` modulo a word with it, after reducing the coefficients.
//
// This is equivalent to `_nmod_poly_evaluate_nmod` from `nmod_poly/evaluate_nmod.c`, FLINT 3.6.0,
// except that rectangular splitting is not used.
#[doc(hidden)]
pub fn mod_evaluate_slice<T: PrimitiveUnsigned>(coefficients: &[T], x: T, m: T) -> T {
    let len = coefficients.len();
    if len == 0 {
        return T::ZERO;
    }
    if len == 1 || x == T::ZERO {
        return coefficients[0];
    }
    // Shoup's method needs the top bit of m clear
    if m.get_highest_bit() || (T::WIDTH > u32::WIDTH && len < MOD_EVALUATE_SHOUP_THRESHOLD) {
        return mod_evaluate_horner(coefficients, x, m);
    }
    let x_precomp = mod_mul_precompute_shoup(x, m);
    // The lazy loop's values are less than 3m - 1, so it is used when those fit: when m <= (2^W +
    // 1) / 3, which is T::MAX / 3 since every width W is even. FLINT calls this bound LAZY_MAX.
    if m <= T::MAX / T::from(3u8) {
        let mut value = mod_evaluate_shoup_lazy(coefficients, x, x_precomp, m);
        // correct the excess
        let two_m = m << 1;
        if value >= two_m {
            value -= two_m;
        } else if value >= m {
            value -= m;
        }
        value
    } else {
        mod_evaluate_shoup(coefficients, x, x_precomp, m)
    }
}

fn mod_evaluate<T: PrimitiveUnsigned>(p: &UnsignedPolynomial<T>, x: T, m: T) -> T {
    assert!(
        p.mod_is_reduced(&m),
        "self must be reduced mod m, but {p} has a coefficient >= {m}"
    );
    assert!(x < m, "x must be reduced mod m, but {x} >= {m}");
    mod_evaluate_slice(&p.coefficients, x, m)
}

impl<T: PrimitiveUnsigned> ModEvaluate<T> for &UnsignedPolynomial<T> {
    type Output = T;

    /// Evaluates an [`UnsignedPolynomial`] at a value of its coefficient type, modulo a value of
    /// that type, taking the polynomial by reference. The coefficients and the value must already
    /// be reduced modulo `m`.
    ///
    /// $$
    /// f(p, x, m) = \sum_{i=0}^{n-1} c_i x^i \bmod m,
    /// $$
    ///
    /// where $c_i$ is the coefficient of $x^i$ in $p$ and $n$ is its length. The zero polynomial
    /// evaluates to 0 everywhere.
    ///
    /// The value is found with Horner's rule, reducing after every step. When the polynomial is
    /// long enough and the top bit of `m` is clear, every multiplication is by the same `x`, so
    /// Shoup's method is used: $\lfloor x 2^W / m \rfloor$, where $W$ is the width of `T`, is
    /// computed once, and each product is then reduced with a multiplication in place of a
    /// division. When $m \leq (2^W - 1) / 3$ the reductions are also lazy, leaving the value below
    /// $3m - 1$ until the end.
    ///
    /// # Worst-case complexity
    /// $T(n) = O(n)$
    ///
    /// $M(n) = O(1)$
    ///
    /// where $T$ is time, $M$ is additional memory, and $n$ is `self.len()`.
    ///
    /// # Panics
    /// Panics if `m` is 0, or if any coefficient of `self` or `x` is greater than or equal to `m`.
    ///
    /// # Examples
    /// ```
    /// use core::str::FromStr;
    /// use malachite_base::polynomial::ModEvaluate;
    /// use malachite_base::unsigned_polynomial::UnsignedPolynomial;
    ///
    /// let p = UnsignedPolynomial::<u8>::from_str("5*x^2+3*x+7").unwrap();
    /// // 5 * 36 + 3 * 6 + 7 = 205, which is 10 mod 13.
    /// assert_eq!((&p).mod_evaluate(6, 13), 10);
    /// assert_eq!((&p).mod_evaluate(0, 13), 7);
    /// // 205 itself, modulo a larger modulus.
    /// assert_eq!((&p).mod_evaluate(6, 211), 205);
    /// ```
    ///
    /// This is equivalent to `nmod_poly_evaluate_nmod` from `nmod_poly/evaluate_nmod.c`, FLINT
    /// 3.6.0, except that the value must be reduced.
    #[inline]
    fn mod_evaluate(self, x: T, m: T) -> T {
        mod_evaluate(self, x, m)
    }
}

impl<T: PrimitiveUnsigned> ModEvaluate<T> for UnsignedPolynomial<T> {
    type Output = T;

    /// Evaluates an [`UnsignedPolynomial`] at a value of its coefficient type, modulo a value of
    /// that type, taking the polynomial by value. The coefficients and the value must already be
    /// reduced modulo `m`.
    ///
    /// $$
    /// f(p, x, m) = \sum_{i=0}^{n-1} c_i x^i \bmod m,
    /// $$
    ///
    /// where $c_i$ is the coefficient of $x^i$ in $p$ and $n$ is its length. The zero polynomial
    /// evaluates to 0 everywhere.
    ///
    /// The value is found with Horner's rule, reducing after every step. When the polynomial is
    /// long enough and the top bit of `m` is clear, every multiplication is by the same `x`, so
    /// Shoup's method is used: $\lfloor x 2^W / m \rfloor$, where $W$ is the width of `T`, is
    /// computed once, and each product is then reduced with a multiplication in place of a
    /// division. When $m \leq (2^W - 1) / 3$ the reductions are also lazy, leaving the value below
    /// $3m - 1$ until the end.
    ///
    /// # Worst-case complexity
    /// $T(n) = O(n)$
    ///
    /// $M(n) = O(1)$
    ///
    /// where $T$ is time, $M$ is additional memory, and $n$ is `self.len()`.
    ///
    /// # Panics
    /// Panics if `m` is 0, or if any coefficient of `self` or `x` is greater than or equal to `m`.
    ///
    /// # Examples
    /// ```
    /// use core::str::FromStr;
    /// use malachite_base::polynomial::ModEvaluate;
    /// use malachite_base::unsigned_polynomial::UnsignedPolynomial;
    ///
    /// let p = UnsignedPolynomial::<u8>::from_str("5*x^2+3*x+7").unwrap();
    /// // 5 * 36 + 3 * 6 + 7 = 205, which is 10 mod 13.
    /// assert_eq!(p.clone().mod_evaluate(6, 13), 10);
    /// assert_eq!(p.clone().mod_evaluate(0, 13), 7);
    /// // 205 itself, modulo a larger modulus.
    /// assert_eq!(p.mod_evaluate(6, 211), 205);
    /// ```
    ///
    /// This is equivalent to `nmod_poly_evaluate_nmod` from `nmod_poly/evaluate_nmod.c`, FLINT
    /// 3.6.0, except that the value must be reduced.
    #[inline]
    fn mod_evaluate(self, x: T, m: T) -> T {
        mod_evaluate(&self, x, m)
    }
}

// The numbers of points evaluated together by `mod_evaluate_many_in_place`, for types at most 32
// bits wide and for wider types. Horner's rule is a chain of dependent multiplications, so running
// several points through one pass over the coefficients keeps more multiplications in flight. Tuned
// on Apple M-series: for `u64` a block of 8 is 3.4 to 4.8 times as fast as one point at a time at
// length 256, and for `u32` a block of 4 is best.
const MOD_EVALUATE_MANY_NARROW_BLOCK: usize = 4;
const MOD_EVALUATE_MANY_WIDE_BLOCK: usize = 8;

// Replaces each of the `N` points in `xs` with the polynomial's value there modulo `m`, with
// Horner's rule for all of them in one pass over the coefficients. `coefficients` must be nonempty,
// and the points and coefficients reduced.
crate_test_fn! {mod_evaluate_horner_block<T: PrimitiveUnsigned, const N: usize>(
    coefficients: &[T],
    xs: &mut [T; N],
    m: T,
) {
    let data = T::precompute_mod_mul_data(&m);
    let (&last, rest) = coefficients.split_last().unwrap();
    let points = *xs;
    let mut values = [last; N];
    for &c in rest.iter().rev() {
        for (value, &x) in values.iter_mut().zip(points.iter()) {
            value.mod_mul_precomputed_assign(x, m, &data);
            value.mod_add_assign(c, m);
        }
    }
    *xs = values;
}}

// Like `mod_evaluate_horner_block`, multiplying by each point with Shoup's method. The top bit of
// `m` must be clear.
crate_test_fn! {mod_evaluate_shoup_block<T: PrimitiveUnsigned, const N: usize>(
    coefficients: &[T],
    xs: &mut [T; N],
    m: T,
) {
    let points = *xs;
    let precomps = points.map(|x| mod_mul_precompute_shoup(x, m));
    let (&last, rest) = coefficients.split_last().unwrap();
    let mut values = [last; N];
    for &c in rest.iter().rev() {
        for ((value, &x), &x_precomp) in values.iter_mut().zip(points.iter()).zip(precomps.iter()) {
            *value = mod_mul_shoup(x, *value, x_precomp, m);
            value.mod_add_assign(c, m);
        }
    }
    *xs = values;
}}

// Like `mod_evaluate_shoup_block`, with lazy reduction as in `mod_evaluate_shoup_lazy`. `m` must be
// at most `T::MAX / 3`. The values are fully reduced at the end.
crate_test_fn! {mod_evaluate_shoup_lazy_block<T: PrimitiveUnsigned, const N: usize>(
    coefficients: &[T],
    xs: &mut [T; N],
    m: T,
) {
    let points = *xs;
    let precomps = points.map(|x| mod_mul_precompute_shoup(x, m));
    let (&last, rest) = coefficients.split_last().unwrap();
    let mut values = [last; N];
    for &c in rest.iter().rev() {
        for ((value, &x), &x_precomp) in values.iter_mut().zip(points.iter()).zip(precomps.iter()) {
            // value is x * value mod m, or that plus m, and then less than 3m - 1
            *value = mod_mul_shoup_lazy(x, *value, x_precomp, m);
            value.wrapping_add_assign(c);
        }
    }
    let two_m = m << 1;
    for value in &mut values {
        if *value >= two_m {
            *value -= two_m;
        } else if *value >= m {
            *value -= m;
        }
    }
    *xs = values;
}}

// Replaces each point in `xs` with the polynomial's value there modulo `m`. The points and the
// coefficients must be reduced; this is not checked. Blocks of points are evaluated together, with
// the method `mod_evaluate_slice` would choose for one point. Evaluates the points in `xs` in
// blocks of `N`, with the method chosen by `mod_evaluate_many_in_place`, and returns the leftover
// points, fewer than `N` of them.
fn mod_evaluate_many_blocks<'a, T: PrimitiveUnsigned, const N: usize>(
    coefficients: &[T],
    xs: &'a mut [T],
    m: T,
    shoup: bool,
    lazy: bool,
) -> &'a mut [T] {
    let (blocks, remainder) = xs.as_chunks_mut::<N>();
    for block in blocks {
        if !shoup {
            mod_evaluate_horner_block(coefficients, block, m);
        } else if lazy {
            mod_evaluate_shoup_lazy_block(coefficients, block, m);
        } else {
            mod_evaluate_shoup_block(coefficients, block, m);
        }
    }
    remainder
}

// Replaces each point in `xs` with the polynomial's value there modulo `m`. The points and the
// coefficients must be reduced; this is not checked. Blocks of points are evaluated together, with
// the method `mod_evaluate_slice` would choose for one point.
crate_test_fn! {mod_evaluate_many_in_place<T: PrimitiveUnsigned>(
    coefficients: &[T],
    xs: &mut [T],
    m: T,
) {
    let len = coefficients.len();
    match len {
        0 => xs.fill(T::ZERO),
        1 => xs.fill(coefficients[0]),
        _ => {
            // As in mod_evaluate_slice: Shoup's method needs the top bit of m clear, and is used
            // for short polynomials only when T is at most 32 bits wide
            let shoup = !m.get_highest_bit()
                && (T::WIDTH <= u32::WIDTH || len >= MOD_EVALUATE_SHOUP_THRESHOLD);
            let lazy = m <= T::MAX / T::from(3u8);
            // Wider types take blocks of the wide size first; the leftover points go through blocks
            // of the narrow size, and then one at a time
            let xs = if T::WIDTH <= u32::WIDTH {
                xs
            } else {
                mod_evaluate_many_blocks::<T, MOD_EVALUATE_MANY_WIDE_BLOCK>(
                    coefficients,
                    xs,
                    m,
                    shoup,
                    lazy,
                )
            };
            let remainder = mod_evaluate_many_blocks::<T, MOD_EVALUATE_MANY_NARROW_BLOCK>(
                coefficients,
                xs,
                m,
                shoup,
                lazy,
            );
            for x in remainder {
                *x = mod_evaluate_slice(coefficients, *x, m);
            }
        }
    }
}}

impl<T: PrimitiveUnsigned> ModEvaluateMany<T> for &UnsignedPolynomial<T> {
    type Output = T;

    /// Evaluates an [`UnsignedPolynomial`] at each of several values of its coefficient type,
    /// modulo a value of that type. The coefficients and the values must already be reduced modulo
    /// `m`.
    ///
    /// $$
    /// f(p, (x_j)_{j=0}^{k-1}, m) = \left ( \sum_{i=0}^{n-1} c_i x_j^i \bmod m
    /// \right )_{j=0}^{k-1},
    /// $$
    ///
    /// where $c_i$ is the coefficient of $x^i$ in $p$ and $n$ is its length.
    ///
    /// The result is the same as calling
    /// [`mod_evaluate`](crate::polynomial::ModEvaluate::mod_evaluate) at each value, with the same
    /// choice between Horner's rule and Shoup's method, but the polynomial is checked once, and
    /// several values are evaluated together in each pass over the coefficients. Horner's rule is a
    /// chain of dependent multiplications, so interleaving independent chains keeps the processor's
    /// multipliers busy.
    ///
    /// # Worst-case complexity
    /// $T(n, k) = O(nk)$
    ///
    /// $M(k) = O(k)$
    ///
    /// where $T$ is time, $M$ is additional memory, $n$ is `self.len()`, and $k$ is `xs.len()`.
    ///
    /// # Panics
    /// Panics if `m` is 0, or if any coefficient of `self` or any value in `xs` is greater than or
    /// equal to `m`.
    ///
    /// # Examples
    /// ```
    /// use core::str::FromStr;
    /// use malachite_base::polynomial::ModEvaluateMany;
    /// use malachite_base::unsigned_polynomial::UnsignedPolynomial;
    ///
    /// let p = UnsignedPolynomial::<u8>::from_str("5*x^2+3*x+7").unwrap();
    /// // 7, 15, 33, 61, 99, and 147, mod 13
    /// assert_eq!(
    ///     (&p).mod_evaluate_many(&[0, 1, 2, 3, 4, 5], 13),
    ///     &[7, 2, 7, 9, 8, 4]
    /// );
    /// ```
    ///
    /// This is equivalent to `nmod_poly_evaluate_nmod_vec_iter` from
    /// `nmod_poly/evaluate_nmod_vec.c`, FLINT 3.6.0, except that the values must be reduced.
    fn mod_evaluate_many(self, xs: &[T], m: T) -> Vec<T> {
        assert!(
            self.mod_is_reduced(&m),
            "self must be reduced mod m, but {self} has a coefficient >= {m}"
        );
        for &x in xs {
            assert!(x < m, "x must be reduced mod m, but {x} >= {m}");
        }
        let mut values = xs.to_vec();
        mod_evaluate_many_in_place(&self.coefficients, &mut values, m);
        values
    }
}

// Whether geometric evaluation of a polynomial of length `n` at `k` points modulo `m` uses
// Bluestein's trick rather than evaluating at each power of q: when the polynomial and the number
// of points are both long enough. Evaluating at each power is fast with Shoup's multiplication, and
// fastest with its lazy form, which needs the top two bits of `m` clear; when the top bit is set,
// Shoup's multiplication is unavailable and Bluestein's trick wins much sooner. The middle product
// accumulates in three words once `m` has more than about W - 8 bits, where W is `T::WIDTH`, which
// slows it. Measured on an Apple M-series machine, 2026-10, for `u64` with moduli of 20, 30, 40,
// 50, 62, 63, and 64 bits; the boundaries for other types are scaled by their width without being
// measured.
fn mod_evaluate_geometric_fast_preferred<T: PrimitiveUnsigned>(n: usize, k: usize, m: T) -> bool {
    let bits = m.significant_bits();
    let width = T::WIDTH;
    let (min_n, min_k) = if bits == width {
        (32, 32)
    } else if bits == width - 1 {
        (256, 256)
    } else if bits + 8 > width {
        (512, 1024)
    } else if bits << 4 > width * 5 {
        (256, 256)
    } else {
        (128, 128)
    };
    n >= min_n && k >= min_k
}

// Evaluates the polynomial with the given coefficients, reduced modulo `m`, at $1, q, q^2, \ldots,
// q^{k-1}$ modulo `m`, generating the powers of q and then evaluating at them all, several points
// at a time.
//
// This is equivalent to `_nmod_poly_evaluate_geometric_nmod_vec_iter` from
// `nmod_poly/evaluate_geometric_nmod_vec.c`, FLINT 3.6.0, with the ratio q in place of FLINT's
// $r^2$.
crate_test_fn! {mod_evaluate_geometric_iter<T: PrimitiveUnsigned>(
    coefficients: &[T],
    q: T,
    k: usize,
    m: T,
) -> Vec<T> {
    let mut values = Vec::with_capacity(k);
    if k != 0 {
        let mut power = T::ONE % m;
        values.push(power);
        if m.get_highest_bit() {
            let data = T::precompute_mod_mul_data(&m);
            for _ in 1..k {
                power.mod_mul_precomputed_assign(q, m, &data);
                values.push(power);
            }
        } else {
            let q_precomp = mod_mul_precompute_shoup(q, m);
            for _ in 1..k {
                power = mod_mul_shoup(q, power, q_precomp, m);
                values.push(power);
            }
        }
    }
    mod_evaluate_many_in_place(coefficients, &mut values, m);
    values
}}

// Evaluates the polynomial with the given coefficients, reduced modulo `m`, at $1, q, q^2, \ldots,
// q^{k-1}$ modulo `m` with Bluestein's trick, where `q_inverse` is the inverse of q modulo `m`.
//
// With $\binom{t}{2} = t(t-1)/2$, every product $ij$ is $\binom{i+j}{2} - \binom{i}{2} -
// \binom{j}{2}$, so
// $$
// \sum_i c_i q^{ij} = q^{-\binom{j}{2}} \sum_i \left ( c_i q^{-\binom{i}{2}} \right )
// q^{\binom{i+j}{2}}.
// $$
// The sum is coefficient $n - 1 + j$ of the product of the reverse of $(c_i q^{-\binom{i}{2}})_i$
// and $(q^{\binom{t}{2}})_t$, so all $k$ values are one middle product, the $k$ coefficients of
// that product from coefficient $n - 1$ on. Leading zero coefficients of the polynomial are
// skipped, shortening $n$.
//
// This is equivalent to `_nmod_poly_evaluate_geometric_nmod_vec_fast` from
// `nmod_poly/evaluate_geometric_nmod_vec.c`, FLINT 3.6.0, with exponents $\binom{t}{2}$ in place of
// FLINT's $t^2/2$, so that no square root of q is needed.
crate_test_fn! {mod_evaluate_geometric_fast<T: PrimitiveUnsigned>(
    coefficients: &[T],
    q: T,
    q_inverse: T,
    k: usize,
    m: T,
) -> Vec<T> {
    if k == 0 {
        return Vec::new();
    }
    let Some(start) = coefficients.iter().position(|&c| c != T::ZERO) else {
        return vec![T::ZERO; k];
    };
    let n = coefficients.len();
    let a_len = n - start;
    let data = T::precompute_mod_mul_data(&m);
    // binomial_powers[t] = q^C(t, 2), for t < n + k - 1.
    let b_len = n + k - 1;
    let mut binomial_powers = Vec::with_capacity(b_len);
    let mut power = T::ONE % m;
    let mut step = T::ONE % m;
    for _ in 0..b_len {
        binomial_powers.push(power);
        power.mod_mul_precomputed_assign(step, m, &data);
        step.mod_mul_precomputed_assign(q, m, &data);
    }
    // inverse_powers[t] = q^-C(t, 2), for t < max(n, k).
    let w_len = max(n, k);
    let mut inverse_powers = Vec::with_capacity(w_len);
    let mut power = T::ONE % m;
    let mut step = T::ONE % m;
    for _ in 0..w_len {
        inverse_powers.push(power);
        power.mod_mul_precomputed_assign(step, m, &data);
        step.mod_mul_precomputed_assign(q_inverse, m, &data);
    }
    // The scaled coefficients, from the first nonzero one on, reversed.
    let scaled: Vec<T> = coefficients[start..]
        .iter()
        .zip(&inverse_powers[start..])
        .rev()
        .map(|(&c, &w)| c.mod_mul_precomputed(w, m, &data))
        .collect();
    let mut product = vec![T::ZERO; k];
    mod_mul_middle_karatsuba(
        &mut product,
        &scaled,
        &binomial_powers[start..start + a_len + k - 1],
        &ModData::new(m, a_len),
    );
    product
        .iter()
        .zip(&inverse_powers)
        .map(|(&z, &w)| z.mod_mul_precomputed(w, m, &data))
        .collect()
}}

// Evaluates like `mod_evaluate_geometric_fast`, but computes the sums with a truncated product of
// length $n + k - 1$ and keeps its last $k$ coefficients, discarding the first $n - 1$, rather than
// with a middle product. It is kept to measure what the middle product saves.
crate_test_fn! {
#[allow(dead_code)]
mod_evaluate_geometric_fast_truncated<T: PrimitiveUnsigned>(
    coefficients: &[T],
    q: T,
    q_inverse: T,
    k: usize,
    m: T,
) -> Vec<T> {
    if k == 0 {
        return Vec::new();
    }
    let Some(start) = coefficients.iter().position(|&c| c != T::ZERO) else {
        return vec![T::ZERO; k];
    };
    let n = coefficients.len();
    let a_len = n - start;
    let data = T::precompute_mod_mul_data(&m);
    // binomial_powers[t] = q^C(t, 2), for t < n + k - 1.
    let b_len = n + k - 1;
    let mut binomial_powers = Vec::with_capacity(b_len);
    let mut power = T::ONE % m;
    let mut step = T::ONE % m;
    for _ in 0..b_len {
        binomial_powers.push(power);
        power.mod_mul_precomputed_assign(step, m, &data);
        step.mod_mul_precomputed_assign(q, m, &data);
    }
    // inverse_powers[t] = q^-C(t, 2), for t < max(n, k).
    let w_len = max(n, k);
    let mut inverse_powers = Vec::with_capacity(w_len);
    let mut power = T::ONE % m;
    let mut step = T::ONE % m;
    for _ in 0..w_len {
        inverse_powers.push(power);
        power.mod_mul_precomputed_assign(step, m, &data);
        step.mod_mul_precomputed_assign(q_inverse, m, &data);
    }
    // The scaled coefficients, from the first nonzero one on, reversed.
    let scaled: Vec<T> = coefficients[start..]
        .iter()
        .zip(&inverse_powers[start..])
        .rev()
        .map(|(&c, &w)| c.mod_mul_precomputed(w, m, &data))
        .collect();
    let mut product = vec![T::ZERO; a_len + k - 1];
    mod_mul_truncated_to_out(
        &mut product,
        &scaled,
        &binomial_powers[start..start + a_len + k - 1],
        m,
    );
    product[a_len - 1..]
        .iter()
        .zip(&inverse_powers)
        .map(|(&z, &w)| z.mod_mul_precomputed(w, m, &data))
        .collect()
}}

impl<T: PrimitiveUnsigned> ModEvaluateGeometric<T> for &UnsignedPolynomial<T> {
    type Output = T;

    /// Evaluates an [`UnsignedPolynomial`] at $1, q, q^2, \ldots, q^{k-1}$, modulo a value of its
    /// coefficient type. The coefficients and `q` must already be reduced modulo `m`.
    ///
    /// $$
    /// f(p, q, k, m) = \left ( \sum_{i=0}^{n-1} c_i q^{ij} \bmod m \right )_{j=0}^{k-1},
    /// $$
    ///
    /// where $c_i$ is the coefficient of $x^i$ in $p$ and $n$ is its length.
    ///
    /// The powers of `q` are computed with Shoup's method when the top bit of `m` is clear, since
    /// every multiplication is by `q`, and are then evaluated as by
    /// [`mod_evaluate_many`](crate::polynomial::ModEvaluateMany::mod_evaluate_many), in place.
    ///
    /// # Worst-case complexity
    /// $T(n, k) = O(nk)$
    ///
    /// $M(k) = O(k)$
    ///
    /// where $T$ is time, $M$ is additional memory, $n$ is `self.len()`, and $k$ is `k`.
    ///
    /// # Panics
    /// Panics if `m` is 0, or if any coefficient of `self` or `q` is greater than or equal to `m`.
    ///
    /// # Examples
    /// ```
    /// use core::str::FromStr;
    /// use malachite_base::polynomial::ModEvaluateGeometric;
    /// use malachite_base::unsigned_polynomial::UnsignedPolynomial;
    ///
    /// let p = UnsignedPolynomial::<u8>::from_str("5*x^2+3*x+7").unwrap();
    /// // At 1, 2, 4, and 8: 15, 33, 99, and 351, mod 13
    /// assert_eq!((&p).mod_evaluate_geometric(2, 4, 13), &[2, 7, 8, 0]);
    /// ```
    ///
    /// This is equivalent to `nmod_poly_evaluate_geometric_nmod_vec_iter` from
    /// `nmod_poly/evaluate_geometric_nmod_vec.c`, FLINT 3.6.0, with `q` in place of FLINT's $r^2$:
    /// FLINT evaluates at the powers of the square of its argument.
    fn mod_evaluate_geometric(self, q: T, k: u64, m: T) -> Vec<T> {
        assert!(
            self.mod_is_reduced(&m),
            "self must be reduced mod m, but {self} has a coefficient >= {m}"
        );
        assert!(q < m, "q must be reduced mod m, but {q} >= {m}");
        let k = usize::exact_from(k);
        let n = self.coefficients.len();
        if mod_evaluate_geometric_fast_preferred(n, k, m)
            && q != T::ZERO
            && let Some(q_inverse) = q.mod_inverse(m)
        {
            mod_evaluate_geometric_fast(&self.coefficients, q, q_inverse, k, m)
        } else {
            mod_evaluate_geometric_iter(&self.coefficients, q, k, m)
        }
    }
}