1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

use std::ops::Deref;

use crate::ring::*;
use crate::seq::*;

pub mod cooley_tuckey;
pub mod bluestein;
pub mod factor_fft;
pub mod complex_fft;

///
/// Trait for objects that can perform a fast fourier transform over some ring. 
/// 
/// Usually fast implementations of FFTs have to store a lot of precomputed data
/// (e.g. powers of roots of unity), hence they should be represented as objects
/// implementing this trait.
/// 
/// # Note on equality
/// 
/// If you choose to implement [`PartialEq`] for an FFTTable, and `F == G`, then
/// `F` and `G` should satisfy the following properties:
///  - `F` and `G` support the same rings
///  - `F.len() == G.len()`
///  - `F.root_of_unity(ring) == G.root_of_unity(ring)` for each supported ring `ring`
///  - `F.unordered_fft_permutation(i) == G.unordered_fft_permutation(i)` for all `i`
/// In other words, `F` and `G` must have exactly the same output for `unordered_fft`
/// (and thus `fft`, `inv_fft`, ...) on same inputs.
/// 
pub trait FFTAlgorithm<R: ?Sized + RingBase> {

    ///
    /// This FFTTable can compute the FFT of arrays of this length.
    /// 
    fn len(&self) -> usize;

    ///
    /// The root of unity used for the FFT. While all primitive `n`-th roots
    /// of unity can be used equally for computing a Fourier transform, the 
    /// concrete one used determines the order of the output values.
    /// 
    /// See also [`FFTAlgorithm::unordered_fft_permutation()`].
    /// 
    fn root_of_unity(&self, ring: &R) -> &R::Element;

    ///
    /// On input `i`, returns `j` such that `unordered_fft(values)[i]` contains the evaluation
    /// at `zeta^j` of values. Here `zeta` is the value returned by [`FFTAlgorithm::root_of_unity()`]
    /// 
    fn unordered_fft_permutation(&self, i: usize) -> usize;

    ///
    /// The inverse of [`FFTAlgorithm::unordered_fft_permutation()`], i.e. for all i, have
    /// `self.unordered_fft_permutation_inv(self.unordered_fft_permutation(i)) == i`.
    /// 
    fn unordered_fft_permutation_inv(&self, i: usize) -> usize;

    ///
    /// Computes the Fourier transform of the given `values` over the specified ring.
    /// The output is in standard order, i.e. the `i`-th output element is the evaluation
    /// of the input at `self.root_of_unity()^-i` (note the `-`, which is standard
    /// convention for Fourier transforms).
    /// 
    /// # Panics
    /// 
    /// This function panics if `values.len() != self.len()`.
    ///
    fn fft<V>(&self, mut values: V, ring: &R)
        where V: SwappableVectorViewMut<R::Element>
    {
        self.unordered_fft(&mut values, ring);
        permute::permute_inv(&mut values, |i| self.unordered_fft_permutation(i));
    }
        
    ///
    /// Computes the Fourier transform of the given `values` over the specified ring.
    /// The output is in standard order, i.e. the `i`-th output element is the evaluation
    /// of the input at `self.root_of_unity()^i`, divided by `self.len()`.
    /// 
    /// # Panics
    /// 
    /// This function panics if `values.len() != self.len()`.
    ///
    fn inv_fft<V>(&self, mut values: V, ring: &R)
        where V: SwappableVectorViewMut<R::Element>
    {
        permute::permute(&mut values, |i| self.unordered_fft_permutation(i));
        self.unordered_inv_fft(&mut values, ring);
    }

    ///
    /// Computes the Fourier transform of the given values, but the output values are arbitrarily permuted
    /// (in a way compatible with [`FFTAlgorithm::unordered_inv_fft()`]).
    /// 
    /// Note that the FFT of a sequence `a_0, ..., a_(N - 1)` is defined as `Fa_k = sum_i a_i z^(-ik)`
    /// where `z` is an N-th root of unity.
    /// 
    fn unordered_fft<V>(&self, values: V, ring: &R)
        where V: SwappableVectorViewMut<R::Element>;
    
    ///
    /// Inverse to [`FFTAlgorithm::unordered_fft()`], with basically the same contract.
    /// 
    fn unordered_inv_fft<V>(&self, values: V, ring: &R)
        where V: SwappableVectorViewMut<R::Element>;
}

impl<T, R: ?Sized + RingBase> FFTAlgorithm<R> for T
    where T: Deref, T::Target: FFTAlgorithm<R>
{
    fn len(&self) -> usize {
        self.deref().len()
    }

    fn root_of_unity(&self, ring: &R) -> &R::Element {
        self.deref().root_of_unity(ring)
    }

    fn unordered_fft_permutation(&self, i: usize) -> usize {
        self.deref().unordered_fft_permutation(i)
    }

    fn unordered_fft_permutation_inv(&self, i: usize) -> usize {
        self.deref().unordered_fft_permutation_inv(i)
    }

    fn fft<V>(&self, values: V, ring: &R)
        where V: SwappableVectorViewMut<R::Element>
    {
        self.deref().fft(values, ring)
    }

    fn inv_fft<V>(&self, values: V, ring: &R)
        where V: SwappableVectorViewMut<R::Element>
    {
        self.deref().inv_fft(values, ring)
    }

    fn unordered_fft<V>(&self, values: V, ring: &R)
        where V: SwappableVectorViewMut<R::Element>
    {
        self.deref().unordered_fft(values, ring)
    }

    fn unordered_inv_fft<V>(&self, values: V, ring: &R)
        where V: SwappableVectorViewMut<R::Element>
    {
        self.deref().unordered_inv_fft(values, ring)
    }
}