抄的 FFT 和 NTT 模板-CSDN博客

本文链接：https://blog.csdn.net/StarrYooSkY/article/details/102484661

洛谷P1919 ： https://www.luogu.org/problem/P1919

FFT ：

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
const int N=1e6+10;
const double Pi=acos(-1.0);
struct cp {
	double x,y;
	cp friend operator + (cp a,cp b) {
		return (cp) {
			a.x+b.x,a.y+b.y
		};
	}
	cp friend operator - (cp a,cp b) {
		return (cp) {
			a.x-b.x,a.y-b.y
		};
	}
	cp friend operator * (cp a,cp b) {
		return (cp) {
			a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x
		};
	}
} A[N],B[N];
int R[N];
int lim,n,m;
int l;
int cnta,cntb;
void FFt(cp *l,double f) {
	for(int i=0; i<lim; i++)
		if(i<R[i])
			std::swap(l[i],l[R[i]]);
	for(int i=1; i<lim; i<<=1) {
		cp Wn;
		Wn=(cp) {
			cos(Pi/i),f*sin(Pi/i)
		};
		for(int j=0; j<lim; j+=(i<<1)) {
			cp W;
			W=(cp) {
				1.00,0.00
			};
			for(int k=0; k<i; k++,W=W*Wn) {
				cp nx,ny;
				nx=l[j+k];
				ny=l[i+j+k]*W;
				l[j+k]=nx+ny;
				l[i+j+k]=nx-ny;
			}
		}
	}
}
char s1[600100],s2[600100];
int ans[600100];
int main() {
	l=0;
	lim=1;
	scanf("%d",&n);--n;
	m=n;
	scanf("%s%s",s1,s2);
	for(int i=0; i<=n; i++) {
		A[i].x=s1[n-i]-'0';
	}
	for(int i=0; i<=m; i++) {
		B[i].x=s2[n-i]-'0';
	}
	while(lim<=(n+m)) {
		lim=lim<<1;
		l++;
	}
	for(int i=0; i<lim; i++) {
		R[i]=(R[i>>1]>>1)|((i&1)<<(l-1));
	}
	FFt(A,1.00);
	FFt(B,1.00);
	for(int i=0; i<=lim; i++) {
		A[i]=A[i]*B[i];
	}
	FFt(A,-1.00);
	for(int i=0; i<=n+m; i++){
		ans[i]+=(int)(A[i].x/lim+0.5);
		ans[i+1]+=ans[i]/10;
		ans[i]%=10;
	}
	int f=0;
	for(int i=n+m+1;i>=0;i--){
		if(ans[i]!=0) f=1;
		if(f) printf("%d",ans[i]);
	}
	return 0;
}

NTT ：

#include<bits/stdc++.h>
#define ll long long 
using namespace std;
const int maxn=2e6+10; 
const int mod=998244353;     //使用NTT需要保证模数mod 为质数
const ll pr=3;  
 //3是998244353的原根，在比赛中请用上面那个算法提前算出....

ll f1[maxn],f2[maxn],U,V;
ll wn[50],R[maxn],N,M,n,m,l;

ll Power(ll bs,ll js){
    ll S = 1 , T = bs;
    while(js){if(js&1)S=S*T%mod; T=T*T%mod; js>>=1;}
    return S;
}

void GetWn(){
    //需要计算floor(log n)个原根
    for(int i = 0; i <= 25; i ++){
        ll tt = 1<<i;
        wn[i] = Power(pr,(mod-1)/tt);
    }return;
}

void NTT(ll P[],int opt){
    for(int i = 0; i < n; i ++)
        if(i < R[i]) swap(P[R[i]],P[i]);
    for(int i = 1,id = 0; i < n; i<<=1){
        id ++;
        for(int j = 0,p = i<<1; j < n; j += p){
            ll w = 1;
            for(int k = 0; k < i; k ++,w = w*wn[id]%mod){
                U = P[j+k]; V = w*P[j+k+i];
                P[j+k] = (U+V)%mod;  P[j+k+i] = ((U-V)%mod+mod)%mod;
            }
        }
    }
    if(opt == -1){
        //caution:反转时是从1开始 for !!!!!
        for(int i = 1; i < n/2; i ++)swap(P[i],P[n-i]);  
        ll inv = Power(n,mod-2);
        for(int i = 0; i < n; i ++)P[i] = P[i]%mod*inv%mod;
    }return;
}
int ans[600100];
int main(){
    //读入数据：
    cin >> N;--N;M=N;
    string s1,s2;
    cin>>s1>>s2;
    for(int i = 0; i <= N; i ++)f1[i]=s1[N-i]-'0';
    for(int i = 0; i <= M; i ++)f2[i]=s2[N-i]-'0';

    //NTT计算：
    m = N+M; l = 0;
    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));
    GetWn();
    NTT(f1,1);  NTT(f2,1);
    for(int i = 0; i < n; i ++)f1[i] = f1[i]*f2[i]%mod;
    NTT(f1,-1);
    
    //转移答案：
    for(int i=0; i<=n+m; i++){
		ans[i]+=f1[i];
		ans[i+1]+=ans[i]/10;
		ans[i]%=10;
	}
    int f=0;
	for(int i=n+m+1;i>=0;i--){
		if(ans[i]!=0) f=1;
		if(f) printf("%d",ans[i]);
	}
    return 0;
}

洛谷P3803 : https://www.luogu.org/problem/P3803

FFT ：

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
const int N=1e7+10;
const double Pi=acos(-1.0);
struct cp {
	double x,y;
	cp friend operator + (cp a,cp b) {
		return (cp) {
			a.x+b.x,a.y+b.y
		};
	}
	cp friend operator - (cp a,cp b) {
		return (cp) {
			a.x-b.x,a.y-b.y
		};
	}
	cp friend operator * (cp a,cp b) {
		return (cp) {
			a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x
		};
	}
} A[N],B[N];
int R[N];
int lim,n,m;
int l;
int cnta,cntb;
void FFt(cp *l,double f) {
	for(int i=0; i<lim; i++)
		if(i<R[i])
			std::swap(l[i],l[R[i]]);
	for(int i=1; i<lim; i<<=1) {
		cp Wn;
		Wn=(cp) {
			cos(Pi/i),f*sin(Pi/i)
		};
		for(int j=0; j<lim; j+=(i<<1)) {
			cp W;
			W=(cp) {
				1.00,0.00
			};
			for(int k=0; k<i; k++,W=W*Wn) {
				cp nx,ny;
				nx=l[j+k];
				ny=l[i+j+k]*W;
				l[j+k]=nx+ny;
				l[i+j+k]=nx-ny;
			}
		}
	}
}
int main() {
	l=0;
	lim=1;
	scanf("%d%d",&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);
	}
	while(lim<=(n+m)) {
		lim=lim<<1;
		l++;
	}
	for(int i=0; i<lim; i++) {
		R[i]=(R[i>>1]>>1)|((i&1)<<(l-1));
	}
	FFt(A,1.00);
	FFt(B,1.00);
	for(int i=0; i<=lim; i++) {
		A[i]=A[i]*B[i];
	}
	FFt(A,-1.00);
	for(int i=0; i<=n+m; i++){
		printf("%d ",(int)(A[i].x/lim+0.5));
	}
	return 0;
}

NTT ：

#include<bits/stdc++.h>
#define ll long long 
using namespace std;
const int maxn=1e7+10; 
const int mod=998244353;     //使用NTT需要保证模数mod 为质数
const ll pr=3;  
 //3是998244353的原根，在比赛中请用上面那个算法提前算出....

ll f1[maxn],f2[maxn],U,V;
ll wn[50],R[maxn],N,M,n,m,l,ans[maxn];

ll Power(ll bs,ll js){
    ll S = 1 , T = bs;
    while(js){if(js&1)S=S*T%mod; T=T*T%mod; js>>=1;}
    return S;
}

void GetWn(){
    //需要计算floor(log n)个原根
    for(int i = 0; i <= 25; i ++){
        ll tt = 1<<i;
        wn[i] = Power(pr,(mod-1)/tt);
    }return;
}

void NTT(ll P[],int opt){
    for(int i = 0; i < n; i ++)
        if(i < R[i]) swap(P[R[i]],P[i]);
    for(int i = 1,id = 0; i < n; i<<=1){
        id ++;
        for(int j = 0,p = i<<1; j < n; j += p){
            ll w = 1;
            for(int k = 0; k < i; k ++,w = w*wn[id]%mod){
                U = P[j+k]; V = w*P[j+k+i];
                P[j+k] = (U+V)%mod;  P[j+k+i] = ((U-V)%mod+mod)%mod;
            }
        }
    }
    if(opt == -1){
        //caution:反转时是从1开始 for !!!!!
        for(int i = 1; i < n/2; i ++)swap(P[i],P[n-i]);  
        ll inv = Power(n,mod-2);
        for(int i = 0; i < n; i ++)P[i] = P[i]%mod*inv%mod;
    }return;
}

int main(){
    //读入数据：
    scanf("%lld%lld",&N,&M);
    for(int i = 0; i <= N; i ++)scanf("%lld",&f1[i]);
    for(int i = 0; i <= M; i ++)scanf("%lld",&f2[i]);

    //NTT计算：
    m = N+M; l = 0;
    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));
    GetWn();
    NTT(f1,1);  NTT(f2,1);
    for(int i = 0; i < n; i ++)f1[i] = f1[i]*f2[i]%mod;
    NTT(f1,-1);
    
    //转移答案：
    for(int i = 0; i <= m; i ++)ans[i] = f1[i];
    for(int i = 0; i <= m; i ++)printf("%lld ",ans[i]);
    return 0;
}

NTT配套的求原根算法：

#include<bits/stdc++.h>
#define ll long long
#define IL inline
#define RG register
using namespace std;

ll prm[1000],tot,N,root;

ll Power(ll bs,ll js,ll MOD){
    ll S = 1,T = bs;
    while(js){
        if(js&1)S = S*T%MOD;
        T = T*T%MOD; 
        js >>= 1; 
    } return S;
}

IL ll GetRoot(RG ll n){
    RG ll tmp = n - 1 , tot = 0;
    for(RG ll i = 2; i <= sqrt(tmp); i ++){
        if(tmp%i==0){
            prm[++tot] = i;
            while(tmp%i==0)tmp /= i;
        }
    }
    if(tmp != 1)prm[++tot] = tmp;            //质因数分解
    for(RG ll g = 2; g <= n-1; g ++){
        bool flag = 1;
        for(RG int i = 1; i <= tot; i ++){     //检测是否符合条件
            if(Power(g,(n-1)/prm[i],n) == 1)
                { flag = 0; break; } 
        }
        if(flag)return g;
    }return 0;                        //无解 
}

int main(){
    cin >> N;
    root = GetRoot(N);
    cout<<root<<endl;
    return 0;
}

HDU 4609： http://acm.hdu.edu.cn/showproblem.php?pid=4609

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=4e5+10;
const double Pi=acos(-1.0);
struct cp {
    double x,y;
    cp friend operator + (cp a,cp b) {
        return (cp) {a.x+b.x,a.y+b.y};
    }
    cp friend operator - (cp a,cp b) {
        return (cp) {a.x-b.x,a.y-b.y};
    }
    cp friend operator * (cp a,cp b) {
        return (cp) {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
    }
} A[N],B[N];
ll R[N];
int lim,n,m;
int l;
int a[N];
ll num[N],sum[N];
void FFt(cp *l,double f) {
    for(int i=0; i<lim; i++)
        if(i<R[i])
            std::swap(l[i],l[R[i]]);
    for(int i=1; i<lim; i<<=1) {
        cp Wn;
        Wn=(cp) {cos(Pi/i),f*sin(Pi/i)};
        for(int j=0; j<lim; j+=(i<<1)) {
            cp W;
            W=(cp) {1.00,0.00};
            for(int k=0; k<i; k++,W=W*Wn) {
                cp nx,ny;
                nx=l[j+k];
                ny=l[i+j+k]*W;
                l[j+k]=nx+ny;
                l[i+j+k]=nx-ny;
            }
        }
    }
}
int main() {
    int t,n;
    scanf("%d",&t);
    while(t--) {
        l=0;
        lim=1;
        scanf("%d", &n);
        memset(sum, 0, sizeof(sum));
        for(int i = 0; i < n; i++) {
            scanf("%d", &a[i]);
            A[a[i]].x++;
        }
        sort(a,a+n);
        int len=a[n-1];
        while(lim<=(len+len)) {
            lim=lim<<1;
            l++;
        }
        for(int i=0; i<lim; i++) {
            R[i]=(R[i>>1]>>1)|((i&1)<<(l-1));
        }
        FFt(A,1.00);
        for(int i=0; i<=lim; i++) {
            A[i]=A[i]*A[i];
        }
        FFt(A,-1.00);
        for(int i=0;i<=len+len;i++){
            num[i]=(ll)(A[i].x/lim+0.5);
        }
        for(int i=0;i<n;i++) num[a[i]+a[i]]--;
        for(int i=0;i<=len+len;i++) num[i]/=2; 
        for(int i=1;i<=len+len;i++) sum[i]=sum[i-1]+num[i];
        ll ans=0;
        for(int i=0;i<n;i++){
            ans+=sum[len+len]-sum[a[i]];
            ans-=1ll*(n-i-1)*i;
            ans-=1ll*(n-i-1)*(n-i-2)/2;
            ans-=n-1;
        }
        printf("%.7lf\n",ans*1.0/(1ll*(n-2)*(n-1)*n/6));
        for(int i=0;i<N;i++){
        	A[i]={0.00,0.00};
		}
    }
    return 0;
}

2019 上海网络赛： https://nanti.jisuanke.com/t/41400

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=4e5+10;
const double Pi=acos(-1.0);
struct cp {
	double x,y;
	cp friend operator + (cp a,cp b) {
		return (cp) {
			a.x+b.x,a.y+b.y
		};
	}
	cp friend operator - (cp a,cp b) {
		return (cp) {
			a.x-b.x,a.y-b.y
		};
	}
	cp friend operator * (cp a,cp b) {
		return (cp) {
			a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x
		};
	}
} A[N],B[N];
ll R[N],sum[N],ans;
int lim,n,m,cnta,cntb,cntc;
int l;
int aa[N],bb[N],cc[N],a[N],b[N],c[N];
void FFt(cp *l,double f) {
	for(int i=0; i<lim; i++)
		if(i<R[i])
			std::swap(l[i],l[R[i]]);
	for(int i=1; i<lim; i<<=1) {
		cp Wn;
		Wn=(cp) {
			cos(Pi/i),f*sin(Pi/i)
		};
		for(int j=0; j<lim; j+=(i<<1)) {
			cp W;
			W=(cp) {
				1.00,0.00
			};
			for(int k=0; k<i; k++,W=W*Wn) {
				cp nx,ny;
				nx=l[j+k];
				ny=l[i+j+k]*W;
				l[j+k]=nx+ny;
				l[i+j+k]=nx-ny;
			}
		}
	}
}
ll solve1(int a[],int b[],int c[]) {
	for(int i = 0; i <= c[n-1]; ++i) sum[i]  = 0;
	for(int i = 0; i < n; ++i){
		for(int j = 0; j < n; ++j){
			sum[a[i]+b[j]]++;
		}
	}
	ans = 0;
	for(int i = 1; i <= c[n-1]; ++i) sum[i] += sum[i-1];
	for(int i = 0; i < n; ++i) ans += sum[c[i]-1];
	return ans;
}
ll solve2(int a[],int b[],int c[]) {
	cnta=cntb=0;
	for(int i=0; i<n; i++) aa[a[i]]++,cnta=max(cnta,a[i]);
	for(int i=0; i<n; i++) bb[b[i]]++,cntb=max(cntb,b[i]);
	l=0;
	lim=1;
	int len=cnta+cntb;
	while(lim<=len) {
		lim=lim<<1;
		l++;
	}
	for(int i=0; i<=cnta; i++) A[i]=(cp) {aa[i],0.00};
	for(int i=cnta+1; i<=lim; i++) A[i]=(cp) {0.00,0.00};
	for(int i=0; i<=cntb; i++) B[i]=(cp) {bb[i],0.00};
	for(int i=cntb+1; i<=lim; i++) B[i]=(cp) {0.00,0.00};
	for(int i=0; i<lim; i++) {
		R[i]=(R[i>>1]>>1)|((i&1)<<(l-1));
	}
	FFt(A,1.00);
	FFt(B,1.00);
	for(int i=0; i<=lim; i++) {
		A[i]=A[i]*B[i];
	}
	FFt(A,-1.00);
	for(int i=1; i<=c[n-1]; i++) sum[i]=(ll)(A[i].x/lim+0.5)+sum[i-1];
	ans=0;
	for(int i=0; i<n; i++) ans+=sum[c[i]-1];
	for(int i=0; i<n; i++) aa[a[i]]=0;
	for(int i=0; i<n; i++) bb[b[i]]=0;
	return ans;
}
int main() {
	int t,sup;
	ll ans;
	scanf("%d",&t);
	for(int Case=1;Case<=t;Case++) {
		ans=0;
		scanf("%d", &n);
		for(int i = 0; i < n; i++) scanf("%d", &a[i]);
		for(int i = 0; i < n; i++) scanf("%d", &b[i]);
		for(int i = 0; i < n; i++) scanf("%d", &c[i]);
		sort(a,a+n),sort(b,b+n),sort(c,c+n);
		printf("Case #%d: ",Case);
		if(n<=1000) {
			ans=1ll*n*n*n-solve1(a,b,c)-solve1(a,c,b)-solve1(b,c,a);
			printf("%lld\n",ans);
		} else {
			ans=1ll*n*n*n-solve2(a,b,c)-solve2(b,c,a)-solve2(a,c,b);
			printf("%lld\n",ans);
		}
	}
	return 0;
}