FFT和FWT原理和代码
曾经写过一篇牛头不对马嘴的傅里叶级数的课程小论文，摘取一点放出来（FFT），FWT只有代码
FFT
原理
傅里叶级数
三角形式
$$f(x)=\frac{a_0}{2}+\sum_{k=1}^{\infty}{(a_k\cos kx+b_k\sin kx)}$$ ，其中
a 0 = 1 π ∫ − π π f ( x ) d x , a n = 1 π ∫ − π π f ( x ) cos n x d x , ( n = 0 , 1 , 2 … ) , b n = 1 π ∫ − π π f ( x ) sin n x d x , ( n = 0 , 1 , 2 … ) . \begin{aligned}
a_0&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x) dx, \\a_n&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\cos nx dx, (n=0,1,2\dots), \\b_n&=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\sin nx dx, (n=0,1,2\dots) .
\end{aligned}
a 0 a n b n = π 1 ∫ − π π f ( x ) d x , = π 1 ∫ − π π f ( x ) cos n x d x , ( n = 0 , 1 , 2 … ) , = π 1 ∫ − π π f ( x ) sin n x d x , ( n = 0 , 1 , 2 … ) .
指数形式
f ( t ) = c 0 + ∑ n = 1 ∞ ( c n e i ω n t + c − n e − i ω n t ) = ∑ n = − ∞ + ∞ c n e i ω n t f(t)=c_{0}+\sum_{n=1}^{\infty}\left(c_{n} e^{i \omega_{n} t}+c_{-n} e^{-i \omega_{n} t}\right)=\sum_{n=-\infty}^{+\infty} c_{n} e^{i \omega_{n} t}
f ( t ) = c 0 + n = 1 ∑ ∞ ( c n e i ω n t + c − n e − i ω n t ) = n = − ∞ ∑ + ∞ c n e i ω n t
傅里叶变换
傅里叶积分定理
对于 $$T$$ 为周期的周期函数 $$f(t)$$ ，考虑 $$T\rightarrow+\infty$$ ，则在 $$[-\frac{T}{2}, \frac{T}{2}]$$ 上的周期函数 $$f(t)$$ 可以表示任意符合傅里叶级数展开条件的函数，此时 $$\omega=\frac{2\pi}{T}\rightarrow 0$$ ，则 $$\omega_n=n\omega=(n+1)\omega-n\omega$$ 是连续变化的，即 $$\omega_n=n\Delta \omega$$ ，则
c n = Δ ω 2 π ∫ − T / 2 T / 2 f ( t ) e − i ω n t d t c_{n}=\frac{\Delta \omega}{2 \pi} \int_{-T / 2}^{T / 2} f(t) e^{-i \omega_{n} t} dt
c n = 2 π Δ ω ∫ − T / 2 T / 2 f ( t ) e − i ω n t d t
代回原式可得，
f ( t ) = 1 2 π ∑ − ∞ + ∞ ( ∫ − T / 2 T / 2 f ( t ) e − i ω n t d t ) e i ω n t Δ ω = 1 2 π ∫ − ∞ + ∞ ( ∫ − ∞ + ∞ f ( t ) e − i ω n t d t ) e i ω n t d ω f(t)=\frac{1}{2 \pi} \sum_{-\infty}^{+\infty}\left(\int_{-T / 2}^{T / 2} f(t) e^{-i \omega_n t} d t\right) e^{i \omega_n t} \Delta \omega=\frac{1}{2 \pi} \int_{-\infty}^{+\infty}\left(\int_{-\infty}^{+\infty} f(t) e^{-i \omega_n t} d t\right) e^{i \omega_n t} d \omega
f ( t ) = 2 π 1 − ∞ ∑ + ∞ ( ∫ − T / 2 T / 2 f ( t ) e − i ω n t d t ) e i ω n t Δ ω = 2 π 1 ∫ − ∞ + ∞ ( ∫ − ∞ + ∞ f ( t ) e − i ω n t d t ) e i ω n t d ω
傅里叶变换
记 $$F(\omega)=\int_{-\infty}^{+\infty} f(t) e^{-i \omega t} d t$$ ，任何一个函数 $$f(t)$$ 都可以通过变换得到 $$F(\omega)$$ ，
ψ [ f ( t ) ] = F ( ω ) = ∫ − ∞ + ∞ f ( t ) e − i ω t d t \psi[f(t)]=F(\omega) =\int_{-\infty}^{+\infty} f(t) e^{-i \omega t} d t
ψ [ f ( t ) ] = F ( ω ) = ∫ − ∞ + ∞ f ( t ) e − i ω t d t
这个对函数的变换称为傅里叶变换，称 $$F(\omega)$$ 是 $$f(t)$$ 的象函数，称 $$f(t)$$ 是 $$F(\omega)$$ 的象原函数，称 $$|F(\omega)|$$ 为 $$f(t)$$ 的振幅谱，$$\arg{F(\omega)} $$ 为函数 $$f(t)$$ 的相位谱。
由以上变换易得傅里叶逆变换
ψ − 1 [ F ( ω ) ] = f ( t ) = 1 2 π ∫ − ∞ + ∞ F ( ω ) e i ω t d ω \psi^{-1}[F(\omega)]=f(t)=\frac{1}{2 \pi} \int_{-\infty}^{+\infty} F(\omega) e^{i \omega t} d \omega
ψ − 1 [ F ( ω ) ] = f ( t ) = 2 π 1 ∫ − ∞ + ∞ F ( ω ) e i ω t d ω
多项式
系数表达法
我们一般将 $$N-1$$ 次多项式记为 $$A(x) =\sum_{n=0}^{N-1}a_n x^n$$ ，它由 $$n$$ 个系数 $$a_0, a_1, \dots, a_n$$ 唯一确定，称系数表达法。
点值表达法
我们还常用另一种记法，令 $$y_i=A(x_i), i=0,1,\dots, N-1$$ ，相当于在多项式的时域上取 $$N$$ 个点： $$(x_0, y_0), (x_1,y_1), \dots, (x_{N-1}, y_{N-1})$$ ，这些点也可以唯一确定该多项式，称点值表达法
傅里叶逆变换（IDFT）
对于 $$N-1$$ 次多项式 $$A(x) =\sum_{n=0}^{N-1}a_n x^n$$ ，取 $$x_k=e^{\frac{2 \pi i k}{N}}=\omega_N^k(k=0,1, \ldots N-1)$$ ，其中 $$i$$ 是虚数单位，设 $$y_{k}=A(x_k)=\sum_{n=0}^{N-1}a_n\omega_N^{kn}$$ ，由前文傅里叶逆变换类比
a n = 1 N ∑ k = 0 N − 1 ( ω N − n ) k y k a_{n}=\frac{1}{N} \sum_{k=0}^{N-1} (\omega_N^{-n})^ky_{k}
a n = N 1 k = 0 ∑ N − 1 ( ω N − n ) k y k
逆变换相当于多项式 $$B(x)=\sum_{k=0}^{N-1}\frac{y_k}{N}x^k$$ 在 $${\omega_N^{-j}}, j=0,1,\cdots,N-1$$ 的点值表示。
单位根
代数中，若 $$w^n=1$$ ，则称 $$w$$ 为 $$n$$ 次单位根。
在复平面上，以原点为圆心，$$1$$ 为半径作单位圆，将单位圆 $$n$$ 等分，从原点以 $$n$$ 等分点为终点作向量，称辐角为正且最小的向量表示的复数 $$\omega_n=e^{\frac{2 \pi i}{n}} $$ 为主 $$n$$ 次单位根，由复数的运算性质得 $$\omega_n^n=1$$ 。
从辐角为正且最小的开始编号，结合复数的运算性质得
ω n k = e 2 π i k n = cos ( 2 π k n ) + i ⋅ sin ( 2 π k n ) , k = 0 , 1 , … n − 1 , \omega_n^k=e^{\frac{2 \pi i k}{n}}=\cos \left(\frac{2 \pi k}{n}\right)+i \cdot \sin \left(\frac{2 \pi k}{n}\right), \quad k=0,1, \ldots n-1,
ω n k = e n 2 π i k = cos ( n 2 π k ) + i ⋅ sin ( n 2 π k ) , k = 0 , 1 , … n − 1 ,
引理1：
ω d n d k = ( e 2 π i d n ) d k = e 2 π i k n = w n k \omega_{dn}^{d k}=\left(e^{\frac{2\pi i}{d n}}\right)^{d k}=e^{\frac{2 \pi i k}{n}}=w_{n}^{k}
ω d n d k = ( e d n 2 π i ) d k = e n 2 π i k = w n k
引理2：
ω n n 2 = ( e 2 π i n ) n 2 = e π i = − 1 \omega_{n}^{\frac{n}{2}}=\left(e^{\frac{2\pi i}{n}}\right)^{\frac{n}{2}}=e^{\pi i}=-1
ω n 2 n = ( e n 2 π i ) 2 n = e π i = − 1
引理3：$$n$$ 为偶数，$$n$$ 个 $$n$$ 次单位复数根组成的集合即为 $$\frac{n}{2}$$ 个 $$\frac{n}{2}$$ 次单位复数根组成的集合。
证明：即证$$0\leq k\leq n-1, 0\leq j\leq \frac{n}{2}-1$$ 且 $$k, j\in N$$ ， $${(\omega_n^{k})^2}={\omega_{\frac{n}{2}}^{j}}$$.
由引理1和引理2易得 $$\omega_n^{j+\frac{n}{2}}=\omega_n^j\omega_n^{\frac{n}{2}}=-\omega_n^j$$ ， 因此 $${(\omega_n^{k})^2}={(\omega_{n}^{j})^2}={\omega_{n}^{2j}}={\omega_{\frac{n}{2}}^{j}}$$ ，证毕。
实现
快速求解多项式
分治算法
令 $$A_{0}(x)=a_{0}+a_{2} x+a_{4} x^{2}+\cdots a_{n-2} x^{n / 2-1},A_{1}(x)=a_{1}+a_{3} x+a_{5} x^{2}+\cdots a_{n-1} x^{n / 2-1}$$,
则 $$ A(x)=a_{0}+a_{1} x+a_{2} x^{2}+a_{3} x^{3} \cdots a_{n-1} x^{n-1} =A_{0}\left(x^{2}\right)+x A_{1}\left(x^{2}\right)$$
求 $$ A(x)$$ 在 $${\omega_N^k}, ~k=0, 1, \cdots, N-1$$ 处的值相当于求 $$A_0$$ 和 $$A_1$$ 在 $${(\omega_N^k)^2}, ~k=0, 1, \cdots, N-1$$ 的值，由引理3得，上述操作相当于求 $$A_0$$ 和 $$A_1$$ 在 $${\omega_{\frac{N}{2}}^{j}}, 0\leq j\leq \frac{N}{2}-1$$ 的值，于是将问题规模缩小呈原来的一半，分别对二者在 $$\frac{N}{2}$$ 上进行FFT，直到只剩一个常数项
算法复杂度为 $$O(n \log n)$$
特别说明，对于 $$N$$ 非 $$2$$ 的整次幂的FFT，只需补充为 $$0$$ 的项即可，不影响结果。
另外，通过逆变换的式子我们能看出，它也可以使用FFT分治算法计算，只不过点值取了倒数。
朴素的快速傅里叶变换分治算法的实现代码如下：
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 typedef complex<double > comp; const int MAXN = 1000005 ;const comp I (0 , 1 ) const double PI = acos (-1 ); comp A[MAXN * 3 ], B[MAXN * 3 ], tmp[MAXN * 3 ], ans[MAXN * 3 ]; void fft (comp F[], int N, int sgn = 1 ) if (N == 1 ) return ; memcpy (tmp, F, sizeof (comp) * N); for (int i = 0 ; i < N; i++) if (i&1 ) *(F + i / 2 + N / 2 ) = tmp[i]; else *(F+i/2 ) = tmp[i]; fft (F, N / 2 , sgn), fft (F + N / 2 , N / 2 , sgn); comp Wn=exp (2 * PI / N * sgn * I), w=1 ; int n=N>>1 ; for (int i=0 ;i<n;i++,w*=Wn) { tmp[i]=F[i]+w*F[i+n], tmp[i+n]=F[i]-w*F[i+n]; } memcpy (F, tmp, sizeof (comp) * N); }
蝴蝶变换
分治算法已经有了理想的复杂度，但是算法实现时有反复的数组复制操作，我们尝试优化常数，考虑能否提前确定序列中每个元素最终的位置。以 $$N=8$$ 为例：
观察易得，最终序列的元素下标是原序列元素下标二进制的翻转，又称蝴蝶变换（暂未找到严格证明），由此我们可以提前确定序列中每个元素的位置，迭代进行快速傅里叶变换，不需要递归写法了，优化了常数。
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 typedef complex<double > comp; const int MAXN = 1000005 ;const comp I (0 , 1 ) const double PI = acos (-1 ); comp A[MAXN * 3 ], B[MAXN * 3 ], tmp[MAXN * 3 ], ans[MAXN * 3 ]; int rev[MAXN * 3 ];void fft (comp F[], int N, int sgn = 1 ) int bit = __lg(N); for (int i = 0 ; i < N; ++i) { rev[i] = (rev[i >> 1 ] >> 1 ) | ((i & 1 ) << (bit - 1 )); if (i < rev[i]) swap (F[i], F[rev[i]]); } for (int l = 1 ; l < N; l <<= 1 ) { comp Wn = exp (sgn * PI / l * I); for (int i = 0 ; i < N; i += (l<<1 ) ) { comp w (1 ,0 ) ; for (int k = i; k < i + l; ++k, w *= Wn) { comp g = F[k], h = F[k + l] * w; F[k] = g + h, F[k + l] = g - h; } } } }
应用：多项式乘法（卷积）
对于$$A(x) =\sum_{n=0}^{N-1}a_n x^n, B(x) =\sum_{n=0}^{N-1}b_n x^n$$ ，有
C ( x ) = A ( x ) B ( x ) = ∑ i = 0 2 n − 1 C i x i , 其中 C i = ( ∑ k = 0 i a k b i − k ) C(x)=A(x)B(x)=\sum_{i=0}^{2 n-1} C_{i} x^{i}, \text { 其中 } C_{i}=\left(\sum_{k=0}^{i} a_{k} b_{i-k}\right)
C ( x ) = A ( x ) B ( x ) = i = 0 ∑ 2 n − 1 C i x i , 其中 C i = ( k = 0 ∑ i a k b i − k )
计算的时间复杂度是 $$O(n^2)$$ ，这种形式系数下，没有优化的方法，考虑改变系数的形式进行优化。在多项式的时域上取 $$N$$ 个点： $$(x_0, A(x_0)), (x_1, A(x_1)), \dots, (x_{N-1}, A(x_{N-1}))$$ 和 $$(x_0, B(x_0)), (x_1, B(x_1)), \dots, (x_{N-1}, B(x_{N-1}))$$ ，可以通过 $$O(n)$$ 的计算得到 $$(x_0, C(x_0)), (x_1, C(x_1)), \dots, (x_{N-1}, C(x_{N-1}))$$ ，因此这 $$N$$ 个点唯一确定多项式 $$C(x)$$ 的系数，由前面知识可得通过傅里叶逆变换（IDFT）可以确定该系数。
实现了快速傅里叶变换，多项式乘法 $$C(x)=A(x)B(x)$$ ，其中C ( x ) C(x) C ( x )
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 int main () int n = read (), m = read (), N = 1 << __lg(n + m + 1 ) + 1 ; for (int i = 0 ; i <= n; ++i) A[i] = read (); for (int i = 0 ; i <= m; ++i) B[i] = read (); fft (A, N), fft (B, N); for (int i = 0 ; i < N; ++i) ans[i] = A[i] * B[i]; fft (ans, N, -1 ); for (int i = 0 ; i <= n + m; ++i) printf ("%d " , int (ans[i].real () / N + 0.1 )); return 0 ; }
FWT
F W T ( C ) = F W T ( A ) × F W T ( B ) FWT(C)=FWT(A)\times FWT(B) F W T ( C ) = F W T ( A ) × F W T ( B )
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 void FWT_or (int *a, int opt) for (int i = 1 ; i < N; i <<= 1 ) for (int p = i << 1 , j = 0 ; j < N; j += p) for (int k = 0 ; k < i; ++k) if (opt == 1 ) a[i + j + k] = (a[j + k] + a[i + j + k]) % MOD; else a[i + j + k] = (a[i + j + k] + MOD - a[j + k]) % MOD; } void FWT_and (int *a, int opt) for (int i = 1 ; i < N; i <<= 1 ) for (int p = i << 1 , j = 0 ; j < N; j += p) for (int k = 0 ; k < i; ++k) if (opt == 1 ) a[j + k] = (a[j + k] + a[i + j + k]) % MOD; else a[j + k] = (a[j + k] + MOD - a[i + j + k]) % MOD; } void FWT_xor (int *a, int opt) for (int i = 1 ; i < N; i <<= 1 ) for (int p = i << 1 , j = 0 ; j < N; j += p) for (int k = 0 ; k < i; ++k) { int X = a[j + k], Y = a[i + j + k]; a[j + k] = (X + Y) % MOD; a[i + j + k] = (X + MOD - Y) % MOD; if (opt == -1 ) a[j + k] = 1ll * a[j + k] * inv2 % MOD, a[i + j + k] = 1ll * a[i + j + k] * inv2 % MOD; } }
其他快速变换（待补）
快速数论变换（NTT）
快速傅里叶变换利用了单位复根的特殊性质，大大减少了运算，但是每个复数系数的实部和虚部是一个正弦及余弦函数，因此大部分系数都是浮点数，计算量会比较大，且浮点数的计算可能会导致误差。在离散正交变换的理论中，已经证明在复数域内，具有循环卷积特性的唯一变换是DFT，所以在复数域中不存在具有循环卷积性质的更简单的离散正交变换，因此提出了以数论为基础的具有循环卷积性质的快速数论变换（NTT）。
MTT
由于NTT仅能针对特定模数，通过数论的相关知识推导，得到找到三个符合要求的模数，在这三个模数意义下分别FFT然后利用中国剩余定理合并一下，即三模数NTT，用发明的人的首字母记作MTT算法。
快速沃尔什变换（FWT）/ 快速莫比乌斯变换（FMT）
快速傅里叶变换可以用来优化卷积，而对于 $$C_{k}=\sum_{i+j=k} A_{i} * B_{j}$$ ，如果换成 $$i|j, i&j, \cdots$$ 等逻辑运算，就不能使用快速傅里叶变换了，但是受其思想启发，对于每种逻辑运算构造不同的运算方法，再还原到结果的系数，且运算中所有系数均为整数，运算速度相对较快。
