多项式乘法 （快速傅里叶变换）

算法导论里基本看懂了傅里叶变换的整个运行流程。

后面剩下的就是如何规划设计程序，使得程序漂亮又高效，暂时没时间，先做个记录，后面慢慢实现。

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform

https://www.zybuluo.com/397915842/note/37965

网上摘录某段代码，可供参考：

#include<cmath>2.#include<cstdio>3.#include<cstring>4.const int MAXN = 2e5 + 5;5.const double PI = acos(-1.0);6.#define max(a, b) (a) > (b) ? (a) : (b)7.8.class Complex9.{10.public:11.    double real, imag;12.13.    Complex(double real = 0.0, double imag = 0.0)14.    {15.        this->real = real, this->imag = imag;16.    }17.18.    Complex operator - (const Complex &elem) const19.    {20.        return Complex(this->real - elem.real, this->imag - elem.imag);21.    }22.23.    Complex operator + (const Complex &elem) const24.    {25.        return Complex(this->real + elem.real, this->imag + elem.imag);26.    }27.28.    Complex operator * (const Complex &elem) const29.    {30.        return Complex(this->real * elem.real - this->imag * elem.imag, this->real * elem.imag + this->imag * elem.real);31.    }32.33.    void setValue(double real = 0.0, double imag = 0.0)34.    {35.        this->real = real, this->imag = imag;36.    }37.};38.39.Complex A[MAXN], B[MAXN];40.int res[MAXN], len, mlen, len1, len2;41.char str1[MAXN >> 1], str2[MAXN >> 1];42.43.void Swap(Complex &a, Complex &b)44.{45.    Complex tmp = a;46.    a = b;47.    b = tmp;48.}49.50.void Prepare()51.{52.    len1 = strlen(str1), len2 = strlen(str2);53.    mlen = max(len1, len2);54.    len = 1;55.56.    // 将 len 扩大到 2 的整数幂57.    while(len < (mlen << 1))58.        len <<= 1;59.60.    //初始化多项式的系数61.    for(int i = 0; i < len1; ++ i)62.        A[i].setValue(str1[len1 - i - 1] - '0', 0);63.64.    for(int i = 0; i < len2; ++ i)65.        B[i].setValue(str2[len2 - i - 1] - '0', 0);66.67.    // 补 068.    for(int i = len1; i < len; ++ i)69.        A[i].setValue();70.71.    for(int i = len2; i < len; ++ i)72.        B[i].setValue();73.}74.75.//雷德算法 位逆序置换76.void Rader(Complex y[])77.{78.    for(int i = 1, j = len >> 1, k; i < len - 1; ++ i)79.    {80.        if(i < j)81.            Swap(y[i], y[j]);82.83.        k = len >> 1;84.85.        while(j >= k)86.        {87.            j -= k;88.            k >>= 1;89.        }90.91.        if(j < k)92.            j += k;93.    }94.}95.96.//DFT : op == 197.//IDFT : op == -198.void FFT(Complex y[], int op)99.{100.    //先位逆序置换101.    Rader(y);102.103.    // h 为每次要处理的长度， h = 1 时不需处理104.    for(int h = 2; h <= len; h <<= 1)105.    {106.        // Wn = e^(2 * PI / n)，如果是插值，那么 Wn = e^(-2 * PI / n)107.        Complex Wn(cos(op * 2 * PI / h), sin(op * 2 * PI / h));108.109.        for(int i = 0; i < len; i += h)110.        {111.            //旋转因子，初始化为 e^0112.            Complex W(1, 0);113.114.            for(int j = i; j < i + h / 2; ++ j)115.            {116.                Complex u = y[j];117.                Complex t = W * y[j + h / 2];118.                //蝴蝶操作119.                y[j] = u + t;120.                y[j + h / 2] = u - t;121.                //每次更新旋转因子122.                W = W * Wn;123.            }124.        }125.    }126.127.    // 插值的时候要除以 len128.    if(op == -1)129.        for(int i = 0; i < len; ++ i)130.            y[i].real /= len;131.}132.133.//DFT 后将 A 和 B 相应点值相乘，将结果放到 res 里面134.void Convolution(Complex *A, Complex *B)135.{136.    //evaluation137.    FFT(A, 1), FFT(B, 1);138.139.    for(int i = 0; i < len; ++ i)140.        A[i] = A[i] * B[i];141.142.    //interpolation143.    FFT(A, -1);144.145.    for(int i = 0; i < len; ++ i)146.        res[i] = (int)(A[i].real + 0.5);147.}148.149.void Adjustment(int *arr)150.{151.    //次数界为 len，所以不用担心进位不会进到第 len 位152.    for(int i = 0; i < len; ++ i)153.    {154.        res[i + 1] += res[i] / 10;155.        res[i] %= 10;156.    }157.158.    //去除多余的 0159.    while(-- len && res[len] == 0);160.}161.162.void Display(int *arr)163.{164.    for(int i = len; i >= 0; -- i)165.        putchar(arr[i] + '0');166.167.    putchar('\n');168.}169.170.int main()171.{172.    while(gets(str1) && gets(str2))173.    {174.        Prepare();175.        Convolution(A, B);176.        Adjustment(res);177.        Display(res);178.    }179.180.    return 0;181.}