FFT(快速傅里叶变换)学习小记_fft cos sin-优快云博客

本文深入讲解了快速傅里叶变换（FFT）算法的基本原理及其在大数乘法中的应用。介绍了单位复数根的概念，FFT的三个重要引理，并详细解释了离散傅里叶变换（DFT）及逆DFT的过程。同时提供了完整的C++代码实现。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

磕了三天，也算是磕了出来，还需要强化一下。

谢谢yl的ppt。

问题引入：

51nod 大数乘法 V2

题目就是大数乘法，n<=10^5。

一个科学的做法：

把a,b视作多项式。

先分别做点值运算，再做插值运算搞回来就行了。

但是我们知道正常的点值和插值都需要O(n^2)的时间去完成。

FFT就是利用特殊的x去优化这个过程。

单位复数根：

n次单位复数根是满足 $w^n=1$ 的复数。

n次单位复数根恰好有n个，它们均匀地分布在一个半径为1的圆上。

怎么分布？

从(1,0)出发，每次逆时针旋转2π/n弧得到一个新的复数根。

第i个n次单位复数根记作： $w_n^i$

其中 $w_n^i=(cos(2πi/n),sin(2πi/n))$

我们知道复数乘法实际上是旋转。

所以 $w_n^{i+j}=w_n^i*w_n^j$

利用这个性质可以推出FFT需要的三个引理：

1.消去引理：

$w_{n*k}^{i*k}=w_n^i$

2.折半引理：

如果n>0为偶数，那么n次单位复数根的平方的集合就是n/2次单位复数根的集合。

3.求和引理：

n不整除k。

则 $\sum_{j = 1}^{n - 1} w_n^{kj}$
$=\sum_{j = 1}^{n - 1} {\sum_{j = 1}^{n - 1} (w_n^{k})}^j$
${w_n^{k})}^n-1}\over{w_n^k-1}$
${w_n^{n})}^k-1}\over{w_n^k-1}$
${1)}^k-1}\over{w_n^k-1}$
$= 0$

当n整除k时，显然有 $\sum_{j = 1}^{n - 1} w_n^{kj}=n$ 。

离散傅里叶变换（DFT）：

定义结果y为：
这里写图片描述
y就是a的离散傅里叶变换，记为 $y=DFT_n(a)$

FFT计算DFT：

大前提：n是2的整数次幂。（不是必须要补全）。

如何求单个数x的函数值A（x）。

差分为两个式子：
这里写图片描述

$A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$

这样拆有什么变化呢？

1.系数因为奇偶性而分开。

2.次数界变成了n/2

回想一下折半引理，单位复数根的平方并不是n个不同的值，而是由n/2个n/2次单位复数根组成，而每个根恰好出现2次。

所以递归分治是可行的。

变换的式子需要自己推一下，这里不展开，可见code。

FFT求逆DFT：

把DFT想象成矩阵乘法。

那么是不是求出友矩阵的逆矩阵，做一遍DFT就是插值了呢？

友矩阵 $v_n$ 的（j,k）= $w_n^{jk}$

而友矩阵的逆矩阵 $v'_n$ 的（j,k）= $w_n^{-jk}\over n$

这里写图片描述

根据求和引理，当k是n的倍数时，此式子为1。

否则此式子为0。

这不就是单位矩阵I吗？

于是逆DFT就诞生了，和DFT没有什么区别。

改法是：
1.交换对象。
2.单位复数根取逆。
3.结果/n

实现优化：

递归数组赋值实现显然很不靠谱，所以牛逼的先人通过找规律的方法搞出了非递归的做法。

具体的找规律过程有些复杂，在这里不讲。

可以通过code和画图来理解规律。

NTT：

代填……

Code(FFT大数乘法)：

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define fo(i, x, y) for(int i = x; i <= y; i ++)
#define ff(i, x, y) for(int i = x; i < y; i ++)
using namespace std;

const double wu = 1e-8, pi = acos(-1);

const int N = 4e5 + 5;

struct Z {
	double x, y;
	Z (double _x = 0, double _y = 0) {x = _x, y = _y;}
};

Z a[N], b[N], t[N];
int a0, b0, n, po;

Z operator +(Z a, Z b) {return Z(a.x + b.x, a.y + b.y);}
Z operator -(Z a, Z b) {return Z(a.x - b.x, a.y - b.y);}
Z operator *(Z a, Z b) {return Z(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

void read(Z *a, int &a0) {
	char c = ' '; for(; c < '0' || c > '9'; c = getchar());
	a0 = 0; for(; c >= '0' && c <= '9'; c = getchar()) a[a0 ++] = c - '0';
	fo(i, 0, (a0 - 1) / 2) swap(a[i], a[a0 - 1 - i]);
}

void dft(Z *a, int sig) {
	ff(i, 0, n) {
		int p = 0;
		for(int tp = i, j = 0; j < po; j ++, tp /= 2) p = p * 2 + tp % 2;
		t[p] = a[i];
	}
	for(int m = 2; m <= n; m *= 2) {
		int h = m / 2;
		ff(i, 0, h) {
			Z w( cos(i * sig * pi / h), sin(i * sig * pi / h));
			for(int j = i; j < n; j += m) {
				int k = j + h;
				Z u = t[j], v = w * t[k];
				t[j] = u + v; t[k] = u - v;
			}
		}
	}
	ff(i, 0, n) a[i] = t[i];
}

int max2(int x) {int p = 1; for(; p < x; p *= 2); return p * 2;}
int ci(int x) {return (x == 1) ? 0 : ci(x / 2) + 1;}

int ans[N];

int main() {、
	read(a, a0); read(b, b0);
	n = max2(max(a0, b0)), po = ci(n);
	dft(a, 1); dft(b, 1);
	ff(i, 0, n) a[i] = a[i] * b[i];
	dft(a, -1);
	ff(i, 0, n) {
		ans[i] += int(a[i].x / n + wu);
		ans[i + 1] += ans[i] / 10;
		ans[i] %= 10;
	}
	while(ans[n] == 0) n --;
	for(int i = n; i >= 0; i --) printf("%d", ans[i]);
}

指针卡常版FFT：

int r[M];
P a[M], b[M], w[M];

void dft(P *a, int tp, int F) {
	int n = 1 << tp;
	ff(i, 0, n) if(r[i] < i) swap(a[r[i]], a[i]);
	for(int h = 1; h < n; h *= 2) {
		int b = n / h / 2;
		for(int j = 0; j < n; j += 2 * h) {
			P A, *l = a + j, *r = a + j + h, *W = w;
			ff(i, 0, h) {
				A = *r * *W, *r = *l - A, *l = *l + A;
				l ++, r ++, W += b;
			}
		}
	}
}
void fft(P *a, P *b, int tp) {
	int n = 1 << tp;
	ff(i, 0, n) r[i] = r[i / 2] / 2 + (i & 1) * n / 2;
	ff(i, 0, n) w[i] = P(cos(2 * pi * i / n), sin(2 * pi * i / n));
	dft(a, tp, 1); dft(b, tp, 1);
	ff(i, 0, n) a[i] = a[i] * b[i];
	fo(i, 1, n / 2) swap(w[i], w[n - i]);
	dft(a, tp, -1); ff(i, 0, n) a[i].x /= n;
}

指针卡常版NTT：

int tp, r[M];
ll a[M], b[M];

void dft(ll *a, int tp, int F) {
	int n = 1 << tp;
	ff(i, 0, n) if(r[i] < i) swap(a[i], a[r[i]]);
	for(int h = 1; h < n; h *= 2) {
		ll wn = ksm(ksm(3, (mo - 1) / (h * 2)), F == 1 ? 1 : mo - 2);
		for(int j = 0; j < n; j += h * 2) {
			ll A, w = 1, *l = a + j, *r = a + j + h;
			ff(i, 0, h) {
				A = *r * w, *r = (*l - A) % mo, *l = (*l + A) % mo;
				l ++, r ++, w = w * wn % mo;
			}
		}
	}
	if(F == -1) {
		ll v = ksm(n, mo - 2);
		ff(i, 0, n) a[i] = (a[i] + mo) * v % mo;
	}
}
void fft(ll *a, ll *b, int tp) {
	int n = 1 << tp;
	ff(i, 0, n) r[i] = r[i / 2] / 2 + (i & 1) * (n / 2);
	dft(a, tp, 1); dft(b, tp, 1);
	ff(i, 0, n) a[i] = a[i] * b[i] % mo;
	dft(a, tp, -1);
}