用途: 计算多项式卷积。
(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} e omega^{j}_{n} forall i e j)
- (omega^{k}_{n}=cos frac{2kpi}{n}+sin frac{2kpi}{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})))
快速傅里叶正变换(系数表示法转化为点表示法)
设:
那么有:
(kin [0,n-1]) , 将值域分为两段 ([0,frac{n}{2}-1] , [frac{n}{2},n-1])
如果 (kin[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 e 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("
");
return 0;
}