多项式乘法

FFT主要用于快速求多项式的乘积。多项式的乘积就叫做卷积

FG来说,显然暴力算法的复杂度是O(nm),而FFT的时间复杂度为O(nlogn)

多项式的性质:用任意n+1个横坐标不同的点,可以唯一确定一个n次多项式。这个性质叫做多项式的点表示法

证明:设这个多项式f=anxn+an1xn1+...+a1x+a0,那么我们代入n+1个点可以获得n+1个方程,这是n+1元一次方程,写出系数行列式就会发现是范德蒙德行列式,由于横坐标不同,所以行列式的值不为0,于是方程有唯一解

H=FG,于是Hn+m次多项式,我们选取n+m+1个点代入F,G之后即可确定H。这就是点表示法的好处,时间复杂度从O(nm)降低到了O(n+m)

于是我们现在要解决的问题是如何快速将一个多项式的系数表示法(也就是我们通常写出来的多项式)与点表示法互相转换

下文的f(x)=an1xn1+...+a1x+a0

系数表示法转点表示法:我们要在取n+1个特殊点上下功夫,我们取复数域上的单位根

假设n=2c(如果不是的话我们补齐),画出复数域上的单位圆:

image

wnk表示将单位圆等分成n份,将x轴逆时针旋转2πkn所得到的复数,其中k的取值为[0,n1]。比如上图,wn0=1wn3=i

wnk具有如下的性质:

1.ij,wniwnj

2.wnk=cos(2kπn)+isin(2kπn)

3.w2n2k=wnk

4.wnk+n2=wnk

现在,我们将f奇数项和偶数项分开,有f=g+h,其中g(x)=an2xn2+an4xn4+...+a2x2+a0,h=an1xn1+an3xn3+...+a1x

g中每个xx2代替可设φ1(x)=an2xn21+an4xn22+...+a2x+a0,将h提取一个x并将剩下的式子的xx2代替可得φ2(x)=an1xn21+an3xn22+...+a1,于是f(x)=φ1(x2)+xφ2(x2)

现在,设k[0,n21],则f(wnk)=φ1(wn2k)+wnkφ2(wn2k)=φ1(wn2k)+wnkφ2(wn2k),f(wnk+n2)=φ1(wn2k)wnkφ2(wn2k)=φ1(wn2k)wnkφ2(wn2k)

也就是说我们要求出将wn0,wn1,...,wnn1代入f的值,只需要递归求解φ1,φ2即可;一共会划分O(logn)层,每层的总复杂度为O(n),所以时间复杂度为O(nlogn)

点表示法转换为系数表示法:设我们现在已经知道了(wn0,f(wn0)),...,(wnn1,f(wnn1)),我们要求f(x),则有ak=i=0n1f(wni)(wnk)in

证明:

i=0n1f(wni)(wnk)i=i=0n1(j=0n1aj(wni)j)(wnk)i=i=0n1(j=0n1aj(wnj)i)(wnk)i=i=0n1j=0n1aj(wnjk)i=j=0n1aji=0n1(wnjk)i=nak

,最后一步成立是因为

i=0n1(wnjk)i={n,j=k1(wnjk)n1wnjk=1(wnn)jk1wnjk=111wnjk=0,jk

于是我们设ρ(x)=i=0n1f(wni)xi,则我们求出ρ(wn0),ρ(wn1),...,ρ(wn(n1))即可,这就是上文讲的傅里叶正变换(系数表示法转化为点表示法)

最后讲一下实现。实践证明,如果用递归来实现的话,常数是非常大的,所以我们一般利用迭代来实现

假设现在n=8,则每一层的划分如下

a0 a1 a2 a3 a4 a5 a6 a7a0 a2 a4 a6 | a1 a3 a5 a7a0 a4|a2 a6|a1 a5|a3 a7a0|a4|a2|a6|a1|a5|a3|a7

我们现在需要快速获得最后一行的顺序。观察,第一行的二进制位与最后一行对应的二进制位刚好为翻转关系。比如第一行的a4对应最后一行的a1,两者的二进制位分别是100001

证明:考虑一个数如何走到最后一层。从第一层到第二层,如果其末尾是1,那么就会被放到右边,否则就会被放到左边;从第二层到第三层,如果其右移一位后末尾是1那么就会被放到右边,否则就会被放到左边;依次类推,每次走完一层之后都会右移一位。于是一个数的每个二进制位决定了其每一步是向左走还是向右走。而每往右走一步,当前下标就会加上当前组数的个数的一半,向左走则不变。比如a5,先向右走,而当前组(指第一层)有8个数,于是下标加上4,再向左走,下标不变,再向右走,而当前组(第三层第三组的{a1,a5})有两个数,于是下标加上1,最终的下标就是0+4+1=5。由上面的过程可知,每次下标加的数就是将当前二进制位放到翻转的位置。于是结论成立

代码见下

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e5+10;
const double PI=acos(-1);
int n,m,len;
complex<double> y[3][N<<2];
int rev[N<<2];
void change(int op)//交换 
{
	for(int i=0;i<len;i++)
	if(i<rev[i]) swap(y[op][i],y[op][rev[i]]);//防止交换两次 
}
void fft(int op,int on)//像BFS一样从下往上一层一层地合并 
{
	change(op);
	for(int h=2;h<=len;h<<=1)//h为当前这一层的区间长度 
	{
		complex<double> wn(cos(2*PI/h),sin(on*2*PI/h));//将单位圆分成h份 
		for(int j=0;j<len;j+=h)//合并起始点 
		{
			complex<double> w(1,0);
			for(int k=j;k<j+h/2;k++)
			{
				complex<double> u=y[op][k],t=w*y[op][k+h/2];
				y[op][k]=u+t,y[op][k+h/2]=u-t;
				//蝴蝶变换,见OI-wiki 
				w*=wn;
			}
		}
	}
}
int main()
{
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n+1;i++)//千万注意这里要加一 
	{
		int a;
		scanf("%d",&a);
		complex<double> u(a,0);//complex是C++自带的一个类 
		y[0][i-1]=u;
	}
	for(int i=1;i<=m+1;i++)
	{
		int a;
		scanf("%d",&a);
		complex<double> u(a,0);
		y[1][i-1]=u;
	}
	len=n+m+1;
	for(int i=0;;i++)//补齐 
	if(len==(1<<i)) break;
	else if(len<(1<<i))
	{
		len=1<<i;
		break;
	}
	for(int i=0;i<len;i++)//这一步操作见OI-wiki的位逆序置换O(n)实现 
	{
		rev[i]=rev[i>>1]>>1;
		if(i&1) rev[i]|=len>>1;
	}
	fft(0,1),fft(1,1);
	for(int i=0;i<len;i++)
	y[2][i]=y[0][i]*y[1][i];
	fft(2,-1);
	for(int i=0;i<=n+m;i++)
	printf("%d ",(int)(real(y[2][i])/len+0.5));//这里都要这么写,不然有精度误差 
    return 0;
}
posted @   最爱丁珰  阅读(21)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
点击右上角即可分享
微信分享提示