FFT、NTT学习笔记

参考资料

picks

miskcoo

menci

胡小兔

unname

自为风月马前卒

上面是FFT的，学完了就来看NTT吧

原根

例题：luogu3803

fft优化后模板

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
int n, m, lim=1, rev[2100005];
const double PI=acos(-1.0);
struct Complex{
	double x, y;
	Complex(double xx=0.0, double yy=0.0){
		x = xx;
		y = yy;
	}
	Complex operator+(const Complex &u)const{
		return Complex(x+u.x, y+u.y);
	}
	Complex operator-(const Complex &u)const{
		return Complex(x-u.x, y-u.y);
	}
	Complex operator*(const Complex &u)const{
		return Complex(x*u.x-y*u.y, x*u.y+y*u.x);
	}
}a[2100005], b[2100005];
template<typename T> void rn(T &x){
	x = 0;
	char ch=getchar();
	while(ch<'0' || ch>'9')	ch = getchar();
	while(ch>='0' && ch<='9'){
		x = x * 10 + ch - '0';
		ch = getchar();
	}
}
void fft(Complex a[], int opt){
	for(int i=0; i<lim; i++)
		if(i<rev[i])
			swap(a[i], a[rev[i]]);
	for(int i=2; i<=lim; i<<=1){
		int tmp=i>>1;
		Complex wn=Complex(cos(PI*2.0/i), opt*sin(PI*2.0/i));
		for(int j=0; j<lim; j+=i){
			Complex w=Complex(1.0, 0.0);
			for(int k=0; k<tmp; k++){
				Complex tmp1=a[j+k], tmp2=w*a[j+k+tmp];
				a[j+k] = tmp1 + tmp2;
				a[j+k+tmp] = tmp1 - tmp2;
				w = w * wn;
			}
		}
	}
	if(opt==-1)
		for(int i=0; i<lim; i++)
			a[i].x /= lim;
}
int main(){
	cin>>n>>m;
	for(int i=0; i<=n; i++)	rn(a[i].x);
	for(int i=0; i<=m; i++)	rn(b[i].x);
	int tmpcnt=0;
	while(lim<=n+m)	lim <<= 1, tmpcnt++;
	for(int i=0; i<lim; i++)
		rev[i] = (rev[i>>1]>>1) | ((i&1)<<(tmpcnt-1));
	fft(a, 1);
	fft(b, 1);
	for(int i=0; i<lim; i++)
		a[i] = a[i] * b[i];
	fft(a, -1);
	for(int i=0; i<=n+m; i++)
		printf("%d ", (int)(a[i].x+0.5));
	printf("
");
	return 0;
}

NTT

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;
int n, m, a[2100005], b[2100005], lim=1, limcnt, rev[2100005];
const int mod=998244353, gg=3, gi=332748118;
void rn(int &x){
	char ch=getchar();
	x = 0;
	while(ch<'0' || ch>'9')	ch = getchar();
	while(ch>='0' && ch<='9'){
		x = x * 10 + ch - '0';
		ch = getchar();
	}
}
int ksm(int a, int b){
	int re=1;
	while(b){
		if(b&1)	re = (ll)re * a % mod;
		a = (ll)a * a % mod;
		b >>= 1;
	}
	return re;
}
void ntt(int a[], int opt){
	for(int i=0; i<lim; i++)
		if(i<rev[i])
			swap(a[i], a[rev[i]]);
	for(int i=2; i<=lim; i<<=1){
		int tmp=i>>1, wn=ksm(opt==1?gg:gi, (mod-1)/i);
		for(int j=0; j<lim; j+=i){
			int w=1;
			for(int k=0; k<tmp; k++){
				int tmp1=a[j+k], tmp2=(ll)w*a[j+k+tmp]%mod;
				a[j+k] = (tmp1 + tmp2) % mod;
				a[j+k+tmp] = (tmp1 - tmp2 + mod) % mod;
				w = (ll)w * wn % mod;
			}
		}
	}
	if(opt==-1){
		int inv=ksm(lim, mod-2);
		for(int i=0; i<lim; i++)
			a[i] = (ll)a[i] * inv % mod;
	}
}
int main(){
	cin>>n>>m;
	for(int i=0; i<=n; i++)	rn(a[i]);
	for(int i=0; i<=m; i++)	rn(b[i]);
	while(lim<=n+m)	lim <<= 1, limcnt++;
	for(int i=0; i<lim; i++)
		rev[i] = (rev[i>>1]>>1) | ((i&1)<<(limcnt-1));
	ntt(a, 1);
	ntt(b, 1);
	for(int i=0; i<lim; i++)
		a[i] = (ll)a[i] * b[i] % mod;
	ntt(a, -1);
	for(int i=0; i<=n+m; i++)
		printf("%d ", a[i]);
	printf("
");
	return 0;
}

递归版裸fft没什么优化

#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
int n, m;
const double PI=acos(-1.0);
struct Complex{
    double x, y;
    Complex(double xx=0.0, double yy=0.0){
        x = xx;
        y = yy;
    }
    Complex operator+(const Complex &u)const{
        return Complex(x+u.x, y+u.y);
    }
    Complex operator-(const Complex &u)const{
        return Complex(x-u.x, y-u.y);
    }
    Complex operator*(const Complex &u)const{
        return Complex(x*u.x-y*u.y, x*u.y+y*u.x);
    }
}a[4000005], b[4000005], buf[4000005];
void fft(Complex a[], int lim, int opt){
	if(lim==1)	return ;
	int tmp=lim/2;
	for(int i=0; i<tmp; i++){
		buf[i] = a[2*i];
		buf[i+tmp] = a[2*i+1];
	}
	for(int i=0; i<lim; i++)
		a[i] = buf[i];
	fft(a, tmp, opt);
	fft(a+tmp, tmp, opt);
	Complex wn=Complex(cos(PI*2.0/lim), opt*sin(PI*2.0/lim)), w=Complex(1.0, 0.0);
    for(int i=0; i<tmp; i++){
        buf[i] = a[i] + w * a[i+tmp];
        buf[i+tmp] = a[i] - w * a[i+tmp];
        w = w * wn;
    }
    for(int i=0; i<lim; i++)
		a[i] = buf[i];
}
int main(){
	cin>>n>>m;
	for(int i=0; i<=n; i++)	scanf("%lf", &a[i].x);
	for(int i=0; i<=m; i++)	scanf("%lf", &b[i].x);
	int lim=1;
	while(lim<=n+m)	lim <<= 1;
	fft(a, lim, 1);
	fft(b, lim, 1);
	for(int i=0; i<=lim; i++)
		a[i] = a[i] * b[i];
	fft(a, lim, -1);
	for(int i=0; i<=n+m; i++)
		printf("%d ", (int)(a[i].x/lim+0.5));
	printf("
");
	return 0;
}