多项式与快速傅里叶变换

前言

快速傅里叶变换（Fast Fourier Transform, FFT）是计算序列的离散傅里叶变换（DFT）或其逆变换（IDFT）的一种算法。

一般用于快速计算多项式乘法。

预备知识

单位复数根

n次单位复数根是满足ωn=1的复数ω

显然n次单位复数根恰好有n个

根据复数的指数形式定义：

eiu=cos(u)+i×sin(u)

可以得到如下引理：

1. 消去引理：对任何整数n≥0,k≥0以及d>0，有ωdkdn=ωkn
2. 折半引理：如果n>0为整数，则n个n次单位复数根的平方构成的集合就是n/2个n/2次单位复数根的平方的集合
3. 求和引理：对任意整数n≥1和不能被n整数的非负整数k，有∑n−1j=0(ωkn)j=0

多项式的表示方法

对于一个次数界为n的多项式

A(x)=∑i=0n−1aixi

我们可以用n维向量(a0,a1……an−1)来表示多项式的系数，这就是多项式的系数表示

或者将n个值xi代入多项式，得到n个值yi

那么点集{(x0,y0),(x1,y1)……(xn−1,yn−1)}就是多项式的点值表示

本文中的点值表示，xi默认取n的单位复数根ωin

那么同样用n维向量(y0,y1……yn−1)作为点值表示

DFT与IDFT

将一个多项式从系数表示转化成点值表示的过程称为DFT

它的逆操作称为IDFT

下面讨论如何使用FFT快速实现DFT与IDFT

为了方便起见，假设所有的n都是2的整次幂

对于多项式A(x)=a0x0+a1x1+a2x2……an−1xn−1

我们构造两个新的多项式A[0](x)和A[1](x)：

A[0](x)=a0x0+a2x1+a4x2+……+an−2xn/2−1

A[1](x)=a1x0+a3x1+a5x2+……+an−1xn/2−1

显然有：

A(x)=A[0](x2)+xA[1](x2)

那么就转化成了一半规模的子问题：

考虑k=0,1,2……n/2−1

A(ωkn)=A[0](ω2kn)+ωknA[1](ω2kn)

=A[0](ωkn/2)+ωknA[1](ωkn/2)

A(ωk+n/2n)=A[0](ω2k+nn)+ωk+n/2nA[1](ω2k+nn)

=A[0](ωkn/2)−ωknA[1](ωkn/2)

也就是说，只需要求出A[0](ωkn/2)和A[1](ωkn/2)就可以求出A(x)

特殊地，当n=1时，A(x)的点值表示就是a0

实际应用中，由于递归处理常数很大，一般采用非递归手段

重点在于确定递归后每个元素所在的位置

发现其实就是把二进制位倒过来，从小到大排序

可以在线性时间内得到这个序列。

再来看IDFT

可以通过范德蒙德矩阵推导出DFT−1n(y)：

aj=1n∑k=0n−1ykω−kjn

这和点值表示的定义非常相似：

yk=∑j=0n−1ajxj=∑j=0n−1ajωjk

于是上述方法完全适用于IDFT的求解

略作改动即可

模板：用FFT实现快速高精数乘

#include<cstdio>#include<cmath>#include<cstring>#include<algorithm>#define cl(x,y) memset(x,y,sizeof(x))using namespace std;const int maxn=300005;const double PI=3.14159265358979323846;int n,lena,lenb,stk[maxn];char a[maxn],b[maxn];struct comp{    double r,i;    comp (double _r=0,double _i=0):r(_r),i(_i) {}    comp operator + (const comp&b) {return comp(r+b.r,i+b.i);}    comp operator - (const comp&b) {return comp(r-b.r,i-b.i);}    comp operator * (const comp&b) {return comp(r*b.r-i*b.i,r*b.i+i*b.r);}    comp operator / (const double&b) {return comp(r/b,i/b);} }A[maxn],B[maxn],C[maxn];int rev[maxn];void get_rev(int n){    rev[0]=0;int l=log2(n);    for (int i=1;i<n;i++)     rev[i]=(rev[i>>1]>>1)|((i&1)<<l-1);}void FFT(comp *a,int n,int d){    for (int i=0;i<n;i++) if (rev[i]<i) swap(a[i],a[rev[i]]);    for (int l=2;l<=n;l<<=1){        comp wl(cos(2*PI/l),d*sin(2*PI/l));        for (int k=0;k<n;k+=l){            comp w(1,0),_t,_T;            for (int j=0,tj=l>>1;j<tj;j++,w=w*wl)             _t=a[k+j],_T=w*a[k+j+tj],a[k+j]=_t+_T,a[k+j+tj]=_t-_T;        }    }    if (d<0) for (int i=0;i<n;i++) a[i]=a[i]/n;}int main(){    while (~scanf("%s%s",a,b)){        lena=strlen(a);lenb=strlen(b);        n=1;cl(A,0);cl(B,0);        while (n<lena+lenb) n<<=1; n<<=1;        for (int i=0;i<lena;i++) A[lena-i-1].r=a[i]-48;        for (int i=0;i<lenb;i++) B[lenb-i-1].r=b[i]-48;        get_rev(n);        FFT(A,n,1);FFT(B,n,1);        for (int i=0;i<n;i++) C[i]=A[i]*B[i];        FFT(C,n,-1);cl(stk,0);        for (int i=0;i<n;i++)         stk[i]+=C[i].r+0.5,         stk[i+1]+=stk[i]/10,         stk[i]%=10;        bool fir=1;        for (int i=n-1;i>=0;i--){            if (stk[i]==0&&fir) continue;            if (stk[i]!=0) fir=0;            putchar(stk[i]+48);        }        if (fir) putchar('0');        putchar('\n');    }    return 0;}