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;
}