Trait NttReuse 

pub trait NttReuse: ConvolveSteps {
    const MULTIPLE: bool = true;

    // Required methods
    fn ntt_doubling(f: Self::F, monic: bool) -> Self::F;
    fn even_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F;
    fn odd_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F;
    fn multiply_prefix(f: &mut Self::F, g: &Self::F);
    fn multiply_add(sum: &mut Self::F, f: &Self::F, g: &Self::F);

    // Provided methods
    fn transform_ntt(t: Self::T, len: usize) -> Self::F { ... }
    fn inverse_transform_ntt(f: Self::F, len: usize) -> Self::T { ... }
    fn max_product_sum_count(_f: &Self::F) -> usize { ... }
    fn power_projection_step(
        p_flat: Self::T,
        q_flat: Self::T,
        n: usize,
        py: usize,
        qy: usize,
    ) -> (Self::T, Self::T) { ... }
}

Provided Associated Constants

const MULTIPLE: bool = true

Required Methods

fn ntt_doubling(f: Self::F, monic: bool) -> Self::F

Extends a value produced by transform_ntt to twice its length. If monic, the input represents a monic degree-n polynomial modulo x^n - 1, where n is the transform length.

fn even_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F

Extracts the even coefficients of a(x) * b(-x) in the usual NTT frequency order.

fn odd_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F

Extracts the odd coefficients of a(x) * b(-x) in the usual NTT frequency order.

fn multiply_prefix(f: &mut Self::F, g: &Self::F)

Multiplies a usual NTT transform by the corresponding prefix of another one.

fn multiply_add(sum: &mut Self::F, f: &Self::F, g: &Self::F)

Adds the pointwise product of two usual NTT transforms to sum.

Provided Methods

fn transform_ntt(t: Self::T, len: usize) -> Self::F

Transforms coefficients into the usual NTT frequency order.

Examples found in repository

crates/competitive/src/math/formal_power_series/berlekamp_massey.rs (line 345)

fn reduced_transform<T, C>(fps: &FormalPowerSeries<T, C>, length: usize) -> C::F
where
    T: FormalPowerSeriesCoefficient,
    C: NttReuse<T = Vec<T>>,
{
    let mut coefficients = vec![T::zero(); length];
    for (i, value) in fps.iter().enumerate() {
        coefficients[i & (length - 1)] += value;
    }
    C::transform_ntt(coefficients, length)
}

fn transform_window<T, C>(fps: &FormalPowerSeries<T, C>, end: isize, length: usize) -> C::F
where
    T: FormalPowerSeriesCoefficient,
    C: NttReuse<T = Vec<T>>,
{
    let start = end - length as isize;
    let coefficients = (start..end).map(|index| coefficient(fps, index)).collect();
    C::transform_ntt(coefficients, length)
}

crates/competitive/src/math/number_theoretic_transform.rs (line 593)

fn convolve_large_ntt<M>(a: Vec<MInt<M>>, b: Vec<MInt<M>>) -> Vec<MInt<M>>
where
    M: Montgomery32NttModulus,
{
    let len = a.len() + b.len() - 1;
    let ntt_len = 1usize << M::RANK;
    let block_len = ntt_len / 2;
    let same = a == b;
    let transform = |a: &[MInt<M>]| {
        let mut f = Vec::with_capacity(ntt_len);
        advise_huge_pages(&mut f);
        f.extend_from_slice(a);
        Convolve::<M>::transform_ntt(f, ntt_len)
    };
    let fa: Vec<_> = a.chunks(block_len).map(transform).collect();
    let fb: Option<Vec<_>> = if same {
        None
    } else {
        Some(b.chunks(block_len).map(transform).collect())
    };
    let b_blocks = fb.as_ref().map_or(fa.len(), Vec::len);
    let mut result = vec![MInt::<M>::zero(); len];
    for diagonal in 0..fa.len() + b_blocks - 1 {
        let mut spectrum = vec![MInt::<M>::zero(); ntt_len];
        let start = diagonal.saturating_sub(b_blocks - 1);
        for i in start..=diagonal.min(fa.len() - 1) {
            let j = diagonal - i;
            let g = if let Some(fb) = &fb { &fb[j] } else { &fa[j] };
            pointwise_multiply_add(&mut spectrum, &fa[i], g);
        }
        spectrum = Convolve::<M>::inverse_transform_ntt(spectrum, ntt_len);
        let offset = diagonal * block_len;
        for (result, value) in result[offset..].iter_mut().zip(spectrum) {
            *result += value;
        }
    }
    result
}

impl<M> ConvolveSteps for Convolve<M>
where
    M: Montgomery32NttModulus,
{
    const CYCLIC: bool = true;

    type T = Vec<MInt<M>>;
    type F = Vec<MInt<M>>;
    fn length(t: &Self::T) -> usize {
        t.len()
    }
    fn transform(mut t: Self::T, len: usize) -> Self::F {
        t.resize_with(len.max(1).next_power_of_two(), Zero::zero);
        #[cfg(target_arch = "x86_64")]
        if use_block_ntt::<M>(t.len()) {
            unsafe { ntt_simd::transform_blocks_avx2(&mut t) };
            return t;
        }
        ntt(&mut t);
        t
    }
    fn inverse_transform(mut f: Self::F, len: usize) -> Self::T {
        #[cfg(target_arch = "x86_64")]
        if use_block_ntt::<M>(f.len()) {
            unsafe { ntt_simd::inverse_transform_blocks_avx2(&mut f) };
            f.truncate(len);
            return f;
        }
        intt(&mut f);
        f.truncate(len);
        f
    }
    fn multiply(f: &mut Self::F, g: &Self::F) {
        assert_eq!(f.len(), g.len());
        #[cfg(target_arch = "x86_64")]
        if use_block_ntt::<M>(f.len()) {
            unsafe { ntt_simd::multiply_blocks_avx2(f, g) };
            return;
        }
        pointwise_multiply(f, g);
    }
    fn square(t: Self::T, len: usize) -> Self::T {
        let mut f = Self::transform(t, len);
        let g = f.clone();
        Self::multiply(&mut f, &g);
        Self::inverse_transform(f, len)
    }
    fn convolve(mut a: Self::T, mut b: Self::T) -> Self::T {
        let (threshold, naive_threshold) = (100, 60);
        #[cfg(target_arch = "x86_64")]
        let (threshold, naive_threshold) = if use_block_ntt::<M>(64) {
            (
                60,
                if M::RANK >= 13 && a.len().max(b.len()) <= 4096 {
                    18
                } else if M::RANK >= 19 && a.len().max(b.len()) <= 262144 {
                    32
                } else {
                    34
                },
            )
        } else {
            (threshold, naive_threshold)
        };
        if Self::length(&a).max(Self::length(&b)) <= threshold {
            return convolve_karatsuba(&a, &b);
        }
        if Self::length(&a).min(Self::length(&b)) <= naive_threshold {
            return convolve_naive(&a, &b);
        }
        let len = (Self::length(&a) + Self::length(&b)).saturating_sub(1);
        let size = len.max(1).next_power_of_two();
        let max_size = 1usize << M::RANK;
        #[cfg(target_arch = "x86_64")]
        let max_size = if use_block_ntt::<M>(size) {
            max_size << 3
        } else {
            max_size
        };
        if size > max_size {
            return convolve_large_ntt(a, b);
        }
        if len <= size / 2 + 2 {
            let xa = a.pop().unwrap();
            let xb = b.pop().unwrap();
            let mut c = vec![MInt::<M>::zero(); len];
            *c.last_mut().unwrap() = xa * xb;
            for (a, c) in a.iter().zip(&mut c[b.len()..]) {
                *c += *a * xb;
            }
            for (b, c) in b.iter().zip(&mut c[a.len()..]) {
                *c += *b * xa;
            }
            let d = Self::convolve(a, b);
            for (d, c) in d.into_iter().zip(&mut c) {
                *c += d;
            }
            return c;
        }
        let same = a == b;
        #[cfg(target_arch = "x86_64")]
        if use_block_ntt::<M>(size) {
            a.reserve(size - a.len());
            b.reserve(size - b.len());
            advise_huge_pages(&mut a);
            advise_huge_pages(&mut b);
            a.resize_with(size, Zero::zero);
            b.resize_with(size, Zero::zero);
            unsafe { ntt_simd::convolve_blocks_avx2(&mut a, &mut b, same) };
            a.truncate(len);
            return a;
        }
        let mut a = Self::transform(a, len);
        if same {
            for a in a.iter_mut() {
                *a *= *a;
            }
        } else {
            let b = Self::transform(b, len);
            Self::multiply(&mut a, &b);
        }
        Self::inverse_transform(a, len)
    }
}

type MVec<M> = Vec<MInt<M>>;

fn convert_crt_input<M, N1, N2, N3>(t: MVec<M>, capacity: usize) -> (MVec<N1>, MVec<N2>, MVec<N3>)
where
    M: MIntConvert<u32>,
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    let mut f = (
        MVec::<N1>::with_capacity(capacity),
        MVec::<N2>::with_capacity(capacity),
        MVec::<N3>::with_capacity(capacity),
    );
    advise_huge_pages(&mut f.0);
    advise_huge_pages(&mut f.1);
    advise_huge_pages(&mut f.2);
    for t in t {
        let t: u32 = t.into();
        f.0.push(t.into());
        f.1.push(t.into());
        f.2.push(t.into());
    }
    f
}

fn reconstruct_mint_crt<M, N1, N2, N3>(f: (MVec<N1>, MVec<N2>, MVec<N3>)) -> MVec<M>
where
    M: MIntConvert + MIntConvert<u32>,
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    let t1 = MInt::<N2>::new(N1::get_mod()).inv();
    let m1_3 = MInt::<N3>::new(N1::get_mod());
    let t2 = (m1_3 * MInt::<N3>::new(N2::get_mod())).inv();
    let modulus = <M as MIntConvert<u32>>::mod_into() as u64;
    let m1 = N1::get_mod() as u64;
    let m2 = m1 * N2::get_mod() as u64 % modulus;
    let fits_u64 = (N1::get_mod() - 1) as u128
        + (N2::get_mod() - 1) as u128 * m1 as u128
        + (N3::get_mod() - 1) as u128 * m2 as u128
        <= u64::MAX as u128;
    f.0.into_iter()
        .zip(f.1)
        .zip(f.2)
        .map(|((c1, c2), c3)| {
            let d1 = c1.inner();
            let d2 = ((c2 - MInt::<N2>::from(d1)) * t1).inner();
            let x = MInt::<N3>::new(d1) + MInt::<N3>::new(d2) * m1_3;
            let d3 = ((c3 - x) * t2).inner();
            let value = if fits_u64 {
                (d1 as u64 + d2 as u64 * m1 + d3 as u64 * m2) % modulus
            } else {
                ((d1 as u128 + d2 as u128 * m1 as u128 + d3 as u128 * m2 as u128) % modulus as u128)
                    as u64
            };
            MInt::<M>::from(value as u32)
        })
        .collect()
}

impl<M, N1, N2, N3> ConvolveSteps for Convolve<(M, (N1, N2, N3))>
where
    M: MIntConvert + MIntConvert<u32>,
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    type T = MVec<M>;
    type F = (MVec<N1>, MVec<N2>, MVec<N3>);
    fn length(t: &Self::T) -> usize {
        t.len()
    }
    fn transform(t: Self::T, len: usize) -> Self::F {
        let npot = len.max(1).next_power_of_two();
        let f = convert_crt_input(t, npot);
        (
            Convolve::<N1>::transform(f.0, npot),
            Convolve::<N2>::transform(f.1, npot),
            Convolve::<N3>::transform(f.2, npot),
        )
    }
    fn inverse_transform(f: Self::F, len: usize) -> Self::T {
        reconstruct_mint_crt((
            Convolve::<N1>::inverse_transform(f.0, len),
            Convolve::<N2>::inverse_transform(f.1, len),
            Convolve::<N3>::inverse_transform(f.2, len),
        ))
    }
    fn multiply(f: &mut Self::F, g: &Self::F) {
        Convolve::<N1>::multiply(&mut f.0, &g.0);
        Convolve::<N2>::multiply(&mut f.1, &g.1);
        Convolve::<N3>::multiply(&mut f.2, &g.2);
    }
    fn convolve(a: Self::T, b: Self::T) -> Self::T {
        let max_len = Self::length(&a).max(Self::length(&b));
        let min_len = Self::length(&a).min(Self::length(&b));
        let (balanced, short) = crate::avx_helper!(@dispatch_avx2_fma (30, 10), (384, 128));
        if max_len <= balanced || min_len <= short {
            return convolve_karatsuba(&a, &b);
        }
        // Limit coefficient growth to leave headroom for FFT roundoff.
        let fft_limit = crate::avx_helper!(@dispatch_avx2_fma
            1usize << ((1u64 << 50) / <M as MIntConvert<u32>>::mod_into() as u64).ilog2().min(20), 0);
        let convolve = |a: Self::T, b: Self::T| {
            let fft_len = (a.len() + b.len() - 1).next_power_of_two();
            if fft_len <= 256 && a.len() * b.len() <= fft_len * 8 {
                return convolve_karatsuba(&a, &b);
            }
            if fft_len <= fft_limit {
                crate::avx_helper!(@dispatch_avx2_fma return unsafe {
                    convolve_mint_avx2(a, b)
                }, ());
            }
            convolve_mint_crt::<M, N1, N2, N3>(a, b)
        };
        let block_len = min_len.next_power_of_two() * 8 - min_len + 1;
        let block_len = if min_len <= fft_limit / 2 {
            block_len.min(fft_limit - min_len + 1)
        } else {
            block_len
        };
        if max_len <= block_len {
            return convolve(a, b);
        }
        let (a, b) = if a.len() >= b.len() { (a, b) } else { (b, a) };
        let mut result = vec![MInt::<M>::zero(); a.len() + b.len() - 1];
        for (i, a) in a.chunks(block_len).enumerate() {
            let product = convolve(a.to_vec(), b.clone());
            for (value, product) in result[i * block_len..].iter_mut().zip(product) {
                *value += product;
            }
        }
        result
    }
}

fn convolve_mint_crt<M, N1, N2, N3>(a: MVec<M>, b: MVec<M>) -> MVec<M>
where
    M: MIntConvert + MIntConvert<u32>,
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    let convolve = |a: MVec<M>, b: MVec<M>| {
        let a_len = a.len();
        let b_len = b.len();
        let a = convert_crt_input(a, a_len);
        let b = convert_crt_input(b, b_len);
        reconstruct_mint_crt((
            Convolve::<N1>::convolve(a.0, b.0),
            Convolve::<N2>::convolve(a.1, b.1),
            Convolve::<N3>::convolve(a.2, b.2),
        ))
    };
    let modulus = <M as MIntConvert<u32>>::mod_into() as u128;
    let capacity = N1::MOD as u128 * N2::MOD as u128 * N3::MOD as u128;
    if a.len().min(b.len()) as u128 * (modulus - 1).pow(2) < capacity {
        return convolve(a, b);
    }
    let block_len = ((capacity - 1) / (modulus - 1).pow(2)) as usize;
    if block_len == 0 {
        return convolve_naive(&a, &b);
    }
    let mut result = vec![MInt::<M>::zero(); a.len() + b.len() - 1];
    for (i, a) in a.chunks(block_len).enumerate() {
        for (j, b) in b.chunks(block_len).enumerate() {
            let product = convolve(a.to_vec(), b.to_vec());
            for (value, product) in result[(i + j) * block_len..].iter_mut().zip(product) {
                *value += product;
            }
        }
    }
    result
}

impl<N1, N2, N3> ConvolveSteps for Convolve<(u64, (N1, N2, N3))>
where
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    type T = Vec<u64>;
    type F = ([MVec<N1>; 3], [MVec<N2>; 3], [MVec<N3>; 3]);

    fn length(t: &Self::T) -> usize {
        t.len()
    }

    fn transform(t: Self::T, len: usize) -> Self::F {
        let npot = len.max(1).next_power_of_two();
        assert!(npot <= 1usize << N1::RANK.min(N2::RANK).min(N3::RANK));
        // The 22-bit fallback needs room for three limb products per coefficient.
        assert!(
            3 * npot as u128 * ((1u128 << 22) - 1).pow(2)
                < N1::MOD as u128 * N2::MOD as u128 * N3::MOD as u128
        );
        let bits = if 2 * npot as u128 * (u32::MAX as u128).pow(2)
            < N1::MOD as u128 * N2::MOD as u128 * N3::MOD as u128
        {
            32
        } else {
            22
        };
        let parts = if bits == 32 && t.iter().all(|&value| value <= u32::MAX as u64) {
            1
        } else {
            64usize.div_ceil(bits)
        };
        fn split<M: Montgomery32NttModulus>(
            t: &[u64],
            len: usize,
            bits: usize,
            parts: usize,
        ) -> [MVec<M>; 3] {
            std::array::from_fn(|part| {
                if part >= parts {
                    return Vec::new();
                }
                Convolve::<M>::transform(
                    t.iter()
                        .map(|&t| MInt::from((t >> (part * bits)) & ((1u64 << bits) - 1)))
                        .collect(),
                    len,
                )
            })
        }
        (
            split(&t, npot, bits, parts),
            split(&t, npot, bits, parts),
            split(&t, npot, bits, parts),
        )
    }

    fn inverse_transform(f: Self::F, len: usize) -> Self::T {
        let bits = if f.0[2].is_empty() { 32 } else { 22 };
        let t1 = MInt::<N2>::new(N1::get_mod()).inv();
        let m1 = N1::get_mod() as u64;
        let m1_3 = MInt::<N3>::new(N1::get_mod());
        let t2 = (m1_3 * MInt::<N3>::new(N2::get_mod())).inv();
        let m2 = m1 * N2::get_mod() as u64;
        let mut result = vec![0u64; len.min(f.0[0].len())];
        for (part, ((f1, f2), f3)) in f.0.into_iter().zip(f.1).zip(f.2).enumerate() {
            if f1.is_empty() {
                continue;
            }
            for (value, ((c1, c2), c3)) in result.iter_mut().zip(
                Convolve::<N1>::inverse_transform(f1, len)
                    .into_iter()
                    .zip(Convolve::<N2>::inverse_transform(f2, len))
                    .zip(Convolve::<N3>::inverse_transform(f3, len)),
            ) {
                let d1 = c1.inner();
                let d2 = ((c2 - MInt::<N2>::from(d1)) * t1).inner();
                let x = MInt::<N3>::new(d1) + MInt::<N3>::new(d2) * m1_3;
                let d3 = ((c3 - x) * t2).inner();
                let limb = (d1 as u64)
                    .wrapping_add((d2 as u64).wrapping_mul(m1))
                    .wrapping_add((d3 as u64).wrapping_mul(m2));
                *value = value.wrapping_add(limb << (part * bits));
            }
        }
        result
    }

    fn multiply(f: &mut Self::F, g: &Self::F) {
        fn multiply<M: Montgomery32NttModulus>(f: &mut [MVec<M>; 3], g: &[MVec<M>; 3]) {
            assert_eq!(f[0].len(), g[0].len());
            if f[1].is_empty() || g[1].is_empty() {
                if f[1].is_empty() && !g[1].is_empty() {
                    f[1] = f[0].clone();
                    Convolve::<M>::multiply(&mut f[1], &g[1]);
                } else if !f[1].is_empty() {
                    Convolve::<M>::multiply(&mut f[1], &g[0]);
                }
                Convolve::<M>::multiply(&mut f[0], &g[0]);
                return;
            }
            #[cfg(target_arch = "x86_64")]
            if use_block_ntt::<M>(f[0].len()) {
                for part in (1..if f[2].is_empty() { 2 } else { 3 }).rev() {
                    let mut sum = f[0].clone();
                    Convolve::<M>::multiply(&mut sum, &g[part]);
                    for left in 1..=part {
                        let mut product = f[left].clone();
                        Convolve::<M>::multiply(&mut product, &g[part - left]);
                        for (value, product) in sum.iter_mut().zip(product) {
                            // Block products contain lazy Montgomery residues.
                            *value = MInt::new(value.inner() + product.inner());
                        }
                    }
                    f[part] = sum;
                }
                Convolve::<M>::multiply(&mut f[0], &g[0]);
                return;
            }
            if f[2].is_empty() {
                for i in 0..f[0].len() {
                    f[1][i] = f[0][i] * g[1][i] + f[1][i] * g[0][i];
                    f[0][i] *= g[0][i];
                }
                return;
            }
            for i in 0..f[0].len() {
                f[2][i] = f[0][i] * g[2][i] + f[1][i] * g[1][i] + f[2][i] * g[0][i];
                f[1][i] = f[0][i] * g[1][i] + f[1][i] * g[0][i];
                f[0][i] *= g[0][i];
            }
        }
        multiply(&mut f.0, &g.0);
        multiply(&mut f.1, &g.1);
        multiply(&mut f.2, &g.2);
    }

    fn square(t: Self::T, len: usize) -> Self::T {
        let mut f = Self::transform(t, len);
        let g = f.clone();
        Self::multiply(&mut f, &g);
        Self::inverse_transform(f, len)
    }

    fn convolve(a: Self::T, b: Self::T) -> Self::T {
        let max_len = Self::length(&a).max(Self::length(&b));
        let min_len = Self::length(&a).min(Self::length(&b));
        let (balanced, short) = crate::avx_helper!(@dispatch_avx2_fma (300, 64), (1536, 512));
        if max_len <= balanced || min_len <= short {
            let a_wrapping: &[Wrapping<u64>] =
                unsafe { std::slice::from_raw_parts(a.as_ptr().cast(), a.len()) };
            let b_wrapping: &[Wrapping<u64>] =
                unsafe { std::slice::from_raw_parts(b.as_ptr().cast(), b.len()) };
            let mut c = std::mem::ManuallyDrop::new(if max_len <= 300 || min_len > 60 {
                convolve_karatsuba(a_wrapping, b_wrapping)
            } else {
                convolve_naive(a_wrapping, b_wrapping)
            });
            return unsafe { Vec::from_raw_parts(c.as_mut_ptr().cast(), c.len(), c.capacity()) };
        }
        let len = (Self::length(&a) + Self::length(&b)).saturating_sub(1);
        let block_len = if min_len >= 1 << 20 {
            1 << 20
        } else {
            (min_len.next_power_of_two() * 8).min(1 << 21) - min_len + 1
        };
        if max_len <= block_len {
            return convolve_u64_fft(a, b);
        }
        let mut result = vec![0u64; len];
        for (i, a) in a.chunks(block_len).enumerate() {
            for (j, b) in b.chunks(block_len).enumerate() {
                if a.len().min(b.len()) <= 60 {
                    for (x, &a) in a.iter().enumerate() {
                        for (y, &b) in b.iter().enumerate() {
                            let value = &mut result[(i + j) * block_len + x + y];
                            *value = value.wrapping_add(a.wrapping_mul(b));
                        }
                    }
                    continue;
                }
                let product = convolve_u64_fft(a.to_vec(), b.to_vec());
                for (value, product) in result[(i + j) * block_len..].iter_mut().zip(product) {
                    *value = value.wrapping_add(product);
                }
            }
        }
        result
    }
}

fn convolve_u64_fft(a: Vec<u64>, b: Vec<u64>) -> Vec<u64> {
    // Keep limb convolutions below 2^47 at the 2^21 FFT limit.
    crate::avx_helper!(@dispatch_avx2_fma return unsafe {
        convolve_u64_avx2(a, b)
    }, ());
    convolve_u64_fft_scalar(a, b)
}

fn convolve_u64_fft_scalar(a: Vec<u64>, b: Vec<u64>) -> Vec<u64> {
    fn split(values: &[u64]) -> [Vec<i64>; 5] {
        let mut result = std::array::from_fn(|_| Vec::with_capacity(values.len()));
        for mut value in values.iter().copied() {
            for part in &mut result {
                let digit = ((value << 51) as i64) >> 51;
                part.push(digit);
                value = (value >> 13).wrapping_add(u64::from(digit < 0));
            }
        }
        result
    }

    let len = a.len() + b.len() - 1;
    let transform = |values: &[u64]| {
        if values.iter().any(|&value| value > u32::MAX as u64) {
            return split(values).map(|part| ConvolveRealFft::transform(part, len));
        }
        let [a, b, c, _, _] = split(values);
        let a = ConvolveRealFft::transform(a, len);
        let size = a.len();
        [
            a,
            ConvolveRealFft::transform(b, len),
            ConvolveRealFft::transform(c, len),
            vec![Zero::zero(); size],
            vec![Zero::zero(); size],
        ]
    };
    let fa = transform(&a);
    drop(a);
    let fb = transform(&b);
    drop(b);
    let values: [Vec<i64>; 5] = std::array::from_fn(|part| {
        let mut sum = fa[0].clone();
        ConvolveRealFft::multiply(&mut sum, &fb[part]);
        for left in 1..=part {
            let mut product = fa[left].clone();
            ConvolveRealFft::multiply(&mut product, &fb[part - left]);
            for (sum, product) in sum.iter_mut().zip(product) {
                *sum += product;
            }
        }
        ConvolveRealFft::inverse_transform(sum, len)
    });
    (0..len)
        .map(|i| {
            (values[0][i] as u64)
                .wrapping_add((values[1][i] as u64) << 13)
                .wrapping_add((values[2][i] as u64) << 26)
                .wrapping_add((values[3][i] as u64) << 39)
                .wrapping_add((values[4][i] as u64) << 52)
        })
        .collect()
}

pub trait NttReuse: ConvolveSteps {
    const MULTIPLE: bool = true;

    /// Transforms coefficients into the usual NTT frequency order.
    fn transform_ntt(t: Self::T, len: usize) -> Self::F {
        Self::transform(t, len)
    }

    /// Inverts a value produced by `transform_ntt`.
    fn inverse_transform_ntt(f: Self::F, len: usize) -> Self::T {
        Self::inverse_transform(f, len)
    }

    /// Extends a value produced by `transform_ntt` to twice its length.
    /// If `monic`, the input represents a monic degree-`n` polynomial modulo
    /// `x^n - 1`, where `n` is the transform length.
    fn ntt_doubling(f: Self::F, monic: bool) -> Self::F;

    /// Extracts the even coefficients of `a(x) * b(-x)` in the usual NTT frequency order.
    fn even_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F;

    /// Extracts the odd coefficients of `a(x) * b(-x)` in the usual NTT frequency order.
    fn odd_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F;

    /// Multiplies a usual NTT transform by the corresponding prefix of another one.
    fn multiply_prefix(f: &mut Self::F, g: &Self::F);

    /// Adds the pointwise product of two usual NTT transforms to `sum`.
    fn multiply_add(sum: &mut Self::F, f: &Self::F, g: &Self::F);

    /// Maximum number of products that can be summed before reconstruction.
    /// Both factors must transform canonical coefficients at the supplied transform's length,
    /// and each cyclic product must itself be reconstructible.
    fn max_product_sum_count(_f: &Self::F) -> usize {
        if Self::MULTIPLE { 1 } else { usize::MAX }
    }

    fn power_projection_step(
        p_flat: Self::T,
        q_flat: Self::T,
        n: usize,
        py: usize,
        qy: usize,
    ) -> (Self::T, Self::T) {
        let base = n * 2;
        let len_p = base * py;
        let len_q = base * qy;
        let len = (len_p + len_q - 1).max(len_q + len_q - 1);
        let half = len.max(1).next_power_of_two() / 2;

        let p_fft = Self::transform_ntt(p_flat, len);
        let q_fft = Self::transform_ntt(q_flat, len);
        let pr_fft = Self::odd_mul_normal_neg(&p_fft, &q_fft);
        let qr_fft = Self::even_mul_normal_neg(&q_fft, &q_fft);
        (
            Self::inverse_transform_ntt(pr_fft, half),
            Self::inverse_transform_ntt(qr_fft, half),
        )
    }
}

thread_local!(
    static BIT_REVERSE: UnsafeCell<Vec<Vec<usize>>> = const { UnsafeCell::new(vec![]) };
);

impl<M> NttReuse for Convolve<M>
where
    M: Montgomery32NttModulus,
{
    const MULTIPLE: bool = false;

    fn transform_ntt(mut t: Self::T, len: usize) -> Self::F {
        t.resize_with(len.max(1).next_power_of_two(), Zero::zero);
        ntt(&mut t);
        t
    }

    fn inverse_transform_ntt(mut f: Self::F, len: usize) -> Self::T {
        intt(&mut f);
        f.truncate(len);
        f
    }

    fn ntt_doubling(mut f: Self::F, monic: bool) -> Self::F {
        let n = f.len();
        let k = n.trailing_zeros() as usize;
        let mut a = Self::inverse_transform_ntt(f.clone(), n);
        if monic {
            a[0] -= MInt::<M>::from(2);
        }
        let zeta = MInt::<M>::new_unchecked(M::INFO.root[k + 1]);
        let zeta2 = zeta * zeta;
        let mut rot = [MInt::one(), zeta, zeta2, zeta2 * zeta];
        let step = zeta2 * zeta2;
        for a in a.chunks_mut(4) {
            for (a, rot) in a.iter_mut().zip(&mut rot) {
                *a *= *rot;
                *rot *= step;
            }
        }
        f.extend(Self::transform_ntt(a, n));
        f
    }

    fn even_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F {
        assert_eq!(f.len(), g.len());
        assert!(f.len().is_power_of_two());
        assert!(f.len() >= 2);
        if std::ptr::eq(f, g) {
            return f.as_chunks::<2>().0.iter().map(|a| a[0] * a[1]).collect();
        }
        let inv2 = MInt::<M>::from(2).inv();
        let n = f.len() / 2;
        (0..n)
            .map(|i| (f[i << 1] * g[i << 1 | 1] + f[i << 1 | 1] * g[i << 1]) * inv2)
            .collect()
    }

    fn odd_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F {
        assert_eq!(f.len(), g.len());
        assert!(f.len().is_power_of_two());
        assert!(f.len() >= 2);
        let mut inv2 = MInt::<M>::from(2).inv();
        let n = f.len() / 2;
        let k = f.len().trailing_zeros() as usize;
        let mut h = vec![MInt::<M>::zero(); n];
        let w = MInt::<M>::new_unchecked(M::INFO.inv_root[k]);
        BIT_REVERSE.with(|br| {
            let br = unsafe { &mut *br.get() };
            if br.len() < k {
                br.resize_with(k, Default::default);
            }
            let k = k - 1;
            if br[k].is_empty() {
                let mut v = vec![0; 1 << k];
                for i in 0..1 << k {
                    v[i] = (v[i >> 1] >> 1) | ((i & 1) << k.saturating_sub(1));
                }
                br[k] = v;
            }
            for &i in &br[k] {
                h[i] = (f[i << 1] * g[i << 1 | 1] - f[i << 1 | 1] * g[i << 1]) * inv2;
                inv2 *= w;
            }
        });
        h
    }

    fn multiply_prefix(f: &mut Self::F, g: &Self::F) {
        pointwise_multiply(f, g);
    }

    fn multiply_add(sum: &mut Self::F, f: &Self::F, g: &Self::F) {
        assert!(sum.len() == f.len() && sum.len() == g.len());
        pointwise_multiply_add(sum, f, g);
    }

    fn power_projection_step(
        p_flat: Vec<MInt<M>>,
        q_flat: Vec<MInt<M>>,
        n: usize,
        py: usize,
        qy: usize,
    ) -> (Vec<MInt<M>>, Vec<MInt<M>>) {
        let high_degree = (qy - 1) * 2;
        let rows = (py + qy - 1).max(high_degree).next_power_of_two();
        let cols = n * 2;
        let size = rows * cols;
        let mut p = p_flat;
        p.resize_with(size, MInt::<M>::zero);
        ntt_rows(&mut p, cols);
        ntt_batch(&mut p, cols);

        let mut q = q_flat;
        q.resize_with(size, MInt::<M>::zero);
        ntt_rows(&mut q, cols);
        let q_high = (rows == high_degree).then(|| q[(qy - 1) * cols..qy * cols].to_vec());
        ntt_batch(&mut q, cols);

        let half = cols / 2;
        let mut odd_factor = vec![MInt::<M>::zero(); half];
        let mut factor = MInt::<M>::from(2).inv();
        let k = cols.trailing_zeros() as usize;
        let w = MInt::<M>::new_unchecked(M::INFO.inv_root[k]);
        BIT_REVERSE.with(|br| {
            let br = unsafe { &mut *br.get() };
            if br.len() < k {
                br.resize_with(k, Default::default);
            }
            let k = k - 1;
            if br[k].is_empty() {
                let mut v = vec![0; 1 << k];
                for i in 0..1 << k {
                    v[i] = (v[i >> 1] >> 1) | ((i & 1) << k.saturating_sub(1));
                }
                br[k] = v;
            }
            for &i in &br[k] {
                odd_factor[i] = factor;
                factor *= w;
            }
        });

        let mut pr = vec![MInt::<M>::zero(); rows * half];
        let mut qr = vec![MInt::<M>::zero(); rows * half];
        for i in 0..pr.len() {
            pr[i] = (p[i << 1] * q[i << 1 | 1] - p[i << 1 | 1] * q[i << 1])
                * odd_factor[i & (half - 1)];
            qr[i] = q[i << 1] * q[i << 1 | 1];
        }
        intt_batch(&mut pr, half);
        intt_rows(&mut pr, half);
        intt_batch(&mut qr, half);
        intt_rows(&mut qr, half);

        if let Some(q_high) = q_high {
            let mut q_high_even = vec![MInt::<M>::zero(); half];
            for i in 0..half {
                q_high_even[i] = q_high[i << 1] * q_high[i << 1 | 1];
            }
            intt(&mut q_high_even);
            for (value, high) in qr.iter_mut().zip(&q_high_even) {
                *value -= *high;
            }
            qr.extend_from_slice(&q_high_even);
        }
        (pr, qr)
    }
}

impl<M, N1, N2, N3> NttReuse for Convolve<(M, (N1, N2, N3))>
where
    M: MIntConvert + MIntConvert<u32>,
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    fn max_product_sum_count(f: &Self::F) -> usize {
        let modulus = <M as MIntConvert<u32>>::mod_into() as u128;
        if modulus == 1 {
            return usize::MAX;
        }
        let capacity = N1::MOD as u128 * N2::MOD as u128 * N3::MOD as u128;
        ((capacity - 1) / ((modulus - 1) * (modulus - 1)) / f.0.len() as u128)
            .clamp(1, usize::MAX as u128) as usize
    }

    fn transform_ntt(t: Self::T, len: usize) -> Self::F {
        let npot = len.max(1).next_power_of_two();
        let f = convert_crt_input(t, npot);
        (
            Convolve::<N1>::transform_ntt(f.0, npot),
            Convolve::<N2>::transform_ntt(f.1, npot),
            Convolve::<N3>::transform_ntt(f.2, npot),
        )
    }

    fn inverse_transform_ntt(f: Self::F, len: usize) -> Self::T {
        reconstruct_mint_crt((
            Convolve::<N1>::inverse_transform_ntt(f.0, len),
            Convolve::<N2>::inverse_transform_ntt(f.1, len),
            Convolve::<N3>::inverse_transform_ntt(f.2, len),
        ))
    }

    fn ntt_doubling(f: Self::F, monic: bool) -> Self::F {
        if monic {
            let n = f.0.len();
            let mut coefficients = Self::inverse_transform_ntt(f, n);
            coefficients[0] -= MInt::<M>::one();
            coefficients.push(MInt::<M>::one());
            Self::transform_ntt(coefficients, n * 2)
        } else {
            (
                Convolve::<N1>::ntt_doubling(f.0, false),
                Convolve::<N2>::ntt_doubling(f.1, false),
                Convolve::<N3>::ntt_doubling(f.2, false),
            )
        }
    }

crates/competitive/src/math/formal_power_series/formal_power_series_impls.rs (line 424)

    fn exp_or_pow(&self, power: Option<T>, deg: usize) -> Self
    where
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        if deg == 1 {
            return Self::one();
        }
        let indices: Vec<_> = (0..=deg).map(T::from).collect();
        let modulus = <T::Base as MIntConvert<usize>>::mod_into();
        let mut inv = vec![T::zero(); deg + 1];
        inv[1] = T::one();
        for i in 2..=deg {
            inv[i] = -T::from(modulus / i) * &inv[modulus % i];
        }
        let block = deg.next_power_of_two() / 16;
        let logarithm = if let Some(rhs) = &power {
            self.prefix_ref(block).log(block) * rhs
        } else {
            self.prefix_ref(block)
        };
        let (kernel, mut kernel_inverse, previous_inverse_fft) =
            logarithm.exp_newton(block, &indices, &inv);
        if power.is_some() {
            kernel_inverse = (self.prefix_ref(block) * &kernel).inv(block);
        } else {
            let mut error_fft = C::transform_ntt(kernel.data.clone(), block);
            C::multiply_prefix(&mut error_fft, &previous_inverse_fft);
            let error = C::inverse_transform_ntt(error_fft, block);
            let mut error_fft =
                C::transform_ntt(error.into_iter().skip(block / 2).collect(), block);
            C::multiply_prefix(&mut error_fft, &previous_inverse_fft);
            let error = C::inverse_transform_ntt(error_fft, block / 2);
            kernel_inverse
                .data
                .extend(error.into_iter().take(block / 2).map(Neg::neg));
        }
        let kernel_data = kernel.data;
        let kernel_inverse_data = kernel_inverse.data;
        let kernel_inverse = C::transform(kernel_inverse_data, block * 2);
        let kernel = C::transform(kernel_data.clone(), block * 2);
        let blocks = deg.div_ceil(block);
        let mut derivative_ffts = Vec::with_capacity(blocks - 1);
        let mut polynomial_ffts = Vec::with_capacity(if power.is_some() { blocks - 1 } else { 0 });
        for q in 1..blocks {
            let mut values = Self::zeros(block * 2);
            for (i, values) in values.data.chunks_mut(block).enumerate() {
                let start = (q - i) * block;
                for (value, x) in values
                    .iter_mut()
                    .zip(self.iter().skip(start).take(deg - start))
                {
                    *value = x.clone();
                }
            }
            if power.is_some() {
                polynomial_ffts.push(C::transform_ntt(values.data.clone(), block * 2));
            }
            for (i, values) in values.data.chunks_mut(block).enumerate() {
                let start = (q - i) * block;
                for (value, index) in values.iter_mut().zip(&indices[start..]) {
                    *value *= index;
                }
            }
            derivative_ffts.push(C::transform_ntt(values.data, block * 2));
        }
        let mut result = kernel_data.clone();
        result.reserve(deg - block);
        let mut result_ffts = Vec::with_capacity(blocks - 1);
        for q in 1..blocks {
            result_ffts.push(C::transform_ntt(
                result[(q - 1) * block..q * block].to_vec(),
                block * 2,
            ));
            let mut values = Self::sum_products(&derivative_ffts[..q], &result_ffts, block);
            if let Some(rhs) = &power {
                let product = Self::sum_products(&polynomial_ffts[..q], &result_ffts, block);
                let factor = rhs.clone() + T::one();
                // The power satisfies f g' = rhs f' g.
                for (i, value) in values.iter_mut().take(deg - q * block).enumerate() {
                    *value = value.clone() * &factor - product[i].clone() * &indices[q * block + i];
                }
            }
            let mut values = C::transform(values, block * 2);
            C::multiply(&mut values, &kernel_inverse);
            let mut values = C::inverse_transform(values, block * 2);
            values.truncate(block);
            let len = block.min(deg - q * block);
            for (i, value) in values.iter_mut().take(len).enumerate() {
                *value *= &inv[q * block + i];
            }
            values[len..].fill(T::zero());
            let mut values = C::transform(values, block * 2);
            C::multiply(&mut values, &kernel);
            let mut values = C::inverse_transform(values, block * 2);
            values.truncate(len);
            result.extend(values);
        }
        Self::from_vec(result)
    }

    fn exp_newton(&self, deg: usize, indices: &[T], inv: &[T]) -> (Self, Self, C::F)
    where
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        if deg == 1 {
            let one = Self::one();
            return (one.clone(), one.clone(), C::transform_ntt(one.data, 1));
        }
        let mut f = Self::from_vec(vec![T::one(), self.coeff(1)]);
        let mut inverse = Self::one();
        let mut inverse_fft = C::transform_ntt(inverse.data.clone(), 2);
        let mut m = 2;
        while m < deg {
            let f_fft = C::transform_ntt(f.data.clone(), 2 * m);

            let previous_inverse_fft = inverse_fft;
            let mut error_fft = previous_inverse_fft.clone();
            C::multiply_prefix(&mut error_fft, &f_fft);
            let mut error = C::inverse_transform_ntt(error_fft, m);
            error[..m / 2].fill(T::zero());
            let mut error_fft = C::transform_ntt(error, m);
            C::multiply_prefix(&mut error_fft, &previous_inverse_fft);
            let error = C::inverse_transform_ntt(error_fft, m);
            inverse
                .data
                .extend(error.into_iter().skip(m / 2).map(Neg::neg));
            inverse_fft = C::transform_ntt(inverse.data.clone(), 2 * m);

            let mut delta = Self::from_vec(
                self.data
                    .iter()
                    .take(m)
                    .enumerate()
                    .skip(1)
                    .map(|(i, value)| value.clone() * &indices[i])
                    .collect(),
            );
            delta.resize(m);
            let mut delta_fft = C::transform_ntt(delta.data, m);
            C::multiply_prefix(&mut delta_fft, &f_fft);
            let mut delta = Self::from_vec(C::inverse_transform_ntt(delta_fft, m));
            for i in 1..f.length() {
                delta[i - 1] -= f[i].clone() * &indices[i];
            }
            delta.resize(2 * m);
            for i in (0..m - 1).rev() {
                delta.data[m + i] = delta.data[i].clone();
            }
            delta.data[..m - 1].fill(T::zero());
            let mut delta_fft = C::transform_ntt(delta.data, 2 * m);
            C::multiply_prefix(&mut delta_fft, &inverse_fft);
            let mut delta = C::inverse_transform_ntt(delta_fft, 2 * m);
            delta.pop();
            delta.push(T::zero());
            let target = (2 * m).min(deg);
            for i in (1..target).rev() {
                delta[i] = delta[i - 1].clone() * &inv[i];
            }
            delta[0] = T::zero();
            delta[target..].fill(T::zero());
            for i in m..(2 * m).min(self.length()) {
                delta[i] += self[i].clone();
            }
            delta[..m].fill(T::zero());
            let mut delta_fft = C::transform_ntt(delta, 2 * m);
            C::multiply_prefix(&mut delta_fft, &f_fft);
            let delta = C::inverse_transform_ntt(delta_fft, 2 * m);
            f.data
                .extend(delta.into_iter().skip(m).take((deg - m).min(m)));
            m *= 2;
        }
        (f, inverse, inverse_fft)
    }
    pub fn log(&self, deg: usize) -> Self {
        if deg == 0 {
            return Self::zero();
        }
        debug_assert!(!self[0].is_zero());
        if deg == 1 {
            return Self::zeros(1);
        }
        if let Some(step) = self.sparse_stride(deg, 2) {
            let pos: Vec<_> = self
                .iter()
                .take(deg)
                .enumerate()
                .skip(1)
                .filter_map(|(i, x)| (!x.is_zero()).then_some(i))
                .collect();
            let mut derivative = Self::zeros(deg);
            let inverse = T::one() / self[0].clone();
            for i in (pos.first().copied().unwrap_or(deg)..deg).step_by(step) {
                let mut value = self.coeff(i) * T::from(i);
                for &j in &pos {
                    if j >= i {
                        break;
                    }
                    value -= self[j].clone() * &derivative[i - j];
                }
                derivative[i] = value * &inverse;
            }
            if pos.is_empty() {
                return derivative;
            }
            derivative.data.remove(0);
            return derivative.integral();
        }
        let n = deg - 1;
        if n <= 64 {
            return (self.inv(deg) * self.prefix_ref(deg).diff())
                .prefix(n)
                .integral();
        }
        let half = n.next_power_of_two() / 2;
        let derivative = self.prefix_ref(deg).diff();
        let inverse = C::transform(self.inv(half).data, half * 2);
        let mut quotient = C::transform(derivative.prefix_ref(half).data, half * 2);
        C::multiply(&mut quotient, &inverse);
        let mut result = C::inverse_transform(quotient, half * 2);
        result.truncate(half);
        if n - half <= 4 {
            let inverse = T::one() / self[0].clone();
            for i in half..n {
                let mut value = derivative.coeff(i);
                for j in 1..=i.min(self.length() - 1) {
                    value -= self[j].clone() * &result[i - j];
                }
                result.push(value * &inverse);
            }
            return Self::from_vec(result).integral();
        }
        let quotient = C::transform(result.clone(), half * 2);
        let mut error = C::transform(self.prefix_ref(n).data, half * 2);
        C::multiply(&mut error, &quotient);
        let mut error = C::inverse_transform(error, half * 2);
        for i in 0..n - half {
            error[i] = derivative.coeff(half + i) - &error[half + i];
        }
        error.truncate(n - half);
        let mut error = C::transform(error, half * 2);
        C::multiply(&mut error, &inverse);
        let error = C::inverse_transform(error, half * 2);
        result.extend(error.into_iter().take(n - half));
        Self::from_vec(result).integral()
    }
    pub fn pow(&self, rhs: usize, deg: usize) -> Self
    where
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        if rhs == 0 {
            return Self::from_vec(
                once(T::one())
                    .chain(repeat_with(T::zero))
                    .take(deg)
                    .collect(),
            );
        }
        if rhs == 1 {
            return self.prefix_ref(deg).resized(deg);
        }
        if let Some(k) = self
            .iter()
            .take(deg.div_ceil(rhs))
            .position(|x| !x.is_zero())
        {
            let deg = deg - k * rhs;
            let x0 = self[k].clone();
            let mut f = (self.prefix_ref(k + deg) >> k) / &x0;
            if let Some(step) = f.sparse_stride(deg, 12) {
                f = f.pow_sparse1(T::from(rhs), deg, step);
            } else if rhs <= 4 {
                let squared = (&f * &f).prefix(deg);
                f = match rhs {
                    2 => squared,
                    3 => (squared * f).prefix(deg),
                    _ => (&squared * &squared).prefix(deg),
                }
                .resized(deg);
            } else {
                f = f.exp_or_pow(Some(T::from(rhs)), deg);
            }
            f *= x0.pow(rhs);
            f <<= k * rhs;
            f
        } else {
            Self::zeros(deg)
        }
    }
    fn pow_sparse1(&self, rhs: T, deg: usize, step: usize) -> Self {
        debug_assert!(!self[0].is_zero());
        let mut pos: Vec<_> = self
            .data
            .iter()
            .take(deg)
            .enumerate()
            .skip(1)
            .filter(|(_, x)| !x.is_zero())
            .map(|(i, x)| (i, x.clone(), T::from(i) * &rhs * x))
            .collect();
        let mut f = Self::zeros(deg);
        f[0] = T::one();
        if pos.is_empty() {
            return f;
        }
        let mf = T::memorized_factorial(deg);
        for (_, coefficient, _) in &mut pos {
            *coefficient *= T::from(step);
        }
        for i in (pos.first().map_or(deg, |x| x.0)..deg).step_by(step) {
            let mut tot = T::zero();
            for (j, coefficient, weight) in &mut pos {
                if *j > i {
                    break;
                }
                tot += weight.clone() * &f[i - *j];
                *weight -= &*coefficient;
            }
            f[i] = tot * T::memorized_inv(&mf, i);
        }
        f
    }

    fn sparse_fold(&self, sparse: impl IntoIterator<Item = (usize, T)>, deg: usize) -> T {
        sparse
            .into_iter()
            .take_while(|&(i, _)| i <= deg)
            .fold(T::zero(), |sum, (i, x)| sum + x * self.coeff(deg - i))
    }

    /// solve: $X(QF)'=\alpha P'(QF)+\beta P(Q'F)$ in $O(deg * max(nz(P), nz(Q), nz(X)))$
    pub fn solve_sparse_differential2(
        p: &Self,
        q: &Self,
        x: &Self,
        alpha: T,
        beta: T,
        deg: usize,
    ) -> Self {
        if deg == 0 {
            return Self::zero();
        }
        let collect_sparse = |p: &Self| -> Vec<(usize, T)> {
            p.iter()
                .enumerate()
                .filter(|&(_, x)| !x.is_zero())
                .map(|(i, x)| (i, x.clone()))
                .collect()
        };
        assert!(q.coeff(0).is_one());
        assert!(x.coeff(0).is_one());
        let p = collect_sparse(p);
        let q = collect_sparse(q);
        let x = collect_sparse(x);
        let diff = |p: &[(usize, T)]| -> Vec<(usize, T)> {
            p.iter()
                .filter(|&&(i, _)| i > 0)
                .map(|&(i, ref x)| (i - 1, x.clone() * T::from(i)))
                .collect()
        };
        let dp = diff(&p);
        let dq = diff(&q);

        let mf = T::memorized_factorial(deg);
        let mut f = Self::zeros(deg);
        let mut qf = Self::zeros(deg);
        let mut dq_f = Self::zeros(deg);
        let mut d_qf = Self::zeros(deg);
        f[0] = T::one();
        for i in 0..deg - 1 {
            qf[i] = f.sparse_fold(q.iter().cloned(), i);
            dq_f[i] = f.sparse_fold(dq.iter().cloned(), i);
            let dp_qf_i = qf.sparse_fold(dp.iter().cloned(), i);
            let p_dq_f_i = dq_f.sparse_fold(p.iter().cloned(), i);
            let x_d_qf_i = d_qf.sparse_fold(
                x.iter()
                    .map(|&(i, ref x)| (i, x.clone() - T::from((i == 0) as usize))),
                i,
            );
            d_qf[i] = alpha.clone() * dp_qf_i + beta.clone() * p_dq_f_i - x_d_qf_i;

            let mut f_ip1 = d_qf[i].clone();
            for &(j, ref q) in q.iter().take_while(|&&(j, _)| j <= i) {
                if j > 0 {
                    f_ip1 -= q.clone() * &f[i - (j - 1)] * T::from(i - (j - 1));
                }
            }
            f[i + 1] = f_ip1 * T::memorized_inv(&mf, i + 1);
        }
        f
    }

    /// P^exp_p * Q^exp_q
    pub fn mul_of_pow_sparse(&self, q: &Self, exp_p: isize, exp_q: isize, deg: usize) -> Self {
        if deg == 0 {
            return Self::zero();
        }
        if exp_p == 0 && exp_q == 0 {
            return Self::from_vec(
                once(T::one())
                    .chain(repeat_with(T::zero))
                    .take(deg)
                    .collect(),
            );
        }
        if exp_p != 0 && self.iter().all(|x| x.is_zero()) {
            assert!(exp_p > 0);
            return Self::zeros(deg);
        }
        if exp_q != 0 && q.iter().all(|x| x.is_zero()) {
            assert!(exp_q > 0);
            return Self::zeros(deg);
        }

        let normalize = |f: &Self, exp: isize| {
            if exp == 0 {
                return (0usize, T::one(), Self::from_vec(vec![T::one()]));
            }
            let k = f.iter().position(|value| !value.is_zero()).unwrap();
            assert!(
                exp >= 0 || k == 0,
                "Negative exponent with zero constant term"
            );
            let c = f[k].clone();
            let f = (f.clone() >> k) / &c;
            (k, c, f)
        };
        let (sp, cp, mut p) = normalize(self, exp_p);
        let (sq, cq, mut q) = normalize(q, exp_q);

        let shift = exp_p
            .saturating_mul(sp as _)
            .saturating_add(exp_q.saturating_mul(sq as _)) as usize;
        if shift >= deg {
            return Self::zeros(deg);
        }
        p.truncate(deg - shift);
        q.truncate(deg - shift);

        let mut f = Self::solve_sparse_differential2(
            &p,
            &q,
            &p,
            T::from(exp_p),
            T::from(exp_q),
            deg - shift,
        );
        f *= cp.signed_pow(exp_p) * cq.signed_pow(exp_q);
        if shift > 0 {
            f <<= shift;
        }
        f.prefix(deg)
    }

    /// exp(P/Q)
    pub fn exp_of_div_sparse(&self, q: &Self, deg: usize) -> Self {
        if deg == 0 {
            return Self::zero();
        }
        let shift_q = q
            .iter()
            .position(|value| !value.is_zero())
            .expect("Zero denominator");
        let shift_p = self.iter().position(|value| !value.is_zero()).unwrap_or(!0);
        assert!(shift_p > shift_q);

        let mut p = self >> shift_q;
        let mut q = q >> shift_q;
        assert!(!q.coeff(0).is_zero());

        let c = q[0].clone();
        p /= c.clone();
        q /= c;

        Self::solve_sparse_differential2(&p, &q, &q, T::one(), -T::one(), deg)
    }
}

impl<T, C> FormalPowerSeries<T, C>
where
    T: FormalPowerSeriesCoefficientSqrt,
    C: ConvolveSteps<T = Vec<T>>,
{
    pub fn sqrt(&self, deg: usize) -> Option<Self> {
        if self[0].is_zero() {
            if let Some(k) = self.iter().position(|x| !x.is_zero()) {
                if k % 2 != 0 {
                    return None;
                } else if deg > k / 2 {
                    return Some((self >> k).sqrt(deg - k / 2)? << (k / 2));
                }
            }
        } else {
            let s = self[0].sqrt_coefficient()?;
            if deg <= 1 {
                return Some(Self::from(s).prefix(deg));
            }
            if let Some(step) = self.sparse_stride(deg, 4) {
                let t = self[0].clone();
                let mut f = self.prefix_ref(deg) / t;
                f = f.pow_sparse1(T::one() / T::from(2usize), deg, step);
                f *= s;
                return Some(f);
            }

            let mut f = Self::from(s);
            let inv2 = T::one() / (T::one() + T::one());
            let inv2s = inv2.clone() / &f[0];
            let extend = |f: &mut Self, end| {
                for i in f.length()..end {
                    let mut value = self.coeff(i);
                    for j in 1..i {
                        value -= f[j].clone() * &f[i - j];
                    }
                    f.data.push(value * &inv2s);
                }
            };
            extend(&mut f, deg.min(32));
            f.truncate(deg);
            if f.length() == deg {
                return Some(f);
            }
            let mut inverse = f.inv(f.length());
            let mut i = f.length();
            while i < deg {
                if deg - i <= 4 {
                    extend(&mut f, deg);
                    break;
                }
                let len = (i * 2).min(deg);
                let factor = C::transform(inverse.data.clone(), i * 2);
                let error = if !C::CYCLIC || i < 128 {
                    (self.prefix_ref(len) - &f * &f) >> i
                } else {
                    let square = C::square(f.data.clone(), i);
                    // The cyclic square folds its high half into the already known low half.
                    Self::from_vec(
                        square
                            .into_iter()
                            .take(len - i)
                            .enumerate()
                            .map(|(j, value)| self.coeff(i + j) + self.coeff(j) - value)
                            .collect(),
                    )
                };
                let mut error_fft = C::transform(error.data, i * 2);
                C::multiply(&mut error_fft, &factor);
                let delta = C::inverse_transform(error_fft, i * 2);
                f.data
                    .extend(delta.into_iter().take(len - i).map(|x| x * &inv2));
                if i * 2 + 4 < deg {
                    let mut error_fft = C::transform(f.data.clone(), i * 2);
                    C::multiply(&mut error_fft, &factor);
                    let error = C::inverse_transform(error_fft, i * 2);
                    let mut error_fft = C::transform(error.into_iter().skip(i).collect(), i * 2);
                    C::multiply(&mut error_fft, &factor);
                    let error = C::inverse_transform(error_fft, i * 2);
                    inverse.data.extend(error.into_iter().take(i).map(Neg::neg));
                }
                i *= 2;
            }
            f.truncate(deg);
            return Some(f);
        }
        Some(Self::zeros(deg))
    }
}

impl<T, C> FormalPowerSeries<T, C>
where
    T: FormalPowerSeriesCoefficient,
    C: ConvolveSteps<T = Vec<T>>,
{
    pub fn count_subset_sum<F>(&self, deg: usize, mut inverse: F) -> Self
    where
        F: FnMut(usize) -> T,
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        let n = self.length();
        let mut f = Self::zeros(n);
        for i in 1..n {
            if !self[i].is_zero() {
                for (j, d) in (0..n).step_by(i).enumerate().skip(1) {
                    if j & 1 != 0 {
                        f[d] += self[i].clone() * &inverse(j);
                    } else {
                        f[d] -= self[i].clone() * &inverse(j);
                    }
                }
            }
        }
        f.exp(deg)
    }
    pub fn count_multiset_sum<F>(&self, deg: usize, mut inverse: F) -> Self
    where
        F: FnMut(usize) -> T,
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        let n = self.length();
        let mut f = Self::zeros(n);
        for i in 1..n {
            if !self[i].is_zero() {
                for (j, d) in (0..n).step_by(i).enumerate().skip(1) {
                    f[d] += self[i].clone() * &inverse(j);
                }
            }
        }
        f.exp(deg)
    }
    /// [x^n] P(x) / Q(x)
    pub fn bostan_mori(mut self, mut rhs: Self, mut n: usize) -> T
    where
        C: NttReuse<T = Vec<T>>,
    {
        let mut res = T::zero();
        rhs.trim_tail_zeros();
        if self.length() >= rhs.length() {
            let r = &self / &rhs;
            if n < r.length() {
                res = r[n].clone();
            }
            self -= r * &rhs;
            self.trim_tail_zeros();
        }
        let mut k = rhs.length().next_power_of_two();
        let mut p = C::transform_ntt(self.data, k * 2);
        let mut q = C::transform_ntt(rhs.data, k * 2);
        while n > 0 {
            let t = C::even_mul_normal_neg(&q, &q);
            p = if n.is_multiple_of(2) {
                C::even_mul_normal_neg(&p, &q)
            } else {
                C::odd_mul_normal_neg(&p, &q)
            };
            q = t;
            n /= 2;
            if n != 0 {
                if n < k / 2 {
                    p = C::transform_ntt(C::inverse_transform_ntt(p, k / 2), k);
                    q = C::transform_ntt(C::inverse_transform_ntt(q, k / 2), k);
                    k /= 2;
                } else if C::MULTIPLE {
                    p = C::transform_ntt(C::inverse_transform_ntt(p, k), k * 2);
                    q = C::transform_ntt(C::inverse_transform_ntt(q, k), k * 2);
                } else {
                    p = C::ntt_doubling(p, false);
                    q = C::ntt_doubling(q, false);
                }
            }
        }
        let p = C::inverse_transform_ntt(p, k);
        let q = C::inverse_transform_ntt(q, k);
        res + p[0].clone() / q[0].clone()
    }
    /// return F(x) where [x^n] P(x) / Q(x) = [x^d-1] P(x) F(x)
    pub fn bostan_mori_msb(self, n: usize) -> Self {
        let d = self.length() - 1;
        if n == 0 {
            return (Self::one() << (d - 1)) / self[0].clone();
        }
        let q = self;
        let mq = q.clone().parity_inversion();
        let w = (q * &mq).even().bostan_mori_msb(n / 2);
        let mut s = Self::zeros(w.length() * 2 - (n % 2));
        for (i, x) in w.iter().enumerate() {
            s[i * 2 + (1 - n % 2)] = x.clone();
        }
        let len = 2 * d + 1;
        let ts = C::transform(s.prefix(len).data, len);
        mq.reversed().middle_product(&ts, len).prefix(d + 1)
    }
    /// x^n mod self
    pub fn pow_mod(self, n: usize) -> Self {
        let d = self.length() - 1;
        let q = self.reversed();
        let u = q.clone().bostan_mori_msb(n);
        let mut f = (u * q).prefix(d).reversed();
        f.trim_tail_zeros();
        f
    }
    fn middle_product(self, other: &C::F, deg: usize) -> Self {
        let n = self.length();
        let mut s = C::transform(self.reversed().data, deg);
        C::multiply(&mut s, other);
        Self::from_vec((C::inverse_transform(s, deg))[n - 1..].to_vec())
    }
    pub fn multipoint_evaluation(self, points: &[T]) -> Vec<T>
    where
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        let n = points.len();
        if n <= 32 || self.length() <= 32 {
            return points.iter().map(|p| self.eval(p.clone())).collect();
        }
        let size = n.next_power_of_two();
        let block = 16;
        let leaves = size / block;
        let mut subproduct_tree = Vec::with_capacity(leaves * 2);
        subproduct_tree.resize_with(leaves * 2, || None);
        let mut leaf_products = Vec::with_capacity(leaves);
        for i in 0..leaves {
            let mut product = vec![T::one()];
            for j in 0..block {
                let x = points.get(i * block + j).cloned().unwrap_or_else(T::zero);
                product.push(T::one());
                for k in (1..=j).rev() {
                    product[k] = product[k - 1].clone() - x.clone() * &product[k];
                }
                product[0] *= -x;
            }
            subproduct_tree[leaves + i] = Some(C::transform_ntt(product.clone(), block * 2));
            leaf_products.push(product);
        }
        for i in (1..leaves).rev() {
            let mut product = subproduct_tree[i * 2].as_ref().unwrap().clone();
            C::multiply_prefix(&mut product, subproduct_tree[i * 2 + 1].as_ref().unwrap());
            if i > 1 {
                product = C::ntt_doubling(product, true);
            }
            subproduct_tree[i] = Some(product);
        }
        let mut product = C::inverse_transform_ntt(subproduct_tree[1].take().unwrap(), size);
        product[0] -= T::one();
        product.push(T::one());
        let mut uptree_t = Vec::with_capacity(leaves * 2);
        uptree_t.resize_with(1, Zero::zero);
        let m = self.length();
        let v = Self::from_vec(product).reversed().resized(m);
        let s = C::transform(self.data, m * 2);
        uptree_t.push(v.inv(m).middle_product(&s, m * 2).resized(size));
        for i in 1..leaves {
            let degree = uptree_t[i].length();
            let spectrum = C::transform_ntt(std::mem::take(&mut uptree_t[i].data), degree);
            let left = subproduct_tree[i * 2].take().unwrap();
            let right = subproduct_tree[i * 2 + 1].take().unwrap();
            let mut child = spectrum.clone();
            C::multiply_prefix(&mut child, &right);
            let mut child = C::inverse_transform_ntt(child, degree);
            child.drain(..degree / 2);
            uptree_t.push(Self::from_vec(child));
            let mut child = spectrum;
            C::multiply_prefix(&mut child, &left);
            let mut child = C::inverse_transform_ntt(child, degree);
            child.drain(..degree / 2);
            uptree_t.push(Self::from_vec(child));
        }
        let mut result = Vec::with_capacity(n);
        for ((values, product), points) in uptree_t[leaves..]
            .iter()
            .zip(leaf_products)
            .zip(points.chunks(block))
        {
            let mut remainder = Self::zeros(block);
            for (j, value) in values.iter().enumerate() {
                for (r, p) in remainder.data[..=j].iter_mut().zip(&product[block - j..]) {
                    *r += value.clone() * p;
                }
            }
            result.extend(points.iter().map(|p| remainder.eval(p.clone())));
        }
        result
    }

fn inverse_transform_ntt(f: Self::F, len: usize) -> Self::T

Inverts a value produced by transform_ntt.

Examples found in repository

crates/competitive/src/math/formal_power_series/berlekamp_massey.rs (line 133)

    fn apply(&self, p: &C::F, q: &C::F, length: usize) -> (Vec<T>, Vec<T>) {
        (
            C::inverse_transform_ntt(Self::product_sum(p, &self.a00, q, &self.a01), length),
            C::inverse_transform_ntt(Self::product_sum(p, &self.a10, q, &self.a11), length),
        )
    }

    fn left_multiply_step(self, quotient: &FormalPowerSeries<T, C>, length: usize) -> Self {
        let negative_quotient = reduced_transform(&(-quotient), length);
        let mut a10 = self.a00;
        C::multiply_add(&mut a10, &negative_quotient, &self.a10);
        let mut a11 = self.a01;
        C::multiply_add(&mut a11, &negative_quotient, &self.a11);
        let result = Self {
            a00: self.a10,
            a01: self.a11,
            a10,
            a11,
        };
        if C::MULTIPLE {
            result.inverse_transform(length).transform(length)
        } else {
            result
        }
    }

    fn inverse_transform(self, length: usize) -> FpsMatrix<T, C> {
        FpsMatrix {
            a00: FormalPowerSeries::from_vec(C::inverse_transform_ntt(self.a00, length)),
            a01: FormalPowerSeries::from_vec(C::inverse_transform_ntt(self.a01, length)),
            a10: FormalPowerSeries::from_vec(C::inverse_transform_ntt(self.a10, length)),
            a11: FormalPowerSeries::from_vec(C::inverse_transform_ntt(self.a11, length)),
        }
    }

crates/competitive/src/math/formal_power_series/formal_power_series_impls.rs (line 387)

    fn sum_products(f: &[C::F], g: &[C::F], len: usize) -> Vec<T>
    where
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        let chunk = C::max_product_sum_count(&f[0]);
        f.rchunks(chunk)
            .zip(g.chunks(chunk))
            .map(|(f, g)| {
                let mut sum = f[f.len() - 1].clone();
                C::multiply_prefix(&mut sum, &g[0]);
                for (f, g) in f.iter().rev().skip(1).zip(&g[1..]) {
                    C::multiply_add(&mut sum, f, g);
                }
                C::inverse_transform_ntt(sum, len)
            })
            .reduce(|mut sum, part| {
                for (sum, value) in sum.iter_mut().zip(part) {
                    *sum += value;
                }
                sum
            })
            .unwrap()
    }

    fn exp_or_pow(&self, power: Option<T>, deg: usize) -> Self
    where
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        if deg == 1 {
            return Self::one();
        }
        let indices: Vec<_> = (0..=deg).map(T::from).collect();
        let modulus = <T::Base as MIntConvert<usize>>::mod_into();
        let mut inv = vec![T::zero(); deg + 1];
        inv[1] = T::one();
        for i in 2..=deg {
            inv[i] = -T::from(modulus / i) * &inv[modulus % i];
        }
        let block = deg.next_power_of_two() / 16;
        let logarithm = if let Some(rhs) = &power {
            self.prefix_ref(block).log(block) * rhs
        } else {
            self.prefix_ref(block)
        };
        let (kernel, mut kernel_inverse, previous_inverse_fft) =
            logarithm.exp_newton(block, &indices, &inv);
        if power.is_some() {
            kernel_inverse = (self.prefix_ref(block) * &kernel).inv(block);
        } else {
            let mut error_fft = C::transform_ntt(kernel.data.clone(), block);
            C::multiply_prefix(&mut error_fft, &previous_inverse_fft);
            let error = C::inverse_transform_ntt(error_fft, block);
            let mut error_fft =
                C::transform_ntt(error.into_iter().skip(block / 2).collect(), block);
            C::multiply_prefix(&mut error_fft, &previous_inverse_fft);
            let error = C::inverse_transform_ntt(error_fft, block / 2);
            kernel_inverse
                .data
                .extend(error.into_iter().take(block / 2).map(Neg::neg));
        }
        let kernel_data = kernel.data;
        let kernel_inverse_data = kernel_inverse.data;
        let kernel_inverse = C::transform(kernel_inverse_data, block * 2);
        let kernel = C::transform(kernel_data.clone(), block * 2);
        let blocks = deg.div_ceil(block);
        let mut derivative_ffts = Vec::with_capacity(blocks - 1);
        let mut polynomial_ffts = Vec::with_capacity(if power.is_some() { blocks - 1 } else { 0 });
        for q in 1..blocks {
            let mut values = Self::zeros(block * 2);
            for (i, values) in values.data.chunks_mut(block).enumerate() {
                let start = (q - i) * block;
                for (value, x) in values
                    .iter_mut()
                    .zip(self.iter().skip(start).take(deg - start))
                {
                    *value = x.clone();
                }
            }
            if power.is_some() {
                polynomial_ffts.push(C::transform_ntt(values.data.clone(), block * 2));
            }
            for (i, values) in values.data.chunks_mut(block).enumerate() {
                let start = (q - i) * block;
                for (value, index) in values.iter_mut().zip(&indices[start..]) {
                    *value *= index;
                }
            }
            derivative_ffts.push(C::transform_ntt(values.data, block * 2));
        }
        let mut result = kernel_data.clone();
        result.reserve(deg - block);
        let mut result_ffts = Vec::with_capacity(blocks - 1);
        for q in 1..blocks {
            result_ffts.push(C::transform_ntt(
                result[(q - 1) * block..q * block].to_vec(),
                block * 2,
            ));
            let mut values = Self::sum_products(&derivative_ffts[..q], &result_ffts, block);
            if let Some(rhs) = &power {
                let product = Self::sum_products(&polynomial_ffts[..q], &result_ffts, block);
                let factor = rhs.clone() + T::one();
                // The power satisfies f g' = rhs f' g.
                for (i, value) in values.iter_mut().take(deg - q * block).enumerate() {
                    *value = value.clone() * &factor - product[i].clone() * &indices[q * block + i];
                }
            }
            let mut values = C::transform(values, block * 2);
            C::multiply(&mut values, &kernel_inverse);
            let mut values = C::inverse_transform(values, block * 2);
            values.truncate(block);
            let len = block.min(deg - q * block);
            for (i, value) in values.iter_mut().take(len).enumerate() {
                *value *= &inv[q * block + i];
            }
            values[len..].fill(T::zero());
            let mut values = C::transform(values, block * 2);
            C::multiply(&mut values, &kernel);
            let mut values = C::inverse_transform(values, block * 2);
            values.truncate(len);
            result.extend(values);
        }
        Self::from_vec(result)
    }

    fn exp_newton(&self, deg: usize, indices: &[T], inv: &[T]) -> (Self, Self, C::F)
    where
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        if deg == 1 {
            let one = Self::one();
            return (one.clone(), one.clone(), C::transform_ntt(one.data, 1));
        }
        let mut f = Self::from_vec(vec![T::one(), self.coeff(1)]);
        let mut inverse = Self::one();
        let mut inverse_fft = C::transform_ntt(inverse.data.clone(), 2);
        let mut m = 2;
        while m < deg {
            let f_fft = C::transform_ntt(f.data.clone(), 2 * m);

            let previous_inverse_fft = inverse_fft;
            let mut error_fft = previous_inverse_fft.clone();
            C::multiply_prefix(&mut error_fft, &f_fft);
            let mut error = C::inverse_transform_ntt(error_fft, m);
            error[..m / 2].fill(T::zero());
            let mut error_fft = C::transform_ntt(error, m);
            C::multiply_prefix(&mut error_fft, &previous_inverse_fft);
            let error = C::inverse_transform_ntt(error_fft, m);
            inverse
                .data
                .extend(error.into_iter().skip(m / 2).map(Neg::neg));
            inverse_fft = C::transform_ntt(inverse.data.clone(), 2 * m);

            let mut delta = Self::from_vec(
                self.data
                    .iter()
                    .take(m)
                    .enumerate()
                    .skip(1)
                    .map(|(i, value)| value.clone() * &indices[i])
                    .collect(),
            );
            delta.resize(m);
            let mut delta_fft = C::transform_ntt(delta.data, m);
            C::multiply_prefix(&mut delta_fft, &f_fft);
            let mut delta = Self::from_vec(C::inverse_transform_ntt(delta_fft, m));
            for i in 1..f.length() {
                delta[i - 1] -= f[i].clone() * &indices[i];
            }
            delta.resize(2 * m);
            for i in (0..m - 1).rev() {
                delta.data[m + i] = delta.data[i].clone();
            }
            delta.data[..m - 1].fill(T::zero());
            let mut delta_fft = C::transform_ntt(delta.data, 2 * m);
            C::multiply_prefix(&mut delta_fft, &inverse_fft);
            let mut delta = C::inverse_transform_ntt(delta_fft, 2 * m);
            delta.pop();
            delta.push(T::zero());
            let target = (2 * m).min(deg);
            for i in (1..target).rev() {
                delta[i] = delta[i - 1].clone() * &inv[i];
            }
            delta[0] = T::zero();
            delta[target..].fill(T::zero());
            for i in m..(2 * m).min(self.length()) {
                delta[i] += self[i].clone();
            }
            delta[..m].fill(T::zero());
            let mut delta_fft = C::transform_ntt(delta, 2 * m);
            C::multiply_prefix(&mut delta_fft, &f_fft);
            let delta = C::inverse_transform_ntt(delta_fft, 2 * m);
            f.data
                .extend(delta.into_iter().skip(m).take((deg - m).min(m)));
            m *= 2;
        }
        (f, inverse, inverse_fft)
    }
    pub fn log(&self, deg: usize) -> Self {
        if deg == 0 {
            return Self::zero();
        }
        debug_assert!(!self[0].is_zero());
        if deg == 1 {
            return Self::zeros(1);
        }
        if let Some(step) = self.sparse_stride(deg, 2) {
            let pos: Vec<_> = self
                .iter()
                .take(deg)
                .enumerate()
                .skip(1)
                .filter_map(|(i, x)| (!x.is_zero()).then_some(i))
                .collect();
            let mut derivative = Self::zeros(deg);
            let inverse = T::one() / self[0].clone();
            for i in (pos.first().copied().unwrap_or(deg)..deg).step_by(step) {
                let mut value = self.coeff(i) * T::from(i);
                for &j in &pos {
                    if j >= i {
                        break;
                    }
                    value -= self[j].clone() * &derivative[i - j];
                }
                derivative[i] = value * &inverse;
            }
            if pos.is_empty() {
                return derivative;
            }
            derivative.data.remove(0);
            return derivative.integral();
        }
        let n = deg - 1;
        if n <= 64 {
            return (self.inv(deg) * self.prefix_ref(deg).diff())
                .prefix(n)
                .integral();
        }
        let half = n.next_power_of_two() / 2;
        let derivative = self.prefix_ref(deg).diff();
        let inverse = C::transform(self.inv(half).data, half * 2);
        let mut quotient = C::transform(derivative.prefix_ref(half).data, half * 2);
        C::multiply(&mut quotient, &inverse);
        let mut result = C::inverse_transform(quotient, half * 2);
        result.truncate(half);
        if n - half <= 4 {
            let inverse = T::one() / self[0].clone();
            for i in half..n {
                let mut value = derivative.coeff(i);
                for j in 1..=i.min(self.length() - 1) {
                    value -= self[j].clone() * &result[i - j];
                }
                result.push(value * &inverse);
            }
            return Self::from_vec(result).integral();
        }
        let quotient = C::transform(result.clone(), half * 2);
        let mut error = C::transform(self.prefix_ref(n).data, half * 2);
        C::multiply(&mut error, &quotient);
        let mut error = C::inverse_transform(error, half * 2);
        for i in 0..n - half {
            error[i] = derivative.coeff(half + i) - &error[half + i];
        }
        error.truncate(n - half);
        let mut error = C::transform(error, half * 2);
        C::multiply(&mut error, &inverse);
        let error = C::inverse_transform(error, half * 2);
        result.extend(error.into_iter().take(n - half));
        Self::from_vec(result).integral()
    }
    pub fn pow(&self, rhs: usize, deg: usize) -> Self
    where
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        if rhs == 0 {
            return Self::from_vec(
                once(T::one())
                    .chain(repeat_with(T::zero))
                    .take(deg)
                    .collect(),
            );
        }
        if rhs == 1 {
            return self.prefix_ref(deg).resized(deg);
        }
        if let Some(k) = self
            .iter()
            .take(deg.div_ceil(rhs))
            .position(|x| !x.is_zero())
        {
            let deg = deg - k * rhs;
            let x0 = self[k].clone();
            let mut f = (self.prefix_ref(k + deg) >> k) / &x0;
            if let Some(step) = f.sparse_stride(deg, 12) {
                f = f.pow_sparse1(T::from(rhs), deg, step);
            } else if rhs <= 4 {
                let squared = (&f * &f).prefix(deg);
                f = match rhs {
                    2 => squared,
                    3 => (squared * f).prefix(deg),
                    _ => (&squared * &squared).prefix(deg),
                }
                .resized(deg);
            } else {
                f = f.exp_or_pow(Some(T::from(rhs)), deg);
            }
            f *= x0.pow(rhs);
            f <<= k * rhs;
            f
        } else {
            Self::zeros(deg)
        }
    }
    fn pow_sparse1(&self, rhs: T, deg: usize, step: usize) -> Self {
        debug_assert!(!self[0].is_zero());
        let mut pos: Vec<_> = self
            .data
            .iter()
            .take(deg)
            .enumerate()
            .skip(1)
            .filter(|(_, x)| !x.is_zero())
            .map(|(i, x)| (i, x.clone(), T::from(i) * &rhs * x))
            .collect();
        let mut f = Self::zeros(deg);
        f[0] = T::one();
        if pos.is_empty() {
            return f;
        }
        let mf = T::memorized_factorial(deg);
        for (_, coefficient, _) in &mut pos {
            *coefficient *= T::from(step);
        }
        for i in (pos.first().map_or(deg, |x| x.0)..deg).step_by(step) {
            let mut tot = T::zero();
            for (j, coefficient, weight) in &mut pos {
                if *j > i {
                    break;
                }
                tot += weight.clone() * &f[i - *j];
                *weight -= &*coefficient;
            }
            f[i] = tot * T::memorized_inv(&mf, i);
        }
        f
    }

    fn sparse_fold(&self, sparse: impl IntoIterator<Item = (usize, T)>, deg: usize) -> T {
        sparse
            .into_iter()
            .take_while(|&(i, _)| i <= deg)
            .fold(T::zero(), |sum, (i, x)| sum + x * self.coeff(deg - i))
    }

    /// solve: $X(QF)'=\alpha P'(QF)+\beta P(Q'F)$ in $O(deg * max(nz(P), nz(Q), nz(X)))$
    pub fn solve_sparse_differential2(
        p: &Self,
        q: &Self,
        x: &Self,
        alpha: T,
        beta: T,
        deg: usize,
    ) -> Self {
        if deg == 0 {
            return Self::zero();
        }
        let collect_sparse = |p: &Self| -> Vec<(usize, T)> {
            p.iter()
                .enumerate()
                .filter(|&(_, x)| !x.is_zero())
                .map(|(i, x)| (i, x.clone()))
                .collect()
        };
        assert!(q.coeff(0).is_one());
        assert!(x.coeff(0).is_one());
        let p = collect_sparse(p);
        let q = collect_sparse(q);
        let x = collect_sparse(x);
        let diff = |p: &[(usize, T)]| -> Vec<(usize, T)> {
            p.iter()
                .filter(|&&(i, _)| i > 0)
                .map(|&(i, ref x)| (i - 1, x.clone() * T::from(i)))
                .collect()
        };
        let dp = diff(&p);
        let dq = diff(&q);

        let mf = T::memorized_factorial(deg);
        let mut f = Self::zeros(deg);
        let mut qf = Self::zeros(deg);
        let mut dq_f = Self::zeros(deg);
        let mut d_qf = Self::zeros(deg);
        f[0] = T::one();
        for i in 0..deg - 1 {
            qf[i] = f.sparse_fold(q.iter().cloned(), i);
            dq_f[i] = f.sparse_fold(dq.iter().cloned(), i);
            let dp_qf_i = qf.sparse_fold(dp.iter().cloned(), i);
            let p_dq_f_i = dq_f.sparse_fold(p.iter().cloned(), i);
            let x_d_qf_i = d_qf.sparse_fold(
                x.iter()
                    .map(|&(i, ref x)| (i, x.clone() - T::from((i == 0) as usize))),
                i,
            );
            d_qf[i] = alpha.clone() * dp_qf_i + beta.clone() * p_dq_f_i - x_d_qf_i;

            let mut f_ip1 = d_qf[i].clone();
            for &(j, ref q) in q.iter().take_while(|&&(j, _)| j <= i) {
                if j > 0 {
                    f_ip1 -= q.clone() * &f[i - (j - 1)] * T::from(i - (j - 1));
                }
            }
            f[i + 1] = f_ip1 * T::memorized_inv(&mf, i + 1);
        }
        f
    }

    /// P^exp_p * Q^exp_q
    pub fn mul_of_pow_sparse(&self, q: &Self, exp_p: isize, exp_q: isize, deg: usize) -> Self {
        if deg == 0 {
            return Self::zero();
        }
        if exp_p == 0 && exp_q == 0 {
            return Self::from_vec(
                once(T::one())
                    .chain(repeat_with(T::zero))
                    .take(deg)
                    .collect(),
            );
        }
        if exp_p != 0 && self.iter().all(|x| x.is_zero()) {
            assert!(exp_p > 0);
            return Self::zeros(deg);
        }
        if exp_q != 0 && q.iter().all(|x| x.is_zero()) {
            assert!(exp_q > 0);
            return Self::zeros(deg);
        }

        let normalize = |f: &Self, exp: isize| {
            if exp == 0 {
                return (0usize, T::one(), Self::from_vec(vec![T::one()]));
            }
            let k = f.iter().position(|value| !value.is_zero()).unwrap();
            assert!(
                exp >= 0 || k == 0,
                "Negative exponent with zero constant term"
            );
            let c = f[k].clone();
            let f = (f.clone() >> k) / &c;
            (k, c, f)
        };
        let (sp, cp, mut p) = normalize(self, exp_p);
        let (sq, cq, mut q) = normalize(q, exp_q);

        let shift = exp_p
            .saturating_mul(sp as _)
            .saturating_add(exp_q.saturating_mul(sq as _)) as usize;
        if shift >= deg {
            return Self::zeros(deg);
        }
        p.truncate(deg - shift);
        q.truncate(deg - shift);

        let mut f = Self::solve_sparse_differential2(
            &p,
            &q,
            &p,
            T::from(exp_p),
            T::from(exp_q),
            deg - shift,
        );
        f *= cp.signed_pow(exp_p) * cq.signed_pow(exp_q);
        if shift > 0 {
            f <<= shift;
        }
        f.prefix(deg)
    }

    /// exp(P/Q)
    pub fn exp_of_div_sparse(&self, q: &Self, deg: usize) -> Self {
        if deg == 0 {
            return Self::zero();
        }
        let shift_q = q
            .iter()
            .position(|value| !value.is_zero())
            .expect("Zero denominator");
        let shift_p = self.iter().position(|value| !value.is_zero()).unwrap_or(!0);
        assert!(shift_p > shift_q);

        let mut p = self >> shift_q;
        let mut q = q >> shift_q;
        assert!(!q.coeff(0).is_zero());

        let c = q[0].clone();
        p /= c.clone();
        q /= c;

        Self::solve_sparse_differential2(&p, &q, &q, T::one(), -T::one(), deg)
    }
}

impl<T, C> FormalPowerSeries<T, C>
where
    T: FormalPowerSeriesCoefficientSqrt,
    C: ConvolveSteps<T = Vec<T>>,
{
    pub fn sqrt(&self, deg: usize) -> Option<Self> {
        if self[0].is_zero() {
            if let Some(k) = self.iter().position(|x| !x.is_zero()) {
                if k % 2 != 0 {
                    return None;
                } else if deg > k / 2 {
                    return Some((self >> k).sqrt(deg - k / 2)? << (k / 2));
                }
            }
        } else {
            let s = self[0].sqrt_coefficient()?;
            if deg <= 1 {
                return Some(Self::from(s).prefix(deg));
            }
            if let Some(step) = self.sparse_stride(deg, 4) {
                let t = self[0].clone();
                let mut f = self.prefix_ref(deg) / t;
                f = f.pow_sparse1(T::one() / T::from(2usize), deg, step);
                f *= s;
                return Some(f);
            }

            let mut f = Self::from(s);
            let inv2 = T::one() / (T::one() + T::one());
            let inv2s = inv2.clone() / &f[0];
            let extend = |f: &mut Self, end| {
                for i in f.length()..end {
                    let mut value = self.coeff(i);
                    for j in 1..i {
                        value -= f[j].clone() * &f[i - j];
                    }
                    f.data.push(value * &inv2s);
                }
            };
            extend(&mut f, deg.min(32));
            f.truncate(deg);
            if f.length() == deg {
                return Some(f);
            }
            let mut inverse = f.inv(f.length());
            let mut i = f.length();
            while i < deg {
                if deg - i <= 4 {
                    extend(&mut f, deg);
                    break;
                }
                let len = (i * 2).min(deg);
                let factor = C::transform(inverse.data.clone(), i * 2);
                let error = if !C::CYCLIC || i < 128 {
                    (self.prefix_ref(len) - &f * &f) >> i
                } else {
                    let square = C::square(f.data.clone(), i);
                    // The cyclic square folds its high half into the already known low half.
                    Self::from_vec(
                        square
                            .into_iter()
                            .take(len - i)
                            .enumerate()
                            .map(|(j, value)| self.coeff(i + j) + self.coeff(j) - value)
                            .collect(),
                    )
                };
                let mut error_fft = C::transform(error.data, i * 2);
                C::multiply(&mut error_fft, &factor);
                let delta = C::inverse_transform(error_fft, i * 2);
                f.data
                    .extend(delta.into_iter().take(len - i).map(|x| x * &inv2));
                if i * 2 + 4 < deg {
                    let mut error_fft = C::transform(f.data.clone(), i * 2);
                    C::multiply(&mut error_fft, &factor);
                    let error = C::inverse_transform(error_fft, i * 2);
                    let mut error_fft = C::transform(error.into_iter().skip(i).collect(), i * 2);
                    C::multiply(&mut error_fft, &factor);
                    let error = C::inverse_transform(error_fft, i * 2);
                    inverse.data.extend(error.into_iter().take(i).map(Neg::neg));
                }
                i *= 2;
            }
            f.truncate(deg);
            return Some(f);
        }
        Some(Self::zeros(deg))
    }
}

impl<T, C> FormalPowerSeries<T, C>
where
    T: FormalPowerSeriesCoefficient,
    C: ConvolveSteps<T = Vec<T>>,
{
    pub fn count_subset_sum<F>(&self, deg: usize, mut inverse: F) -> Self
    where
        F: FnMut(usize) -> T,
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        let n = self.length();
        let mut f = Self::zeros(n);
        for i in 1..n {
            if !self[i].is_zero() {
                for (j, d) in (0..n).step_by(i).enumerate().skip(1) {
                    if j & 1 != 0 {
                        f[d] += self[i].clone() * &inverse(j);
                    } else {
                        f[d] -= self[i].clone() * &inverse(j);
                    }
                }
            }
        }
        f.exp(deg)
    }
    pub fn count_multiset_sum<F>(&self, deg: usize, mut inverse: F) -> Self
    where
        F: FnMut(usize) -> T,
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        let n = self.length();
        let mut f = Self::zeros(n);
        for i in 1..n {
            if !self[i].is_zero() {
                for (j, d) in (0..n).step_by(i).enumerate().skip(1) {
                    f[d] += self[i].clone() * &inverse(j);
                }
            }
        }
        f.exp(deg)
    }
    /// [x^n] P(x) / Q(x)
    pub fn bostan_mori(mut self, mut rhs: Self, mut n: usize) -> T
    where
        C: NttReuse<T = Vec<T>>,
    {
        let mut res = T::zero();
        rhs.trim_tail_zeros();
        if self.length() >= rhs.length() {
            let r = &self / &rhs;
            if n < r.length() {
                res = r[n].clone();
            }
            self -= r * &rhs;
            self.trim_tail_zeros();
        }
        let mut k = rhs.length().next_power_of_two();
        let mut p = C::transform_ntt(self.data, k * 2);
        let mut q = C::transform_ntt(rhs.data, k * 2);
        while n > 0 {
            let t = C::even_mul_normal_neg(&q, &q);
            p = if n.is_multiple_of(2) {
                C::even_mul_normal_neg(&p, &q)
            } else {
                C::odd_mul_normal_neg(&p, &q)
            };
            q = t;
            n /= 2;
            if n != 0 {
                if n < k / 2 {
                    p = C::transform_ntt(C::inverse_transform_ntt(p, k / 2), k);
                    q = C::transform_ntt(C::inverse_transform_ntt(q, k / 2), k);
                    k /= 2;
                } else if C::MULTIPLE {
                    p = C::transform_ntt(C::inverse_transform_ntt(p, k), k * 2);
                    q = C::transform_ntt(C::inverse_transform_ntt(q, k), k * 2);
                } else {
                    p = C::ntt_doubling(p, false);
                    q = C::ntt_doubling(q, false);
                }
            }
        }
        let p = C::inverse_transform_ntt(p, k);
        let q = C::inverse_transform_ntt(q, k);
        res + p[0].clone() / q[0].clone()
    }
    /// return F(x) where [x^n] P(x) / Q(x) = [x^d-1] P(x) F(x)
    pub fn bostan_mori_msb(self, n: usize) -> Self {
        let d = self.length() - 1;
        if n == 0 {
            return (Self::one() << (d - 1)) / self[0].clone();
        }
        let q = self;
        let mq = q.clone().parity_inversion();
        let w = (q * &mq).even().bostan_mori_msb(n / 2);
        let mut s = Self::zeros(w.length() * 2 - (n % 2));
        for (i, x) in w.iter().enumerate() {
            s[i * 2 + (1 - n % 2)] = x.clone();
        }
        let len = 2 * d + 1;
        let ts = C::transform(s.prefix(len).data, len);
        mq.reversed().middle_product(&ts, len).prefix(d + 1)
    }
    /// x^n mod self
    pub fn pow_mod(self, n: usize) -> Self {
        let d = self.length() - 1;
        let q = self.reversed();
        let u = q.clone().bostan_mori_msb(n);
        let mut f = (u * q).prefix(d).reversed();
        f.trim_tail_zeros();
        f
    }
    fn middle_product(self, other: &C::F, deg: usize) -> Self {
        let n = self.length();
        let mut s = C::transform(self.reversed().data, deg);
        C::multiply(&mut s, other);
        Self::from_vec((C::inverse_transform(s, deg))[n - 1..].to_vec())
    }
    pub fn multipoint_evaluation(self, points: &[T]) -> Vec<T>
    where
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        let n = points.len();
        if n <= 32 || self.length() <= 32 {
            return points.iter().map(|p| self.eval(p.clone())).collect();
        }
        let size = n.next_power_of_two();
        let block = 16;
        let leaves = size / block;
        let mut subproduct_tree = Vec::with_capacity(leaves * 2);
        subproduct_tree.resize_with(leaves * 2, || None);
        let mut leaf_products = Vec::with_capacity(leaves);
        for i in 0..leaves {
            let mut product = vec![T::one()];
            for j in 0..block {
                let x = points.get(i * block + j).cloned().unwrap_or_else(T::zero);
                product.push(T::one());
                for k in (1..=j).rev() {
                    product[k] = product[k - 1].clone() - x.clone() * &product[k];
                }
                product[0] *= -x;
            }
            subproduct_tree[leaves + i] = Some(C::transform_ntt(product.clone(), block * 2));
            leaf_products.push(product);
        }
        for i in (1..leaves).rev() {
            let mut product = subproduct_tree[i * 2].as_ref().unwrap().clone();
            C::multiply_prefix(&mut product, subproduct_tree[i * 2 + 1].as_ref().unwrap());
            if i > 1 {
                product = C::ntt_doubling(product, true);
            }
            subproduct_tree[i] = Some(product);
        }
        let mut product = C::inverse_transform_ntt(subproduct_tree[1].take().unwrap(), size);
        product[0] -= T::one();
        product.push(T::one());
        let mut uptree_t = Vec::with_capacity(leaves * 2);
        uptree_t.resize_with(1, Zero::zero);
        let m = self.length();
        let v = Self::from_vec(product).reversed().resized(m);
        let s = C::transform(self.data, m * 2);
        uptree_t.push(v.inv(m).middle_product(&s, m * 2).resized(size));
        for i in 1..leaves {
            let degree = uptree_t[i].length();
            let spectrum = C::transform_ntt(std::mem::take(&mut uptree_t[i].data), degree);
            let left = subproduct_tree[i * 2].take().unwrap();
            let right = subproduct_tree[i * 2 + 1].take().unwrap();
            let mut child = spectrum.clone();
            C::multiply_prefix(&mut child, &right);
            let mut child = C::inverse_transform_ntt(child, degree);
            child.drain(..degree / 2);
            uptree_t.push(Self::from_vec(child));
            let mut child = spectrum;
            C::multiply_prefix(&mut child, &left);
            let mut child = C::inverse_transform_ntt(child, degree);
            child.drain(..degree / 2);
            uptree_t.push(Self::from_vec(child));
        }
        let mut result = Vec::with_capacity(n);
        for ((values, product), points) in uptree_t[leaves..]
            .iter()
            .zip(leaf_products)
            .zip(points.chunks(block))
        {
            let mut remainder = Self::zeros(block);
            for (j, value) in values.iter().enumerate() {
                for (r, p) in remainder.data[..=j].iter_mut().zip(&product[block - j..]) {
                    *r += value.clone() * p;
                }
            }
            result.extend(points.iter().map(|p| remainder.eval(p.clone())));
        }
        result
    }

crates/competitive/src/math/number_theoretic_transform.rs (line 611)

fn convolve_large_ntt<M>(a: Vec<MInt<M>>, b: Vec<MInt<M>>) -> Vec<MInt<M>>
where
    M: Montgomery32NttModulus,
{
    let len = a.len() + b.len() - 1;
    let ntt_len = 1usize << M::RANK;
    let block_len = ntt_len / 2;
    let same = a == b;
    let transform = |a: &[MInt<M>]| {
        let mut f = Vec::with_capacity(ntt_len);
        advise_huge_pages(&mut f);
        f.extend_from_slice(a);
        Convolve::<M>::transform_ntt(f, ntt_len)
    };
    let fa: Vec<_> = a.chunks(block_len).map(transform).collect();
    let fb: Option<Vec<_>> = if same {
        None
    } else {
        Some(b.chunks(block_len).map(transform).collect())
    };
    let b_blocks = fb.as_ref().map_or(fa.len(), Vec::len);
    let mut result = vec![MInt::<M>::zero(); len];
    for diagonal in 0..fa.len() + b_blocks - 1 {
        let mut spectrum = vec![MInt::<M>::zero(); ntt_len];
        let start = diagonal.saturating_sub(b_blocks - 1);
        for i in start..=diagonal.min(fa.len() - 1) {
            let j = diagonal - i;
            let g = if let Some(fb) = &fb { &fb[j] } else { &fa[j] };
            pointwise_multiply_add(&mut spectrum, &fa[i], g);
        }
        spectrum = Convolve::<M>::inverse_transform_ntt(spectrum, ntt_len);
        let offset = diagonal * block_len;
        for (result, value) in result[offset..].iter_mut().zip(spectrum) {
            *result += value;
        }
    }
    result
}

impl<M> ConvolveSteps for Convolve<M>
where
    M: Montgomery32NttModulus,
{
    const CYCLIC: bool = true;

    type T = Vec<MInt<M>>;
    type F = Vec<MInt<M>>;
    fn length(t: &Self::T) -> usize {
        t.len()
    }
    fn transform(mut t: Self::T, len: usize) -> Self::F {
        t.resize_with(len.max(1).next_power_of_two(), Zero::zero);
        #[cfg(target_arch = "x86_64")]
        if use_block_ntt::<M>(t.len()) {
            unsafe { ntt_simd::transform_blocks_avx2(&mut t) };
            return t;
        }
        ntt(&mut t);
        t
    }
    fn inverse_transform(mut f: Self::F, len: usize) -> Self::T {
        #[cfg(target_arch = "x86_64")]
        if use_block_ntt::<M>(f.len()) {
            unsafe { ntt_simd::inverse_transform_blocks_avx2(&mut f) };
            f.truncate(len);
            return f;
        }
        intt(&mut f);
        f.truncate(len);
        f
    }
    fn multiply(f: &mut Self::F, g: &Self::F) {
        assert_eq!(f.len(), g.len());
        #[cfg(target_arch = "x86_64")]
        if use_block_ntt::<M>(f.len()) {
            unsafe { ntt_simd::multiply_blocks_avx2(f, g) };
            return;
        }
        pointwise_multiply(f, g);
    }
    fn square(t: Self::T, len: usize) -> Self::T {
        let mut f = Self::transform(t, len);
        let g = f.clone();
        Self::multiply(&mut f, &g);
        Self::inverse_transform(f, len)
    }
    fn convolve(mut a: Self::T, mut b: Self::T) -> Self::T {
        let (threshold, naive_threshold) = (100, 60);
        #[cfg(target_arch = "x86_64")]
        let (threshold, naive_threshold) = if use_block_ntt::<M>(64) {
            (
                60,
                if M::RANK >= 13 && a.len().max(b.len()) <= 4096 {
                    18
                } else if M::RANK >= 19 && a.len().max(b.len()) <= 262144 {
                    32
                } else {
                    34
                },
            )
        } else {
            (threshold, naive_threshold)
        };
        if Self::length(&a).max(Self::length(&b)) <= threshold {
            return convolve_karatsuba(&a, &b);
        }
        if Self::length(&a).min(Self::length(&b)) <= naive_threshold {
            return convolve_naive(&a, &b);
        }
        let len = (Self::length(&a) + Self::length(&b)).saturating_sub(1);
        let size = len.max(1).next_power_of_two();
        let max_size = 1usize << M::RANK;
        #[cfg(target_arch = "x86_64")]
        let max_size = if use_block_ntt::<M>(size) {
            max_size << 3
        } else {
            max_size
        };
        if size > max_size {
            return convolve_large_ntt(a, b);
        }
        if len <= size / 2 + 2 {
            let xa = a.pop().unwrap();
            let xb = b.pop().unwrap();
            let mut c = vec![MInt::<M>::zero(); len];
            *c.last_mut().unwrap() = xa * xb;
            for (a, c) in a.iter().zip(&mut c[b.len()..]) {
                *c += *a * xb;
            }
            for (b, c) in b.iter().zip(&mut c[a.len()..]) {
                *c += *b * xa;
            }
            let d = Self::convolve(a, b);
            for (d, c) in d.into_iter().zip(&mut c) {
                *c += d;
            }
            return c;
        }
        let same = a == b;
        #[cfg(target_arch = "x86_64")]
        if use_block_ntt::<M>(size) {
            a.reserve(size - a.len());
            b.reserve(size - b.len());
            advise_huge_pages(&mut a);
            advise_huge_pages(&mut b);
            a.resize_with(size, Zero::zero);
            b.resize_with(size, Zero::zero);
            unsafe { ntt_simd::convolve_blocks_avx2(&mut a, &mut b, same) };
            a.truncate(len);
            return a;
        }
        let mut a = Self::transform(a, len);
        if same {
            for a in a.iter_mut() {
                *a *= *a;
            }
        } else {
            let b = Self::transform(b, len);
            Self::multiply(&mut a, &b);
        }
        Self::inverse_transform(a, len)
    }
}

type MVec<M> = Vec<MInt<M>>;

fn convert_crt_input<M, N1, N2, N3>(t: MVec<M>, capacity: usize) -> (MVec<N1>, MVec<N2>, MVec<N3>)
where
    M: MIntConvert<u32>,
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    let mut f = (
        MVec::<N1>::with_capacity(capacity),
        MVec::<N2>::with_capacity(capacity),
        MVec::<N3>::with_capacity(capacity),
    );
    advise_huge_pages(&mut f.0);
    advise_huge_pages(&mut f.1);
    advise_huge_pages(&mut f.2);
    for t in t {
        let t: u32 = t.into();
        f.0.push(t.into());
        f.1.push(t.into());
        f.2.push(t.into());
    }
    f
}

fn reconstruct_mint_crt<M, N1, N2, N3>(f: (MVec<N1>, MVec<N2>, MVec<N3>)) -> MVec<M>
where
    M: MIntConvert + MIntConvert<u32>,
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    let t1 = MInt::<N2>::new(N1::get_mod()).inv();
    let m1_3 = MInt::<N3>::new(N1::get_mod());
    let t2 = (m1_3 * MInt::<N3>::new(N2::get_mod())).inv();
    let modulus = <M as MIntConvert<u32>>::mod_into() as u64;
    let m1 = N1::get_mod() as u64;
    let m2 = m1 * N2::get_mod() as u64 % modulus;
    let fits_u64 = (N1::get_mod() - 1) as u128
        + (N2::get_mod() - 1) as u128 * m1 as u128
        + (N3::get_mod() - 1) as u128 * m2 as u128
        <= u64::MAX as u128;
    f.0.into_iter()
        .zip(f.1)
        .zip(f.2)
        .map(|((c1, c2), c3)| {
            let d1 = c1.inner();
            let d2 = ((c2 - MInt::<N2>::from(d1)) * t1).inner();
            let x = MInt::<N3>::new(d1) + MInt::<N3>::new(d2) * m1_3;
            let d3 = ((c3 - x) * t2).inner();
            let value = if fits_u64 {
                (d1 as u64 + d2 as u64 * m1 + d3 as u64 * m2) % modulus
            } else {
                ((d1 as u128 + d2 as u128 * m1 as u128 + d3 as u128 * m2 as u128) % modulus as u128)
                    as u64
            };
            MInt::<M>::from(value as u32)
        })
        .collect()
}

impl<M, N1, N2, N3> ConvolveSteps for Convolve<(M, (N1, N2, N3))>
where
    M: MIntConvert + MIntConvert<u32>,
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    type T = MVec<M>;
    type F = (MVec<N1>, MVec<N2>, MVec<N3>);
    fn length(t: &Self::T) -> usize {
        t.len()
    }
    fn transform(t: Self::T, len: usize) -> Self::F {
        let npot = len.max(1).next_power_of_two();
        let f = convert_crt_input(t, npot);
        (
            Convolve::<N1>::transform(f.0, npot),
            Convolve::<N2>::transform(f.1, npot),
            Convolve::<N3>::transform(f.2, npot),
        )
    }
    fn inverse_transform(f: Self::F, len: usize) -> Self::T {
        reconstruct_mint_crt((
            Convolve::<N1>::inverse_transform(f.0, len),
            Convolve::<N2>::inverse_transform(f.1, len),
            Convolve::<N3>::inverse_transform(f.2, len),
        ))
    }
    fn multiply(f: &mut Self::F, g: &Self::F) {
        Convolve::<N1>::multiply(&mut f.0, &g.0);
        Convolve::<N2>::multiply(&mut f.1, &g.1);
        Convolve::<N3>::multiply(&mut f.2, &g.2);
    }
    fn convolve(a: Self::T, b: Self::T) -> Self::T {
        let max_len = Self::length(&a).max(Self::length(&b));
        let min_len = Self::length(&a).min(Self::length(&b));
        let (balanced, short) = crate::avx_helper!(@dispatch_avx2_fma (30, 10), (384, 128));
        if max_len <= balanced || min_len <= short {
            return convolve_karatsuba(&a, &b);
        }
        // Limit coefficient growth to leave headroom for FFT roundoff.
        let fft_limit = crate::avx_helper!(@dispatch_avx2_fma
            1usize << ((1u64 << 50) / <M as MIntConvert<u32>>::mod_into() as u64).ilog2().min(20), 0);
        let convolve = |a: Self::T, b: Self::T| {
            let fft_len = (a.len() + b.len() - 1).next_power_of_two();
            if fft_len <= 256 && a.len() * b.len() <= fft_len * 8 {
                return convolve_karatsuba(&a, &b);
            }
            if fft_len <= fft_limit {
                crate::avx_helper!(@dispatch_avx2_fma return unsafe {
                    convolve_mint_avx2(a, b)
                }, ());
            }
            convolve_mint_crt::<M, N1, N2, N3>(a, b)
        };
        let block_len = min_len.next_power_of_two() * 8 - min_len + 1;
        let block_len = if min_len <= fft_limit / 2 {
            block_len.min(fft_limit - min_len + 1)
        } else {
            block_len
        };
        if max_len <= block_len {
            return convolve(a, b);
        }
        let (a, b) = if a.len() >= b.len() { (a, b) } else { (b, a) };
        let mut result = vec![MInt::<M>::zero(); a.len() + b.len() - 1];
        for (i, a) in a.chunks(block_len).enumerate() {
            let product = convolve(a.to_vec(), b.clone());
            for (value, product) in result[i * block_len..].iter_mut().zip(product) {
                *value += product;
            }
        }
        result
    }
}

fn convolve_mint_crt<M, N1, N2, N3>(a: MVec<M>, b: MVec<M>) -> MVec<M>
where
    M: MIntConvert + MIntConvert<u32>,
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    let convolve = |a: MVec<M>, b: MVec<M>| {
        let a_len = a.len();
        let b_len = b.len();
        let a = convert_crt_input(a, a_len);
        let b = convert_crt_input(b, b_len);
        reconstruct_mint_crt((
            Convolve::<N1>::convolve(a.0, b.0),
            Convolve::<N2>::convolve(a.1, b.1),
            Convolve::<N3>::convolve(a.2, b.2),
        ))
    };
    let modulus = <M as MIntConvert<u32>>::mod_into() as u128;
    let capacity = N1::MOD as u128 * N2::MOD as u128 * N3::MOD as u128;
    if a.len().min(b.len()) as u128 * (modulus - 1).pow(2) < capacity {
        return convolve(a, b);
    }
    let block_len = ((capacity - 1) / (modulus - 1).pow(2)) as usize;
    if block_len == 0 {
        return convolve_naive(&a, &b);
    }
    let mut result = vec![MInt::<M>::zero(); a.len() + b.len() - 1];
    for (i, a) in a.chunks(block_len).enumerate() {
        for (j, b) in b.chunks(block_len).enumerate() {
            let product = convolve(a.to_vec(), b.to_vec());
            for (value, product) in result[(i + j) * block_len..].iter_mut().zip(product) {
                *value += product;
            }
        }
    }
    result
}

impl<N1, N2, N3> ConvolveSteps for Convolve<(u64, (N1, N2, N3))>
where
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    type T = Vec<u64>;
    type F = ([MVec<N1>; 3], [MVec<N2>; 3], [MVec<N3>; 3]);

    fn length(t: &Self::T) -> usize {
        t.len()
    }

    fn transform(t: Self::T, len: usize) -> Self::F {
        let npot = len.max(1).next_power_of_two();
        assert!(npot <= 1usize << N1::RANK.min(N2::RANK).min(N3::RANK));
        // The 22-bit fallback needs room for three limb products per coefficient.
        assert!(
            3 * npot as u128 * ((1u128 << 22) - 1).pow(2)
                < N1::MOD as u128 * N2::MOD as u128 * N3::MOD as u128
        );
        let bits = if 2 * npot as u128 * (u32::MAX as u128).pow(2)
            < N1::MOD as u128 * N2::MOD as u128 * N3::MOD as u128
        {
            32
        } else {
            22
        };
        let parts = if bits == 32 && t.iter().all(|&value| value <= u32::MAX as u64) {
            1
        } else {
            64usize.div_ceil(bits)
        };
        fn split<M: Montgomery32NttModulus>(
            t: &[u64],
            len: usize,
            bits: usize,
            parts: usize,
        ) -> [MVec<M>; 3] {
            std::array::from_fn(|part| {
                if part >= parts {
                    return Vec::new();
                }
                Convolve::<M>::transform(
                    t.iter()
                        .map(|&t| MInt::from((t >> (part * bits)) & ((1u64 << bits) - 1)))
                        .collect(),
                    len,
                )
            })
        }
        (
            split(&t, npot, bits, parts),
            split(&t, npot, bits, parts),
            split(&t, npot, bits, parts),
        )
    }

    fn inverse_transform(f: Self::F, len: usize) -> Self::T {
        let bits = if f.0[2].is_empty() { 32 } else { 22 };
        let t1 = MInt::<N2>::new(N1::get_mod()).inv();
        let m1 = N1::get_mod() as u64;
        let m1_3 = MInt::<N3>::new(N1::get_mod());
        let t2 = (m1_3 * MInt::<N3>::new(N2::get_mod())).inv();
        let m2 = m1 * N2::get_mod() as u64;
        let mut result = vec![0u64; len.min(f.0[0].len())];
        for (part, ((f1, f2), f3)) in f.0.into_iter().zip(f.1).zip(f.2).enumerate() {
            if f1.is_empty() {
                continue;
            }
            for (value, ((c1, c2), c3)) in result.iter_mut().zip(
                Convolve::<N1>::inverse_transform(f1, len)
                    .into_iter()
                    .zip(Convolve::<N2>::inverse_transform(f2, len))
                    .zip(Convolve::<N3>::inverse_transform(f3, len)),
            ) {
                let d1 = c1.inner();
                let d2 = ((c2 - MInt::<N2>::from(d1)) * t1).inner();
                let x = MInt::<N3>::new(d1) + MInt::<N3>::new(d2) * m1_3;
                let d3 = ((c3 - x) * t2).inner();
                let limb = (d1 as u64)
                    .wrapping_add((d2 as u64).wrapping_mul(m1))
                    .wrapping_add((d3 as u64).wrapping_mul(m2));
                *value = value.wrapping_add(limb << (part * bits));
            }
        }
        result
    }

    fn multiply(f: &mut Self::F, g: &Self::F) {
        fn multiply<M: Montgomery32NttModulus>(f: &mut [MVec<M>; 3], g: &[MVec<M>; 3]) {
            assert_eq!(f[0].len(), g[0].len());
            if f[1].is_empty() || g[1].is_empty() {
                if f[1].is_empty() && !g[1].is_empty() {
                    f[1] = f[0].clone();
                    Convolve::<M>::multiply(&mut f[1], &g[1]);
                } else if !f[1].is_empty() {
                    Convolve::<M>::multiply(&mut f[1], &g[0]);
                }
                Convolve::<M>::multiply(&mut f[0], &g[0]);
                return;
            }
            #[cfg(target_arch = "x86_64")]
            if use_block_ntt::<M>(f[0].len()) {
                for part in (1..if f[2].is_empty() { 2 } else { 3 }).rev() {
                    let mut sum = f[0].clone();
                    Convolve::<M>::multiply(&mut sum, &g[part]);
                    for left in 1..=part {
                        let mut product = f[left].clone();
                        Convolve::<M>::multiply(&mut product, &g[part - left]);
                        for (value, product) in sum.iter_mut().zip(product) {
                            // Block products contain lazy Montgomery residues.
                            *value = MInt::new(value.inner() + product.inner());
                        }
                    }
                    f[part] = sum;
                }
                Convolve::<M>::multiply(&mut f[0], &g[0]);
                return;
            }
            if f[2].is_empty() {
                for i in 0..f[0].len() {
                    f[1][i] = f[0][i] * g[1][i] + f[1][i] * g[0][i];
                    f[0][i] *= g[0][i];
                }
                return;
            }
            for i in 0..f[0].len() {
                f[2][i] = f[0][i] * g[2][i] + f[1][i] * g[1][i] + f[2][i] * g[0][i];
                f[1][i] = f[0][i] * g[1][i] + f[1][i] * g[0][i];
                f[0][i] *= g[0][i];
            }
        }
        multiply(&mut f.0, &g.0);
        multiply(&mut f.1, &g.1);
        multiply(&mut f.2, &g.2);
    }

    fn square(t: Self::T, len: usize) -> Self::T {
        let mut f = Self::transform(t, len);
        let g = f.clone();
        Self::multiply(&mut f, &g);
        Self::inverse_transform(f, len)
    }

    fn convolve(a: Self::T, b: Self::T) -> Self::T {
        let max_len = Self::length(&a).max(Self::length(&b));
        let min_len = Self::length(&a).min(Self::length(&b));
        let (balanced, short) = crate::avx_helper!(@dispatch_avx2_fma (300, 64), (1536, 512));
        if max_len <= balanced || min_len <= short {
            let a_wrapping: &[Wrapping<u64>] =
                unsafe { std::slice::from_raw_parts(a.as_ptr().cast(), a.len()) };
            let b_wrapping: &[Wrapping<u64>] =
                unsafe { std::slice::from_raw_parts(b.as_ptr().cast(), b.len()) };
            let mut c = std::mem::ManuallyDrop::new(if max_len <= 300 || min_len > 60 {
                convolve_karatsuba(a_wrapping, b_wrapping)
            } else {
                convolve_naive(a_wrapping, b_wrapping)
            });
            return unsafe { Vec::from_raw_parts(c.as_mut_ptr().cast(), c.len(), c.capacity()) };
        }
        let len = (Self::length(&a) + Self::length(&b)).saturating_sub(1);
        let block_len = if min_len >= 1 << 20 {
            1 << 20
        } else {
            (min_len.next_power_of_two() * 8).min(1 << 21) - min_len + 1
        };
        if max_len <= block_len {
            return convolve_u64_fft(a, b);
        }
        let mut result = vec![0u64; len];
        for (i, a) in a.chunks(block_len).enumerate() {
            for (j, b) in b.chunks(block_len).enumerate() {
                if a.len().min(b.len()) <= 60 {
                    for (x, &a) in a.iter().enumerate() {
                        for (y, &b) in b.iter().enumerate() {
                            let value = &mut result[(i + j) * block_len + x + y];
                            *value = value.wrapping_add(a.wrapping_mul(b));
                        }
                    }
                    continue;
                }
                let product = convolve_u64_fft(a.to_vec(), b.to_vec());
                for (value, product) in result[(i + j) * block_len..].iter_mut().zip(product) {
                    *value = value.wrapping_add(product);
                }
            }
        }
        result
    }
}

fn convolve_u64_fft(a: Vec<u64>, b: Vec<u64>) -> Vec<u64> {
    // Keep limb convolutions below 2^47 at the 2^21 FFT limit.
    crate::avx_helper!(@dispatch_avx2_fma return unsafe {
        convolve_u64_avx2(a, b)
    }, ());
    convolve_u64_fft_scalar(a, b)
}

fn convolve_u64_fft_scalar(a: Vec<u64>, b: Vec<u64>) -> Vec<u64> {
    fn split(values: &[u64]) -> [Vec<i64>; 5] {
        let mut result = std::array::from_fn(|_| Vec::with_capacity(values.len()));
        for mut value in values.iter().copied() {
            for part in &mut result {
                let digit = ((value << 51) as i64) >> 51;
                part.push(digit);
                value = (value >> 13).wrapping_add(u64::from(digit < 0));
            }
        }
        result
    }

    let len = a.len() + b.len() - 1;
    let transform = |values: &[u64]| {
        if values.iter().any(|&value| value > u32::MAX as u64) {
            return split(values).map(|part| ConvolveRealFft::transform(part, len));
        }
        let [a, b, c, _, _] = split(values);
        let a = ConvolveRealFft::transform(a, len);
        let size = a.len();
        [
            a,
            ConvolveRealFft::transform(b, len),
            ConvolveRealFft::transform(c, len),
            vec![Zero::zero(); size],
            vec![Zero::zero(); size],
        ]
    };
    let fa = transform(&a);
    drop(a);
    let fb = transform(&b);
    drop(b);
    let values: [Vec<i64>; 5] = std::array::from_fn(|part| {
        let mut sum = fa[0].clone();
        ConvolveRealFft::multiply(&mut sum, &fb[part]);
        for left in 1..=part {
            let mut product = fa[left].clone();
            ConvolveRealFft::multiply(&mut product, &fb[part - left]);
            for (sum, product) in sum.iter_mut().zip(product) {
                *sum += product;
            }
        }
        ConvolveRealFft::inverse_transform(sum, len)
    });
    (0..len)
        .map(|i| {
            (values[0][i] as u64)
                .wrapping_add((values[1][i] as u64) << 13)
                .wrapping_add((values[2][i] as u64) << 26)
                .wrapping_add((values[3][i] as u64) << 39)
                .wrapping_add((values[4][i] as u64) << 52)
        })
        .collect()
}

pub trait NttReuse: ConvolveSteps {
    const MULTIPLE: bool = true;

    /// Transforms coefficients into the usual NTT frequency order.
    fn transform_ntt(t: Self::T, len: usize) -> Self::F {
        Self::transform(t, len)
    }

    /// Inverts a value produced by `transform_ntt`.
    fn inverse_transform_ntt(f: Self::F, len: usize) -> Self::T {
        Self::inverse_transform(f, len)
    }

    /// Extends a value produced by `transform_ntt` to twice its length.
    /// If `monic`, the input represents a monic degree-`n` polynomial modulo
    /// `x^n - 1`, where `n` is the transform length.
    fn ntt_doubling(f: Self::F, monic: bool) -> Self::F;

    /// Extracts the even coefficients of `a(x) * b(-x)` in the usual NTT frequency order.
    fn even_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F;

    /// Extracts the odd coefficients of `a(x) * b(-x)` in the usual NTT frequency order.
    fn odd_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F;

    /// Multiplies a usual NTT transform by the corresponding prefix of another one.
    fn multiply_prefix(f: &mut Self::F, g: &Self::F);

    /// Adds the pointwise product of two usual NTT transforms to `sum`.
    fn multiply_add(sum: &mut Self::F, f: &Self::F, g: &Self::F);

    /// Maximum number of products that can be summed before reconstruction.
    /// Both factors must transform canonical coefficients at the supplied transform's length,
    /// and each cyclic product must itself be reconstructible.
    fn max_product_sum_count(_f: &Self::F) -> usize {
        if Self::MULTIPLE { 1 } else { usize::MAX }
    }

    fn power_projection_step(
        p_flat: Self::T,
        q_flat: Self::T,
        n: usize,
        py: usize,
        qy: usize,
    ) -> (Self::T, Self::T) {
        let base = n * 2;
        let len_p = base * py;
        let len_q = base * qy;
        let len = (len_p + len_q - 1).max(len_q + len_q - 1);
        let half = len.max(1).next_power_of_two() / 2;

        let p_fft = Self::transform_ntt(p_flat, len);
        let q_fft = Self::transform_ntt(q_flat, len);
        let pr_fft = Self::odd_mul_normal_neg(&p_fft, &q_fft);
        let qr_fft = Self::even_mul_normal_neg(&q_fft, &q_fft);
        (
            Self::inverse_transform_ntt(pr_fft, half),
            Self::inverse_transform_ntt(qr_fft, half),
        )
    }
}

thread_local!(
    static BIT_REVERSE: UnsafeCell<Vec<Vec<usize>>> = const { UnsafeCell::new(vec![]) };
);

impl<M> NttReuse for Convolve<M>
where
    M: Montgomery32NttModulus,
{
    const MULTIPLE: bool = false;

    fn transform_ntt(mut t: Self::T, len: usize) -> Self::F {
        t.resize_with(len.max(1).next_power_of_two(), Zero::zero);
        ntt(&mut t);
        t
    }

    fn inverse_transform_ntt(mut f: Self::F, len: usize) -> Self::T {
        intt(&mut f);
        f.truncate(len);
        f
    }

    fn ntt_doubling(mut f: Self::F, monic: bool) -> Self::F {
        let n = f.len();
        let k = n.trailing_zeros() as usize;
        let mut a = Self::inverse_transform_ntt(f.clone(), n);
        if monic {
            a[0] -= MInt::<M>::from(2);
        }
        let zeta = MInt::<M>::new_unchecked(M::INFO.root[k + 1]);
        let zeta2 = zeta * zeta;
        let mut rot = [MInt::one(), zeta, zeta2, zeta2 * zeta];
        let step = zeta2 * zeta2;
        for a in a.chunks_mut(4) {
            for (a, rot) in a.iter_mut().zip(&mut rot) {
                *a *= *rot;
                *rot *= step;
            }
        }
        f.extend(Self::transform_ntt(a, n));
        f
    }

    fn even_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F {
        assert_eq!(f.len(), g.len());
        assert!(f.len().is_power_of_two());
        assert!(f.len() >= 2);
        if std::ptr::eq(f, g) {
            return f.as_chunks::<2>().0.iter().map(|a| a[0] * a[1]).collect();
        }
        let inv2 = MInt::<M>::from(2).inv();
        let n = f.len() / 2;
        (0..n)
            .map(|i| (f[i << 1] * g[i << 1 | 1] + f[i << 1 | 1] * g[i << 1]) * inv2)
            .collect()
    }

    fn odd_mul_normal_neg(f: &Self::F, g: &Self::F) -> Self::F {
        assert_eq!(f.len(), g.len());
        assert!(f.len().is_power_of_two());
        assert!(f.len() >= 2);
        let mut inv2 = MInt::<M>::from(2).inv();
        let n = f.len() / 2;
        let k = f.len().trailing_zeros() as usize;
        let mut h = vec![MInt::<M>::zero(); n];
        let w = MInt::<M>::new_unchecked(M::INFO.inv_root[k]);
        BIT_REVERSE.with(|br| {
            let br = unsafe { &mut *br.get() };
            if br.len() < k {
                br.resize_with(k, Default::default);
            }
            let k = k - 1;
            if br[k].is_empty() {
                let mut v = vec![0; 1 << k];
                for i in 0..1 << k {
                    v[i] = (v[i >> 1] >> 1) | ((i & 1) << k.saturating_sub(1));
                }
                br[k] = v;
            }
            for &i in &br[k] {
                h[i] = (f[i << 1] * g[i << 1 | 1] - f[i << 1 | 1] * g[i << 1]) * inv2;
                inv2 *= w;
            }
        });
        h
    }

    fn multiply_prefix(f: &mut Self::F, g: &Self::F) {
        pointwise_multiply(f, g);
    }

    fn multiply_add(sum: &mut Self::F, f: &Self::F, g: &Self::F) {
        assert!(sum.len() == f.len() && sum.len() == g.len());
        pointwise_multiply_add(sum, f, g);
    }

    fn power_projection_step(
        p_flat: Vec<MInt<M>>,
        q_flat: Vec<MInt<M>>,
        n: usize,
        py: usize,
        qy: usize,
    ) -> (Vec<MInt<M>>, Vec<MInt<M>>) {
        let high_degree = (qy - 1) * 2;
        let rows = (py + qy - 1).max(high_degree).next_power_of_two();
        let cols = n * 2;
        let size = rows * cols;
        let mut p = p_flat;
        p.resize_with(size, MInt::<M>::zero);
        ntt_rows(&mut p, cols);
        ntt_batch(&mut p, cols);

        let mut q = q_flat;
        q.resize_with(size, MInt::<M>::zero);
        ntt_rows(&mut q, cols);
        let q_high = (rows == high_degree).then(|| q[(qy - 1) * cols..qy * cols].to_vec());
        ntt_batch(&mut q, cols);

        let half = cols / 2;
        let mut odd_factor = vec![MInt::<M>::zero(); half];
        let mut factor = MInt::<M>::from(2).inv();
        let k = cols.trailing_zeros() as usize;
        let w = MInt::<M>::new_unchecked(M::INFO.inv_root[k]);
        BIT_REVERSE.with(|br| {
            let br = unsafe { &mut *br.get() };
            if br.len() < k {
                br.resize_with(k, Default::default);
            }
            let k = k - 1;
            if br[k].is_empty() {
                let mut v = vec![0; 1 << k];
                for i in 0..1 << k {
                    v[i] = (v[i >> 1] >> 1) | ((i & 1) << k.saturating_sub(1));
                }
                br[k] = v;
            }
            for &i in &br[k] {
                odd_factor[i] = factor;
                factor *= w;
            }
        });

        let mut pr = vec![MInt::<M>::zero(); rows * half];
        let mut qr = vec![MInt::<M>::zero(); rows * half];
        for i in 0..pr.len() {
            pr[i] = (p[i << 1] * q[i << 1 | 1] - p[i << 1 | 1] * q[i << 1])
                * odd_factor[i & (half - 1)];
            qr[i] = q[i << 1] * q[i << 1 | 1];
        }
        intt_batch(&mut pr, half);
        intt_rows(&mut pr, half);
        intt_batch(&mut qr, half);
        intt_rows(&mut qr, half);

        if let Some(q_high) = q_high {
            let mut q_high_even = vec![MInt::<M>::zero(); half];
            for i in 0..half {
                q_high_even[i] = q_high[i << 1] * q_high[i << 1 | 1];
            }
            intt(&mut q_high_even);
            for (value, high) in qr.iter_mut().zip(&q_high_even) {
                *value -= *high;
            }
            qr.extend_from_slice(&q_high_even);
        }
        (pr, qr)
    }
}

impl<M, N1, N2, N3> NttReuse for Convolve<(M, (N1, N2, N3))>
where
    M: MIntConvert + MIntConvert<u32>,
    N1: Montgomery32NttModulus,
    N2: Montgomery32NttModulus,
    N3: Montgomery32NttModulus,
{
    fn max_product_sum_count(f: &Self::F) -> usize {
        let modulus = <M as MIntConvert<u32>>::mod_into() as u128;
        if modulus == 1 {
            return usize::MAX;
        }
        let capacity = N1::MOD as u128 * N2::MOD as u128 * N3::MOD as u128;
        ((capacity - 1) / ((modulus - 1) * (modulus - 1)) / f.0.len() as u128)
            .clamp(1, usize::MAX as u128) as usize
    }

    fn transform_ntt(t: Self::T, len: usize) -> Self::F {
        let npot = len.max(1).next_power_of_two();
        let f = convert_crt_input(t, npot);
        (
            Convolve::<N1>::transform_ntt(f.0, npot),
            Convolve::<N2>::transform_ntt(f.1, npot),
            Convolve::<N3>::transform_ntt(f.2, npot),
        )
    }

    fn inverse_transform_ntt(f: Self::F, len: usize) -> Self::T {
        reconstruct_mint_crt((
            Convolve::<N1>::inverse_transform_ntt(f.0, len),
            Convolve::<N2>::inverse_transform_ntt(f.1, len),
            Convolve::<N3>::inverse_transform_ntt(f.2, len),
        ))
    }

    fn ntt_doubling(f: Self::F, monic: bool) -> Self::F {
        if monic {
            let n = f.0.len();
            let mut coefficients = Self::inverse_transform_ntt(f, n);
            coefficients[0] -= MInt::<M>::one();
            coefficients.push(MInt::<M>::one());
            Self::transform_ntt(coefficients, n * 2)
        } else {
            (
                Convolve::<N1>::ntt_doubling(f.0, false),
                Convolve::<N2>::ntt_doubling(f.1, false),
                Convolve::<N3>::ntt_doubling(f.2, false),
            )
        }
    }

fn max_product_sum_count(_f: &Self::F) -> usize

Maximum number of products that can be summed before reconstruction. Both factors must transform canonical coefficients at the supplied transform’s length, and each cyclic product must itself be reconstructible.

Examples found in repository

crates/competitive/src/math/formal_power_series/formal_power_series_impls.rs (line 378)

    fn sum_products(f: &[C::F], g: &[C::F], len: usize) -> Vec<T>
    where
        C: NttReuse<T = Vec<T>>,
        C::F: Clone,
    {
        let chunk = C::max_product_sum_count(&f[0]);
        f.rchunks(chunk)
            .zip(g.chunks(chunk))
            .map(|(f, g)| {
                let mut sum = f[f.len() - 1].clone();
                C::multiply_prefix(&mut sum, &g[0]);
                for (f, g) in f.iter().rev().skip(1).zip(&g[1..]) {
                    C::multiply_add(&mut sum, f, g);
                }
                C::inverse_transform_ntt(sum, len)
            })
            .reduce(|mut sum, part| {
                for (sum, value) in sum.iter_mut().zip(part) {
                    *sum += value;
                }
                sum
            })
            .unwrap()
    }

fn power_projection_step( p_flat: Self::T, q_flat: Self::T, n: usize, py: usize, qy: usize, ) -> (Self::T, Self::T)

Examples found in repository

crates/competitive/src/math/formal_power_series/formal_power_series_impls.rs (line 1290)

    pub fn power_projection(&self, w: &[T], m: usize) -> Self
    where
        C: NttReuse<T = Vec<T>>,
    {
        if w.is_empty() {
            return Self::zeros(m);
        }
        if m <= 1 {
            return Self::from_vec(vec![w[0].clone(); m]);
        }

        let n0 = w.len();
        let mut n = n0.next_power_of_two();
        let mut f = self.prefix_ref(n);
        f.resize(n);

        let base = n * 2;
        let mut p_flat = vec![T::zero(); base];
        for (i, wi) in w.iter().enumerate() {
            p_flat[n - 1 - i] = wi.clone();
        }
        let mut q_flat = vec![T::zero(); base * 2];
        q_flat[0] = T::one();
        let q_offset = base;
        for (i, fi) in f.iter().enumerate() {
            q_flat[q_offset + i] = -fi.clone();
        }
        let mut py = 1usize;
        let mut qy = 2usize;

        let y_limit = m;
        while n > 1 {
            let (mut p, mut q) = C::power_projection_step(p_flat, q_flat, n, py, qy);
            let new_py = (py + qy - 1).min(y_limit);
            let new_qy = (qy + qy - 1).min(y_limit);
            p.resize_with(n * new_py, T::zero);
            q.resize_with(n * new_qy, T::zero);

            let n2 = n / 2;
            for row in p.chunks_exact_mut(n) {
                row[n2..].fill_with(T::zero);
            }
            for row in q.chunks_exact_mut(n) {
                row[n2..].fill_with(T::zero);
            }
            p_flat = p;
            q_flat = q;
            py = new_py;
            qy = new_qy;
            n = n2;
        }

        let base = 2;
        let mut p_y = Vec::with_capacity(py);
        for y in 0..py {
            p_y.push(p_flat[base * y].clone());
        }
        let mut q_y = Vec::with_capacity(qy);
        for y in 0..qy {
            q_y.push(q_flat[base * y].clone());
        }
        (Self::from_vec(p_y) * Self::from_vec(q_y).inv(m)).prefix(m)
    }

Dyn Compatibility

This trait is not dyn compatible.

In older versions of Rust, dyn compatibility was called "object safety".

Implementors

impl<M, N1, N2, N3> NttReuse for Convolve<(M, (N1, N2, N3))>

where M: MIntConvert + MIntConvert<u32>, N1: Montgomery32NttModulus, N2: Montgomery32NttModulus, N3: Montgomery32NttModulus,

impl<M> NttReuse for Convolve<M>

where M: Montgomery32NttModulus,

const MULTIPLE: bool = false