FFT & NTT / 板砖的形象很优雅
普通 FFT,有点小慢
#include<bits/stdc++.h>
#define db double
#define up(i,l,r) for(int i=l; i<=r; ++i)
#define dn(i,r,l) for(int i=r; i>=l; --i)
using namespace std;
const db pi=acos(-1);
const int N=2145141;
int n, m, tr[N];
struct cp { db x, y; } f[N], g[N], sav[N];
typedef const cp ccp;
cp operator+(ccp &a,ccp &b) { return (cp){a.x+b.x,a.y+b.y}; }
cp operator-(ccp &a,ccp &b) { return (cp){a.x-b.x,a.y-b.y}; }
cp operator*(ccp &a,ccp &b) { return (cp){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x}; }
void fft(cp *f,bool op) {
up(i,0,n-1) if(i<tr[i]) swap(f[i],f[tr[i]]);
for(int p=2; p<=n; p<<=1) {
int len=p>>1;
cp tG=(cp){cos(2*pi/p),sin(2*pi/p)};
if(!op) tG.y*=-1;
for(int k=0; k<n; k+=p) {
cp buf=(cp){1,0};
up(l,k,k+len-1) {
cp tt=buf*f[len+l];
f[len+l]=f[l]-tt, f[l]=f[l]+tt;
buf=buf*tG;
}
}
}
}
signed main() {
ios::sync_with_stdio(0);
cin.tie(0);
cin >> n >> m;
up(i,0,n) cin >> f[i].x;
up(i,0,m) cin >> g[i].x;
for(m+=n, n=1; n<=m; n<<=1);
up(i,0,n-1) tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
fft(f,1), fft(g,1);
up(i,0,n-1) f[i]=f[i]*g[i];
fft(f,0);
up(i,0,m) cout << (int)(f[i].x/n+0.49) << ' ';
return 0;
}
三次变两次 FFT,好写,比 fft 快一点,掉精度严重
#include<bits/stdc++.h>
#define up(i,l,r) for(int i=l; i<=r; ++i)
#define dn(i,r,l) for(int i=r; i>=l; --i)
#define db double
using namespace std;
const db pi=acos(-1);
struct cp { db x, y; };
typedef const cp ccp;
cp operator+(ccp &a,ccp &b) { return (cp){a.x+b.x,a.y+b.y}; }
cp operator-(ccp &a,ccp &b) { return (cp){a.x-b.x,a.y-b.y}; }
cp operator*(ccp &a,ccp &b) { return (cp){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x}; }
vector<cp> w[20];
void init(int m) {
for(int i=1; (1<<i)<=m; ++i) {
int n=(1<<i); w[i].resize(n>>1);
up(j,0,(n>>1)-1) w[i][j]=(cp){cos(2*j*pi/n),sin(2*j*pi/n)};
}
}
void dft(vector<cp> &f) {
int n=f.size();
if(n==1) return;
vector<cp> f0(n>>1), f1(n>>1);
up(i,0,n-1) if(i&1) f1[i>>1]=f[i]; else f0[i>>1]=f[i];
dft(f0), dft(f1);
int id=0; while((1<<id)<n) ++id;
up(i,0,(n>>1)-1) {
cp tmp=w[id][i]*f1[i];
f[i]=f0[i]+tmp, f[i+(n>>1)]=f0[i]-tmp;
}
}
void idft(vector<cp> &f) {
int n=f.size();
dft(f), reverse(&f[1],&f[n]);
up(i,0,n-1) f[i].x/=n, f[i].y/=n;
}
signed main() {
ios::sync_with_stdio(0);
cin.tie(0);
int n, m, N=1;
cin >> n >> m;
while(N<=n+m) N<<=1;
init(N); vector<cp> f(N);
up(i,0,n) cin >> f[i].x;
up(i,0,m) cin >> f[i].y;
dft(f);
up(i,0,N-1) f[i]=f[i]*f[i];
idft(f);
up(i,0,n+m) {
int res=f[i].y*0.5+0.1;
cout << res << ' ';
}
return 0;
}
NTT 取模,做不了浮点
#include<bits/stdc++.h>
#define int long long
#define db long double
#define up(i,l,r) for(int i=l; i<=r; ++i)
#define dn(i,r,l) for(int i=r; i>=l; --i)
using namespace std;
const int P=998244353, G=3, N=1350000;
int ksm(int a,int b=P-2) {
int res=1;
for( ; b; b>>=1) {
if(b&1) res=res*a%P;
a=a*a%P;
}
return res;
}
int n, m, tr[N<<1], f[N<<1], g[N<<1], invn, invG=ksm(G);
void NTT(int *f,bool op) {
up(i,0,n-1) if(i<tr[i]) swap(f[i],f[tr[i]]);
for(int p=2; p<=n; p<<=1) {
int len=p>>1, tG=ksm(op?G:invG,(P-1)/p);
for(int k=0; k<n; k+=p) {
int buf=1;
up(l,k,k+len-1) {
int tt=buf*f[len+l]%P;
f[len+l]=f[l]-tt;
if(f[len+l]<0) f[len+l]+=P;
f[l]=f[l]+tt;
if(f[l]>P) f[l]-=P;
buf=buf*tG%P;
}
}
}
}
signed main() {
ios::sync_with_stdio(0);
cin.tie(0);
cin >> n >> m, ++n, ++m;
up(i,0,n-1) cin >> f[i];
up(i,0,m-1) cin >> g[i];
for(m+=n, n=1; n<m; n<<=1);
up(i,0,n-1) tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
NTT(f,1), NTT(g,1);
up(i,0,n-1) f[i]=f[i]*g[i]%P;
NTT(f,0), invn=ksm(n);
up(i,0,m-2) cout << f[i]*invn%P << ' ';
return 0;
}