多项式和快速傅里叶变换

多项式和快速傅里叶变换

FFT算法

多项式用点值表示

首先设一个多项式 $A(x)$

$A(x)=\sum_{i=0}^{n-1}a_ix^i=a_0+a_1x+a_2x^2+…+a_{n-1}x^{n-1}$

按$A(x)$下标的奇偶性把 $A(x)$ 分成两半，右边再提一个$x$

${A(x)=(a_0+a_2x^2+…+a_{n-2}x^{n-2})+(a_1x+a_3x^3+…+a_{n-1}x^{n-1})}$
${=(a_0+a_2x^2+…+a_{n-2}x^{n-2})+x(a_1+a_3x^2+…+a_{n-1}x^{n-2})}$

两边好像非常相似，于是再设两个多项式 $A_1(x),A_2(x)$，令

$A_1(x)=a_0+a_2x+a_4x^2+…+a_{n-2}x^{​{n \over 2}-1}$

$A_2(x)=a_1+a_3x+a_5x^2+…+a_{n-1}x^{​{n \over 2}-1}$

比较明显得出

$A(x)=A_1(x^2)+xA_2(x^2)$

单位根：
复数$\omega$满足$\omega^n=1$，称$\omega$是n次单位根
$\omega_n^k=\cos{k\over n}2\pi +i\sin{k\over n} 2\pi$
再设 $k< {n\over 2}$ ，把 $\omega^k_n$ 作为 $x$ 代入 $A(x)$（接下来几步变换要多想想单位根的性质）

$A(\omega^k_n)$
$=A_1( (\omega^k_n)^2 )+\omega^k_nA_2( ( \omega^k_n)^2 )$
$=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})$
$=A_1( \omega^k_{n\over2} )+\omega^k_nA_2( \omega^k_{n \over 2} )$

那么对于${A(\omega^{k+{n\over2}}_n)}$的话，代入
${\omega^{k+{n \over 2}}_n}$

$A(\omega^{k+{n\over 2}}_n)=A_1(\omega^{2k+n}_n)+\omega^{k+{n\over 2}}_nA_2(\omega^{2k+n}_n)$

$=A_1(\omega^{2k}_n\omega^n_n)-\omega^k_nA_2(\omega^{2k}_n\omega^n_n)$

$=A_1( \omega_n^{2k} )-\omega_n^kA_2(\omega_n^{2k})$

$=A_1( \omega_{n \over 2}^k )-\omega^k_nA_2(\omega_{n \over 2}^k)$

如果已知 $A_1(\omega^k_{n\over 2})和 A_2(\omega^k_{n\over 2})$的值，我们就可以同时知道 $A(\omega^k_n)$和 $A(\omega^{k+{n\over2}}_n)$的值

现在我们就可以递归分治来搞 FFT 了

FFT逆变换

一个多项式在分治的过程中乘上单位根的共轭复数（即$\omega_n^k=\cos{k\over n}2\pi-i\sin{k\over n} 2\pi$），分治完的每一项除以 n 即为原多项式的每一项系数

注意，fft只能处理$2^k$的项，所以每次fft前，需要找到相乘后大于最高系数的最小的$2^k$作为limit

具体代码实现

递归版
时间复杂度为$nlogn$

#include <bits/stdc++.h>

using namespace std;
const int N=pow(2,18)+100;
const double pi = acos(-1.0);

struct Complex{
 double x, y;
 Complex(double xx = 0, double yy = 0){x = xx, y = yy;}
}a[N];
Complex operator + (Complex a, Complex b){ return Complex(a.x + b.x, a.y + b.y);}
Complex operator - (Complex a, Complex b){ return Complex(a.x - b.x, a.y - b.y);}
Complex operator * (Complex a, Complex b){ return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);}

/**
 * 递归实现快速傅里叶变换（FFT）或逆变换（IFFT）
 * @param limit 当前处理的数组长度（必须为2的幂次）
 * @param a     复数数组指针（输入/输出）
 * @param type  变换类型：1为FFT（正变换），-1为IFFT（逆变换）
 */

void fft(int limit, Complex* a, int type) {
    if (limit == 1) return; 

    int mid = limit >> 1;   
    Complex* a1 = (Complex*)malloc(sizeof(Complex) * mid);  // 存储偶数索引元素
    Complex* a2 = (Complex*)malloc(sizeof(Complex) * mid);  // 存储奇数索引元素

    // 分治：将数组按奇偶索引拆分
    for (int i = 0; i < limit; i += 2) {
        a1[i >> 1] = a[i];      // 偶数位置元素（索引0, 2, 4...）
        a2[i >> 1] = a[i + 1];  // 奇数位置元素（索引1, 3, 5...）
    }

    fft(mid, a1, type);
    fft(mid, a2, type);

    // 合并结果
    Complex Wn(cos(2.0 * pi / limit), type * sin(2.0 * pi / limit)); // 计算旋转因子
    Complex w(1, 0);  // 初始旋转角度为0，对应单位复数1
    for (int i = 0; i < mid; i++, w = w * Wn) {
        Complex tmp = w * a2[i];  // 计算旋转后的奇数部分
        a[i] = a1[i] + tmp;       // 前一半结果：偶项 + 旋转后的奇项
        a[i + mid] = a1[i] - tmp;  // 后一半结果：偶项 - 旋转后的奇项
    }

    free(a1);
    free(a2);
}

非递归版（效率高，但是难理解）

#include <bits/stdc++.h>
using namespace std;
const int N=100200;
const double pi = acos(-1.0);
struct Complex{
    double x, y;
    Complex(double xx = 0, double yy = 0){x = xx, y = yy;}
}a[N],b[N];

int rev[N*2];

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

int limit=1;

void fft(Complex *a, int type){
    Complex Wn, w;
    for(int i = 0; i < limit; i++){
    if(i < rev[i]) swap(a[i], a[rev[i]]);//得到后序列，i < rev[i] 是保证只交换一次 
    } 
    for(int mid = 1; mid < limit; mid <<= 1){//mid 是每次要处理序列长度的一半 
    Wn = Complex{cos(pi / mid), type * sin(pi / mid)};//得到单位根，角度是 2 * pi / 2 * mid， 2被约掉了 
        for(int R = (mid << 1), j = 0; j < limit; j += R){//枚举序列左端点 
            w = Complex{1, 0};//得到单位根的 0 次方 
            for(int k = 0; k < mid; k++, w = w * Wn){//j + k是在序列中的位置，同时得到单位根的 k 次方 
                Complex x = a[j + k], y = a[j + mid + k] * w;//由单位根的性质(1),(2) 推导而来 
                a[j + k] = x + y;
                a[j + mid + k] = x - y;
            }
        }
    }
}