FFT - 快速傅里叶变换

    · 博主语文体育老师教的

    · 本文年龄限定 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]);
}