前置知识: F F T FFT FFT
介绍
一般的NTT的模数 P P P都要满足 P = r × 2 k + 1 P=r\times 2^k+1 P=r×2k+1。但如果不满足这个条件,那就不能直接用NTT了,需要用到任意模数NTT。
可以了解一下三模数NTT。
下面就来介绍一种任意模数NTT——拆系数FFT。
拆系数FFT
对于多项式 f ( x ) f(x) f(x)的系数 A A A,可以将 A A A分为高位 A 1 A_1 A1和低位 A 0 A_0 A0,比如令 A = A 1 × 2 15 + A 0 A=A_1\times 2^{15}+A_0 A=A1×215+A0。那么,可以将 A 1 A_1 A1作为系数构造多项式 f 1 ( x ) f_1(x) f1(x),将 A 0 A_0 A0作为系数构造多项式 f 0 ( x ) f_0(x) f0(x)。同理,将多项式 g ( x ) g(x) g(x)的系数 B B B分为高位 B 1 B_1 B1和低位 B 0 B_0 B0,得到 g 1 ( x ) g_1(x) g1(x)和 g 0 ( x ) g_0(x) g0(x)。
那么, f ( x ) g ( x ) = f 1 ( x ) g 1 ( x ) × 2 30 + ( f 1 ( x ) g 0 ( x ) + f 0 ( x ) g 1 ( x ) ) × 2 15 + f 0 ( x ) g 0 ( x ) f(x)g(x)=f_1(x)g_1(x)\times 2^{30}+(f_1(x)g_0(x)+f_0(x)g_1(x))\times 2^{15}+f_0(x)g_0(x) f(x)g(x)=f1(x)g1(x)×230+(f1(x)g0(x)+f0(x)g1(x))×215+f0(x)g0(x)
用 F F T FFT FFT来处理,正变换需要四次。本来逆变换也需要四次,但因为 f 1 ( x ) g 0 ( x ) f_1(x)g_0(x) f1(x)g0(x)和 f 0 ( x ) g 1 ( x ) f_0(x)g_1(x) f0(x)g1(x)是同阶的,可以相加之后在进行逆变换,所以逆变换只需要三次。
在得出三个部分之后,在乘各个部分的阶数的时候,
整个过程中,总共有 7 7 7次 D F T DFT DFT或 I D F T IDFT IDFT。
时间复杂度为 O ( n log n ) O(n\log n) O(nlogn)。
相对于三模数NTT说,拆系数FFT好理解、好实现、跑得快,但有时候精度可能会出问题,而且如果结果多项式中的系数较大,long double存不下的话,也不能用,在做题时要有所注意。
例题
#include<bits/stdc++.h>
using namespace std;
const int N=500000;
const long double pi=acos(-1.0);
int n,m,l=1;
long long p,f[N+5],g[N+5],w1[N+5],w2[N+5],w3[N+5],ans[N+5];
struct cp{
long double a,b;
cp operator +(const cp ax)const{
return (cp){a+ax.a,b+ax.b};
}
cp operator -(const cp ax)const{
return (cp){a-ax.a,b-ax.b};
}
cp operator *(const cp ax)const{
return (cp){a*ax.a-b*ax.b,b*ax.a+a*ax.b};
}
}w,wn,v1[N+5],v2[N+5],v3[N+5],v4[N+5],h1[N+5],h2[N+5],h3[N+5];
void ch(cp *a){
for(int i=1,j=l/2,k;i<l-1;i++){
if(i<j) swap(a[i],a[j]);
k=l/2;
while(j>=k){
j-=k;k>>=1;
}
j+=k;
}
}
void fft(cp *a,int fl){
ch(a);
for(int i=2;i<=l;i<<=1){
wn=(cp){cos(fl*2*pi/i),sin(fl*2*pi/i)};
for(int j=0;j<l;j+=i){
w=(cp){1,0};
for(int k=j;k<j+i/2;k++,w=w*wn){
cp t=a[k],u=w*a[k+i/2];
a[k]=t+u;
a[k+i/2]=t-u;
}
}
}
}
int main()
{
scanf("%d%d%lld",&n,&m,&p);
for(int i=0;i<=n;i++){
scanf("%lld",&f[i]);f[i]%=p;
}
for(int i=0;i<=m;i++){
scanf("%lld",&g[i]);g[i]%=p;
}
while(l<n+m+1) l<<=1;
for(int i=0;i<l;i++){
v1[i].a=f[i]>>15;
v2[i].a=f[i]&((1<<15)-1);
v3[i].a=g[i]>>15;
v4[i].a=g[i]&((1<<15)-1);
}
fft(v1,1);fft(v2,1);fft(v3,1);fft(v4,1);
for(int i=0;i<l;i++){
h1[i]=v1[i]*v3[i];
h2[i]=v1[i]*v4[i]+v2[i]*v3[i];
h3[i]=v2[i]*v4[i];
}
fft(h1,-1);fft(h2,-1);fft(h3,-1);
for(int i=0;i<l;i++){
w1[i]=(long long)(h1[i].a/l+0.5)%p;
w2[i]=(long long)(h2[i].a/l+0.5)%p;
w3[i]=(long long)(h3[i].a/l+0.5)%p;
ans[i]=(w1[i]*(1<<30)%p+w2[i]*(1<<15)%p+w3[i])%p;
}
for(int i=0;i<=n+m;i++){
printf("%lld ",ans[i]);
}
return 0;
}