Function convolve_i64_naive 

fn convolve_i64_naive(a: Vec<i64>, b: Vec<i64>, len: usize) -> Vec<i64>

Examples found in repository

crates/competitive/src/math/fast_fourier_transform.rs (line 330)

    pub unsafe fn convolve_i64_avx2(a: Vec<i64>, b: Vec<i64>, len: usize) -> Vec<i64> {
        super::convolve_i64_naive(a, b, len)
    }
    #[target_feature(enable = "avx512f,avx512dq,avx512cd,avx512bw,avx512vl")]
    pub unsafe fn convolve_i64_avx512(a: Vec<i64>, b: Vec<i64>, len: usize) -> Vec<i64> {
        super::convolve_i64_naive(a, b, len)
    }
}

fn bit_reverse<T>(f: &mut [T]) {
    let mut ip = vec![0u32];
    let mut k = f.len();
    let mut m = 1;
    while 2 * m < k {
        k /= 2;
        for j in 0..m {
            ip.push(ip[j] + k as u32);
        }
        m *= 2;
    }
    if m == k {
        for i in 1..m {
            for j in 0..i {
                let ji = j + ip[i] as usize;
                let ij = i + ip[j] as usize;
                f.swap(ji, ij);
            }
        }
    } else {
        for i in 1..m {
            for j in 0..i {
                let ji = j + ip[i] as usize;
                let ij = i + ip[j] as usize;
                f.swap(ji, ij);
                f.swap(ji + m, ij + m);
            }
        }
    }
}

fn real_twiddles(n: usize, inverse: bool, mut f: impl FnMut(usize, Complex<f64>)) {
    const BLOCK: usize = 256;
    let sign = if inverse { 1.0 } else { -1.0 };
    let step = Complex::primitive_nth_root_of_unity(sign * n as f64);
    for start in (1..n / 4).step_by(BLOCK) {
        let mut w = Complex::polar(1.0, sign * std::f64::consts::TAU * start as f64 / n as f64);
        for k in start..(start + BLOCK).min(n / 4) {
            f(k, w);
            w *= step;
        }
    }
}

pub fn transform_real(t: impl IntoIterator<Item = f64>, len: usize) -> Vec<Complex<f64>> {
    let n = len.max(4).next_power_of_two();
    let mut f = vec![Complex::zero(); n / 2];
    for (i, t) in t.into_iter().enumerate() {
        if i & 1 == 0 {
            f[i / 2].re = t;
        } else {
            f[i / 2].im = t;
        }
    }
    fft(&mut f);
    bit_reverse(&mut f);
    f[0] = Complex::new(f[0].re + f[0].im, f[0].re - f[0].im);
    f[n / 4] = f[n / 4].conjugate();
    real_twiddles(n, false, |k, wk| {
        let c = wk.conjugate().transpose() + 1.;
        let d = c * (f[k] - f[n / 2 - k].conjugate()) * 0.5;
        f[k] -= d;
        f[n / 2 - k] += d.conjugate();
    });
    f
}

pub fn inverse_transform_real(mut f: Vec<Complex<f64>>, len: usize) -> Vec<f64> {
    let n = len.max(4).next_power_of_two();
    assert_eq!(f.len(), n / 2);
    f[0] = Complex::new((f[0].re + f[0].im) * 0.5, (f[0].re - f[0].im) * 0.5);
    f[n / 4] = f[n / 4].conjugate();
    real_twiddles(n, true, |k, wk| {
        let c = wk.transpose().conjugate() + 1.;
        let d = c * (f[k] - f[n / 2 - k].conjugate()) * 0.5;
        f[k] -= d;
        f[n / 2 - k] += d.conjugate();
    });
    bit_reverse(&mut f);
    ifft(&mut f);
    let inv = 1. / (n / 2) as f64;
    (0..len)
        .map(|i| inv * if i & 1 == 0 { f[i / 2].re } else { f[i / 2].im })
        .collect()
}

#[inline(always)]
fn convolve_i64_naive(a: Vec<i64>, b: Vec<i64>, len: usize) -> Vec<i64> {
    let (a, b) = if a.len() < b.len() { (b, a) } else { (a, b) };
    if b.len() == 1 {
        return a.into_iter().map(|a| a * b[0]).collect();
    }
    let mut c = vec![0; len];
    for (i, a) in a.chunks(1024).enumerate() {
        for (j, b) in b.iter().enumerate() {
            let start = i * 1024 + j;
            for (c, a) in c[start..start + a.len()].iter_mut().zip(a) {
                *c += *a * *b;
            }
        }
    }
    c
}

impl ConvolveSteps for ConvolveRealFft {
    type T = Vec<i64>;
    type F = Vec<Complex<f64>>;
    fn length(t: &Self::T) -> usize {
        t.len()
    }
    fn transform(t: Self::T, len: usize) -> Self::F {
        transform_real(t.into_iter().map(|t| t as f64), len)
    }
    fn inverse_transform(f: Self::F, len: usize) -> Self::T {
        inverse_transform_real(f, len)
            .into_iter()
            .map(|value| value.round() as i64)
            .collect()
    }
    fn convolve(a: Self::T, b: Self::T) -> Self::T {
        let len = (a.len() + b.len()).saturating_sub(1);
        // Keep accumulation overflow-free and exact in the FFT's f64 representation.
        if (a.len().min(b.len()) <= 32 || {
            let size = len.next_power_of_two();
            let log = size.ilog2();
            let limit = crate::avx_helper!(@dispatch simd_backend, SimdBackend;
                3 * log + 16, 2 * log + 16, 4 * log + 16
            );
            2 * a.len() as u128 * b.len() as u128 <= limit as u128 * size as u128
        }) && a.iter().map(|x| x.unsigned_abs()).max().unwrap_or(0) as u128
            * b.iter().map(|x| x.unsigned_abs()).max().unwrap_or(0) as u128
            <= (1u128 << 53) / a.len().min(b.len()).max(1) as u128
        {
            return crate::avx_helper!(@dispatch simd_backend, SimdBackend;
                unsafe {
                    if a.len().max(b.len()) < 32 {
                        simd::convolve_i64_avx2(a, b, len)
                    } else {
                        simd::convolve_i64_avx512(a, b, len)
                    }
                },
                unsafe { simd::convolve_i64_avx2(a, b, len) },
                convolve_i64_naive(a, b, len)
            );
        }
        if !a.is_empty() && !b.is_empty() {
            crate::avx_helper!(@dispatch_avx2_fma return unsafe {
                simd::convolve_f64_avx2(
                    a.into_iter().map(|x| x as f64),
                    b.into_iter().map(|x| x as f64),
                    0..len,
                )
                .into_iter()
                .map(|x| x.round() as i64)
                .collect()
            }, ());
        }
        let mut a = Self::transform(a, len);
        let b = Self::transform(b, len);
        Self::multiply(&mut a, &b);
        Self::inverse_transform(a, len)
    }