您的位置: 首页  >  IT文章  >  【Uoj34】多项式乘法（NTT，FFT） 【Uoj34】多项式乘法（NTT，FFT）

【Uoj34】多项式乘法（NTT，FFT） 【Uoj34】多项式乘法（NTT，FFT）

分类: IT文章 • 2024-01-20 07:40:48

题面

uoj

题解

首先多项式乘法用(FFT)是一个很久很久以前就写过的东西
直接贴一下代码吧。。

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<queue>
#include<complex>
using namespace std;
#define ll long long
#define RG register
#define MAX 300000
inline int read()
{
    RG int x=0,t=1;RG char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*t;
}
const double Pi=acos(-1);
complex<double> a[MAX],b[MAX];
int r[MAX],n,m,l;
void FFT(complex<double> *P,int opt)
{
	for(int i=0;i<n;++i)if(i<r[i])swap(P[i],P[r[i]]);
	for(int i=1;i<n;i<<=1)
	{
		complex<double> W(cos(Pi/i),opt*sin(Pi/i));
		for(int p=i<<1,j=0;j<n;j+=p)
		{
			complex<double> w(1,0);
			for(int k=0;k<i;w*=W,++k)
			{
				complex<double> X=P[j+k],Y=w*P[i+j+k];
				P[j+k]=X+Y;P[i+j+k]=X-Y;
			}
		}
	}
}
int main()
{
	n=read();m=read();
	for(int i=0;i<=n;++i)a[i]=read();
	for(int i=0;i<=m;++i)b[i]=read();
	m+=n;
	for(n=1;n<=m;n<<=1)++l;
	for(int i=0;i<n;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	FFT(a,1);FFT(b,1);
	for(int i=0;i<n;++i)a[i]*=b[i];
	FFT(a,-1);
	for(int i=0;i<=m;++i)printf("%d ",(int)(a[i].real()/n+0.5));
	puts("");
	return 0;
}

我们知道(FFT)中使用单位复根
满足两个引理

[(W_{2n}^k)^2=W_{n}^{k} ]

[W_{n}^k=-W_{n}^{k+n/2} ]

单位复根在算的过程中很容易出现精度的问题
现在要找到一个拥有相同性质的东西能够代替单位复根就好了

主要是第二个性质难找
因为(W_n)是(n)次单位复根
所以：((W_n)^n=1,(W_n)^{n/2}=-1)

其实，这个性质可以被原根满足：
假设(p)的原根是(g)
再膜(p)意义下：
(g^{varphi(p)}=1 o g^{varphi(p)/2}=sqrt {1})
因为原根不存在一个比(varphi(p))小的数使得(g^k=1)
所以(g^{varphi(p)/2}=-1)
我们发现上面的性质也可以满足

所以，把(n)次单位复根可以替换成原根的(varphi(p)/(2^n))来做
这样就解决了小数精度的问题
当然也是用来解决卷积取膜的问题

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<queue>
using namespace std;
#define ll long long
#define RG register
#define MAX 3000000
inline int read()
{
    RG int x=0,t=1;RG char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*t;
}
const int pr=3;
const int MOD=998244353;
const int phi=MOD-1;
int n,m,r[MAX],l;
int a[MAX],b[MAX];
int fpow(int a,int b)
{
	int s=1;
	while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}
	return s;
}
void NTT(int *P,int opt)
{
	for(int i=0;i<n;++i)if(i<r[i])swap(P[i],P[r[i]]);
	for(int i=1;i<n;i<<=1)
	{
		int W=fpow(pr,phi/(i<<1));
		for(int p=i<<1,j=0;j<n;j+=p)
		{
			int w=1;
			for(int k=0;k<i;++k,w=1ll*w*W%MOD)
			{
				int X=P[j+k],Y=1ll*w*P[i+j+k]%MOD;
				P[j+k]=(X+Y)%MOD;P[i+j+k]=(X-Y+MOD)%MOD;
			}
		}
	}
	if(opt==-1)reverse(&P[1],&P[n]);
}
int main()
{
	n=read();m=read();
	for(int i=0;i<=n;++i)a[i]=read();
	for(int i=0;i<=m;++i)b[i]=read();
	m+=n;
	for(n=1;n<=m;n<<=1)++l;
	for(int i=0;i<n;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	NTT(a,1);NTT(b,1);
	for(int i=0;i<n;++i)a[i]=1ll*a[i]*b[i]%MOD;
	NTT(a,-1);
	int inv=fpow(n,MOD-2);
	for(int i=0;i<n;++i)a[i]=1ll*a[i]*inv%MOD;
	for(int i=0;i<=m;++i)printf("%d ",a[i]);puts("");
	return 0;
}

相关推荐