FFT(快速傅里叶变换)
用途: 计算多项式卷积。
\(A(x)=a_{0}+a_{1}x+a_{2}x^{2}+···+a_{n-1}x^{n-1}\)
\(B(x)=b_{0}+b_{1}x+b_{2}x^{2}+···+b_{n-1}x^{n-1}\)
求 : \(A(x)B(x)\)
前置芝士
多项式性质: 用任意 \(n\) 个多项式函数图像上的点，可唯一确定一个 \(n-1\) 次多项式。
复数,不会的话右转数学选修2-2.
做法
将 \(A(x)B(x)\) 的系数表示法转化为点表示法，
即 \((x_{1},A(x_{1}))\ \ \ (x_{2},A(x_{2}))\ ···\ (x_{n},A(x_{n}))\)
和 \((x_{1},B(x_{1}))\ \ \ (x_{2},B(x_{2}))\ ···\ (x_{n},B(x_{n}))\)
那么 \(A(x)B(x)\) 就可以表示为
再将其转化为系数表示法。
如何转化？
点的横坐标选择复数域上的单位根。
将复平面中以原点为圆心的单位圆等分为 \(n\) 份，以 \(\omega _{n}^{k}\) 表示从 \(x\) 正半轴一次逆时针取 \(n\) 份中的 \(k\) 份， \(\omega _{n}^{k}\) 就被称为 \(n\) 次单位根。
\(n\) 次单位根的性质:
- \(\omega ^{i}_{n} \ne \omega^{j}_{n} \ \ \ \forall i\ne j\)
- \(\omega^{k}_{n}=\cos \frac{2k\pi}{n}+\sin \frac{2k\pi}{n}i\)
- \(\omega^{0}_{n}=\omega^{n}_{n}=1\)
- \(\omega^{2k}_{2n}=\omega^{k}_{n}\)
- \(\omega^{k+\frac{n}{2}}_{n}=-\omega^{k}_{n}\)
那么点表示法选取的横坐标为: \(\omega^{0}_{n}\ \ \omega^{1}_{n}\ \ \omega^{2}_{n}···\omega^{n-1}_{n}\)
现在，需要快速计算:
\((\omega^{0}_{n},A(\omega^{0}_{n})\ \ \ (\omega^{1}_{n},A(\omega^{1}_{n}))\ ···\ (\omega^{n-1}_{n},A(\omega^{n-1}_{n}))\)
快速傅里叶正变换(系数表示法转化为点表示法)
设:
那么有:
\(k\in [0,n-1]\) , 将值域分为两段 \([0,\frac{n}{2}-1]\ ,\ [\frac{n}{2},n-1]\)
如果 \(k\in[0,\frac{n}{2}-1]\) , 那么有 \((\omega^{k}_{n})^{2}=\omega^{2k}_{n}\)
如果 \([\frac{n}{2},n-1]\) , 将这一段全部减去 \(\frac{n}{2}\) , \(k\) 还是遍历 \([0,\frac{n}{2}-1]\) 将 \(k\) 加上 \(\frac{n}{2}\) 就行了。
发现 \(A(\omega^{k}_{n})\ ,\ A(\omega^{k+\frac{n}{2}}_{n})\) 计算非常相似，只要求出 \(A_{1}(\omega^{k}_{\frac{n}{2}})\ ,\ A_{2}(\omega^{k}_{\frac{n}{2}})\) 即可。
这个过程可以迭代 ，每次将区间分为两段 ，最多会有 \(log\ n\) 层 ，总复杂度 \(O(nlog\ n)\) .
迭代划分(假设有 \(2^3\) 项)
初始序列 \(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [x_{0},x_{1},x_{2},x_{3},x_{4},x_{5},x_{6},x_{7}]\)
奇偶划分一次 \(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [x_{0},x_{2},x_{4},x_{6}],[x_{1},x_{3},x_{5},x_{7}]\)
将这两段继续奇偶划分 \(\ \ \ \ [x_{0},x_{4}],[x_{2},x_{6}],[x_{1},x_{5}],[x_{3},x_{7}]\)
接着分 \(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ [x_{0}],[x_{4}],[x_{2}],[x_{6}],[x_{1}],[x_{5}],[x_{3}],[x_{7}]\)
发现最后一层每一位下标是第一层下标二进制反过来(比如第二项 \((1)_{10}=(001)_2\ ,\ (4)_{10}=(100)_2\))
我们可以递推计算有 \(2^{bit}\) 项的多项式最后一层的第 \(i\) 位下标 \(chg_{i}\)。
也就是
chg[i]=(chg[i>>1]>>1)|((i&1)<<(bit-1));
理解一下就是: 对于一个数\(i\)，找到没有 \(i\) 的第 \(0\) 位的数的反转之后的结果，将这个结果右移一位去掉最高位，再把 \(i\) 的第 \(0\) 位放到最高位，就实现了二进制反转.
\(for\ \ example\) : 计算到 \(6\) 时，将 \((6)_{10}=(110)_{2}\) 的二进制表示去掉第 \(0\) 位，也就是 \((3)_{10}=(011)_{2}\) ，\(3\) 的结果之前处理过了，是 \((110)_{2}=(6)_{10}\)，(这个例子好失败啊)，将反转后的结果右移 \(1\) 位 变成了 \((011)_{2}\) ,然后把 \(6\) 的第 \(0\) 位(还是 \(0\)) ，放到这个数的最高位，就变成了 \((011)_2=(3)_{10}\)。
快速傅里叶逆变换(点表示法转化为系数表示法)
有点 \((\omega^{0}_{n},A(\omega^{0}_{n}))\ \ (\omega^{1}_{n},A(\omega^{1}_{n}))···(\omega^{n-1}_{n},A(\omega^{n-1}_{n}))\).
求 \(A(x)=a_{0}+a_{1}x+a_{2}x^2+···+a_{n-1}x^{n-1}\)
设 \(y_{k}=\omega^{k}_{n}\) 则有
证明:
设 \(S(x)=1+x+x^2+···+x^{n-1}\) .
\(k\ne 0\) 时 , 有：
两柿相减,有：
故 \(S(\omega^{k}_{n})=0\)
\(k=0\) 时 , \(S(\omega^{k}_{n})=n\) ,代入显然
综上 : 把 \(S\) 代入柿子可知：
证毕.
设 \(S(x)=1+x+x^2+···+x^{n-1}\) .
设 \(B(x)=y_{0}+y_{1}x^1+y_{2}x^2+···+y_{n-1}x^{n-1}\).
则有 \(c_{i}=B(\omega^{-i}_{n})\) ,然后用上面说的正变换求一下就行。
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define N 3000005
#define LL long long
#define LB long double
using namespace std;
int n,m,bit,tot,chg[N];
const LB pi=3.1415926535897932384626433832795028841;//圆周率也能写 acos(-1)
struct complex
{
LB x,y;
complex operator+ (const complex &tmp) const//复数加
{
return (complex){x+tmp.x,y+tmp.y};
}
complex operator- (const complex &tmp) const//复数减
{
return (complex){x-tmp.x,y-tmp.y};
}
complex operator* (const complex &tmp) const//复数乘
{
return (complex){x*tmp.x-y*tmp.y,x*tmp.y+y*tmp.x};
}
}a[N],b[N];
inline int qr()
{
char ch=0;int w=1,x=0;
while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}
while(ch<='9'&&ch>='0'){x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
return x*w;
}
inline void FFT(complex a[],int opt)
{
for(register int i=0;i<tot;i++)
if(i<chg[i])
swap(a[i],a[chg[i]]);//找到a[i]奇偶位置转化后的值
for(register int mid=1;mid<tot;mid<<=1)//mid表示现在区间长度的一半，区间长度len=mid<<1
{
complex w1=(complex){cos((LB)pi/mid),(LB)opt*sin((LB)pi/mid)};//计算w1;
for(register int i=0;i<tot;i+=(mid<<1))//计算这一层所有长度为len的区间
{
complex wk=(complex){1,0};//wk初始值为w0=(1,0),w[k]=w[k-1]*w1,每次乘上w1就是要求的wk.
for(register int j=0;j<mid;j++,wk=wk*w1)
{
complex x=a[i+j];//所有偶数项
complex y=wk*a[i+mid+j];//所有奇数项乘wk
a[i+j]=x+y;//迭代
a[i+j+mid]=x-y;
}
}
}
}
int main()
{
n=qr();
m=qr();
for(register int i=0;i<=n;i++)
a[i].x=qr();
for(register int i=0;i<=m;i++)
b[i].x=qr();
while((1<<bit)<m+n+1)//计算一下 log n+m 向上取整
bit++;
tot=(1<<bit);
for(register int i=0;i<tot;i++)//递推第bit层奇偶项转换位置
chg[i]=(chg[i>>1]>>1)|((i&1)<<(bit-1));
FFT(a,1);//系数表示法转化为点表示法
FFT(b,1);
for(register int i=0;i<tot;i++)
a[i]=a[i]*b[i];//计算点表示法纵坐标
FFT(a,-1);//点表示法转化为系数表示法
for(register int i=0;i<m+n+1;i++)
printf("%d ",(int)(a[i].x/tot+0.5));//为了防止精度问题加上0.5再向下取整
printf("\n");
return 0;
}
