fn real_twiddles(n: usize, inverse: bool, f: impl FnMut(usize, Complex<f64>))Examples found in repository?
crates/competitive/src/math/fast_fourier_transform.rs (lines 396-401)
382pub fn transform_real(t: impl IntoIterator<Item = f64>, len: usize) -> Vec<Complex<f64>> {
383 let n = len.max(4).next_power_of_two();
384 let mut f = vec![Complex::zero(); n / 2];
385 for (i, t) in t.into_iter().enumerate() {
386 if i & 1 == 0 {
387 f[i / 2].re = t;
388 } else {
389 f[i / 2].im = t;
390 }
391 }
392 fft(&mut f);
393 bit_reverse(&mut f);
394 f[0] = Complex::new(f[0].re + f[0].im, f[0].re - f[0].im);
395 f[n / 4] = f[n / 4].conjugate();
396 real_twiddles(n, false, |k, wk| {
397 let c = wk.conjugate().transpose() + 1.;
398 let d = c * (f[k] - f[n / 2 - k].conjugate()) * 0.5;
399 f[k] -= d;
400 f[n / 2 - k] += d.conjugate();
401 });
402 f
403}
404
405pub fn inverse_transform_real(mut f: Vec<Complex<f64>>, len: usize) -> Vec<f64> {
406 let n = len.max(4).next_power_of_two();
407 assert_eq!(f.len(), n / 2);
408 f[0] = Complex::new((f[0].re + f[0].im) * 0.5, (f[0].re - f[0].im) * 0.5);
409 f[n / 4] = f[n / 4].conjugate();
410 real_twiddles(n, true, |k, wk| {
411 let c = wk.transpose().conjugate() + 1.;
412 let d = c * (f[k] - f[n / 2 - k].conjugate()) * 0.5;
413 f[k] -= d;
414 f[n / 2 - k] += d.conjugate();
415 });
416 bit_reverse(&mut f);
417 ifft(&mut f);
418 let inv = 1. / (n / 2) as f64;
419 (0..len)
420 .map(|i| inv * if i & 1 == 0 { f[i / 2].re } else { f[i / 2].im })
421 .collect()
422}