FFT - 快速傅里叶变换

来源：互联网发布：java获取python返回值编辑：程序博客网时间：2024/04/29 01:34

· 博主语文体育老师教的

· 本文年龄限定 16+

· 吐槽上面 2 条的都是⑨

老实说我一点也不想写这东西……

但Delayyy君都写掉了那你说我能不写么> <~

需要深入了解的话推荐：

1、《算法合集 —— 多项式乘法》( 张嘉琳 )

2、《算法导论》关于快速傅里叶变换的章节.

3、《一切皆整数的 FFT》关于点值 X 取整数的方法 (百度空间 AekdyCoin)

自己就把关于复数的证明写了……纠结了半天的东西.

设 W ( k, n ) 表示 1 的 N 次单位根 ( 注意在复数情况下是正好有 n 个解的 ). 则 W (k, n) 具有以下性质：

W (k, n) = cos (2 * π * k / n) + i * sin(2 * π * k / n)

折半性 : W (k, n) = W (2k, 2n) (代入上式系数2 直接约掉)

周期性 : W (k, n) = W ( k + n, n) (三角函数的周期是 n )

　　半反性 : W (k, n) = - W(k + n / 2, n) (同上, 自己YY个三角函数就懂了)

则对于 FFT 主算法中的式子证明 :

S(x) = G(x) + x · F(x)

S(w(k, n)) = G(w(k, n)²) + w(k, n) · F(w(k, n)²)

S(w(k, n)) = G(w(k, n / 2)) + w(k, n) · F(w(k, n / 2))

并且:

S(w(k + n / 2, n)) = G(w(k, n / 2)) + w(k + n / 2, n) · F(w(k, n / 2))

S(w(k + n / 2, n)) = G(w(k, n / 2)) - w(k + n / 2, n) · F(w(k, n / 2))

另外, 弄完主算法后的插值的证明就不放了……如此巨坑还是请自行Google

以上证明完毕. FFT的分治还是很强大的......

参考： CSDN evil live 《FFT高精乘法》

CSDN delayyy 《FFT快速傅里叶变换》

以下是很挫的Code ( 用 FFT 优化高精度乘法 )

#include <cmath>#include <cstdio>#include <cstdlib>#include <algorithm>const double pi = acos((double)-1);using namespace std;char s[10010];long long c[10010];int len_a, len_b, n, t;struct line{    double x, y;}   a[10010], b[10010], w[10010], y[10010];line operator - (line &a, line &b) { line c; c.x = a.x - b.x;  c.y = a.y - b.y;  return c; }line operator + (line &a, line &b) { line c; c.x = a.x + b.x;  c.y = a.y + b.y;  return c; }line operator * (line &a, line &b) { line c; c.x = a.x * b.x - a.y * b.y;  c.y = a.x * b.y + a.y * b.x; return c; }void fft(line *a, int s, int t){    int len = n >> t;    if (len == 1)  return;  line wk;    fft(a, s, t + 1);  fft(a, s + (1 << t), t + 1);    for (int k = 0; k < (len >> 1); ++k)      {         int p = s + (k << (t + 1));          wk = w[k << t] * a[p + (1 << t)];          y[k] = a[p] + wk,          y[k + (len >> 1)] = a[p] - wk;      }        for (int k = 0; k < len; ++k)  a[s + (k << t)] = y[k];}void set(char *s, line *a, int &len)                                                //读入, 四位一压{    int lenk, sta, i;    gets(s + 1);      lenk = strlen(s + 1);  sta = lenk % 4;  len = (lenk + 4) / 4 - 1;     for (i = 1; i <= sta; ++i)  a[len].x = a[len].x * 10 + s[i] - '0';    for (; i <= lenk; i += 4)       a[--len].x = (s[i] - '0') * 1000 + (s[i + 1] - '0') * 100 + (s[i + 2] - '0') * 10 + s[i + 3] - '0';    len = (lenk + 3) / 4;}int main(){    freopen("input.txt", "r", stdin);    freopen("output.txt", "w", stdout);    set(s, a, len_a);   set(s, b, len_b);    n = len_a + len_b;    for (t = 1; t < n; t <<= 1);  n = t;    for (int i = 0; i < n; ++i)  w[i].x = cos(pi * i * 2 / n), w[i].y = sin(pi * 2 * i / n);    fft(a, 0, 0);                                                                  //a 化点值过程    fft(b, 0, 0);                                                                  //b 化点值过程    for (int i = 0; i < n; ++i)  a[i] = a[i] * b[i], w[i].y = -w[i].y;             //逆矩阵处理    fft(a, 0, 0);                                                                  //插值过程    for (int i = 0; i < n; ++i)  c[i] += (long long) round(a[i].x / n), c[i + 1] = c[i] / 10000,  c[i] %= 10000;       int top = n;    while (!c[top]  &&  top > 0)  --top;    printf("%d", c[top]);    for (int i = top - 1; i >= 0; --i)  printf("%04d", c[i]);}