FFT/NTT(板子整理)

小衣同学

已于 2022-02-26 17:25:07 修改

阅读量1.3k

点赞数

分类专栏：知识点总结文章标签： FFT NTT

于 2019-07-19 21:57:06 首次发布

本文链接：https://blog.youkuaiyun.com/Code92007/article/details/96361115

版权

知识点总结专栏收录该内容

90 篇文章

订阅专栏

思路来源

快速傅里叶变换 & 快速数论变换 - Candy? - 博客园

FTT板子整理

//fft
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N=(1<<18)+5, INF=1e9;
const double PI=acos(-1);
inline int read(){
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}

struct meow{
    double x, y;
    meow(double a=0, double b=0):x(a), y(b){}
};
meow operator +(meow a, meow b) {return meow(a.x+b.x, a.y+b.y);}
meow operator -(meow a, meow b) {return meow(a.x-b.x, a.y-b.y);}
meow operator *(meow a, meow b) {return meow(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x);}
meow conj(meow a) {return meow(a.x, -a.y);}
typedef meow cd;

struct FastFourierTransform {
    int n, rev[N];
    cd omega[N], omegaInv[N];
    void ini(int lim) {
        n=1; int k=0;
        while(n<lim) n<<=1, k++;
        for(int i=0; i<n; i++) rev[i] = (rev[i>>1]>>1) | ((i&1)<<(k-1));
        
        for(int k=0; k<n; k++) {
            omega[k] = cd(cos(2*PI/n*k), sin(2*PI/n*k));
            omegaInv[k] = conj(omega[k]);
        }
    }
    void fft(cd *a, cd *w) {
        for(int i=0; i<n; i++) if(i<rev[i]) swap(a[i], a[rev[i]]);
        for(int l=2; l<=n; l<<=1) {
            int m=l>>1;
            for(cd *p=a; p!=a+n; p+=l) 
                for(int k=0; k<m; k++) {
                    cd t = w[n/l*k] * p[k+m];
                    p[k+m]=p[k]-t;
                    p[k]=p[k]+t;
                }
        }
    }
    void dft(cd *a, int flag) {
        if(flag==1) fft(a, omega);
        else {
            fft(a, omegaInv);
            for(int i=0; i<n; i++) a[i].x/=n;
        }
    }
    void mul(cd *a, cd *b, int m) {
        ini(m);
        dft(a, 1); dft(b, 1);
        for(int i=0; i<n; i++) a[i]=a[i]*b[i];
        dft(a, -1);
    }
}f;

int n1, n2, m, c[N];
cd a[N], b[N];
char s1[N], s2[N];
int main() {
    //freopen("in","r",stdin);
    scanf("%s%s",s1,s2);
    n1=strlen(s1); n2=strlen(s2);
    for(int i=0; i<n1; i++) a[i].x = s1[n1-i-1]-'0';
    for(int i=0; i<n2; i++) b[i].x = s2[n2-i-1]-'0';
    m=n1+n2-1;
    f.mul(a, b, m);
    for(int i=0; i<m; i++) c[i]=floor(a[i].x+0.5);
    for(int i=0; i<m; i++) c[i+1]+=c[i]/10, c[i]%=10;
    if(c[m]) m++;
    for(int i=m-1; i>=0; i--) printf("%d",c[i]);
}

NTT板子整理

~~什么，都0202年了你的ntt的wn还要现求，快换一个不用现求的板子吧~~

以hdu6900为例，Z是模数类，本题是一道分治/启发式NTT的经典题目

Q：为什么998244353（7*17*2^23+1）模数下N最大为2^23？

A：

原根也符合复数域的单位原根的性质，

但区别在于，fft可以不断细分，因为复数无穷，而ntt的原根不能再细分，

好比说，最后递归长度=2，但是我只剩了一个点值，

求不出这个原函数，dft和idft就不可逆了，即失真

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned ll
const int N = 1<<18, P = 998244353;
const int Primitive_root = 3;//先用Get_root求出来原根然后当const用 
struct Z{
    int x;
    Z(const int _x=0):x(_x){}
    Z operator +(const Z &r)const{ return x+r.x<P?x+r.x:x+r.x-P;}
    Z operator -(const Z &r)const{ return x<r.x?x-r.x+P:x-r.x;}
    Z operator -()const{ return x?P-x:0;}
    Z operator *(const Z &r)const{ return static_cast<ull>(x)*r.x%P;}
    Z operator +=(const Z &r){ return x=x+r.x<P?x+r.x:x+r.x-P, *this;}
    Z operator -=(const Z &r){ return x=x<r.x?x-r.x+P:x-r.x, *this;}
    Z operator *=(const Z &r){ return x=static_cast<ull>(x)*r.x%P, *this;}
    friend Z Pow(Z, int);
    pair<Z,Z> Mul(pair<Z,Z> x, pair<Z,Z> y, Z f)const{
        return make_pair(
            x.first*y.first+x.second*y.second*f,
            x.second*y.first+x.first*y.second
        );
    }
    Z Quadratic_residue()const{
        if(x<=1) return x;
        if(Pow((Z)x, (P-1)/2).x!=1) return -1;
        Z y, f;
        mt19937 rng(20030226);
        do y=rng()%(x-1)+1; while(Pow(f=y*y-x, (P-1)/2).x==1);
        pair<Z,Z> ans=make_pair(1, 0), t=make_pair(y, 1);
        for(int i=(P+1)/2; i; i>>=1, t=Mul(t, t, f)) if(i&1) ans=Mul(ans, t, f);
        return min(ans.first.x, P-ans.first.x);
    }
};
Z Pow(Z x, int y=P-2){
    Z ans=1;
    for(; y; y>>=1, x=x*x) if(y&1) ans=ans*x;
    return ans;
}
namespace Poly{
    Z w[N];
    Z Inv[N];
    vector<Z> ans;
    vector<vector<Z>> p;
    ull F[N];
    int Get_root(){
		static int pr[N],cnt;
		int n=P-1,sz=(int)(sqrt(n)),root=-1;
		for(int i=2;i<=sz;++i){if(n%i==0)pr[cnt++]=i;while(n%i==0)n/=i;}
		if(n>1)pr[cnt++]=n;
		for(int i=1;i<P;++i){
			if(Pow((Z)i,P-1).x==1){
				bool fl=true;
				for(int j=0;j<cnt;++j){
					if(Pow(i,(P-1)/pr[j]).x==1){
						fl=false;break;
					}
				}
				if(fl){root=i;break;}
			}
		}
		return root;
	}
    void Init(){
    	//printf("root:%d\n",Primitive_root=Get_root()); 先求出来原根然后当const用 
        for(int i=1; i<N; i<<=1){
            w[i]=1;
            Z t=Pow((Z)Primitive_root, (P-1)/i/2);
            for(int j=1; j<i; ++j) w[i+j]=w[i+j-1]*t;
        }
        Inv[1]=1;
        for(int i=2; i<N; ++i) Inv[i]=Inv[P%i]*(P-P/i);
    }
    int Get(int x){ int n=1; while(n<=x) n<<=1; return n;}
    int Mod(int x){ return x<P?x:x-P;}
    void DFT(vector<Z> &f, int n){
        if((int)f.size()!=n) f.resize(n);
        for(int i=0, j=0; i<n; ++i){
            F[i]=f[j].x;
            for(int k=n>>1; (j^=k)<k; k>>=1);
        }
        if(n<=4){
            for(int i=1; i<n; i<<=1) for(int j=0; j<n; j+=i<<1){
                Z *W=w+i;
                ull *F0=F+j, *F1=F+j+i;
                for(int k=j; k<j+i; ++k, ++W, ++F0, ++F1){
                    ull t=(*F1)*(W->x)%P;
                    (*F1)=*F0+P-t, (*F0)+=t;
                }
            }
        }
        else{
            for(int j=0; j<n; j+=2){
                int t=F[j+1];
                F[j+1]=Mod(F[j]+P-t), F[j]=Mod(F[j]+t);
            }
            for(int j=0; j<n; j+=4){
                int t0=F[j+2], t1=F[j+3]*w[3].x%P;
                F[j+2]=F[j]+P-t0, F[j]+=t0;
                F[j+3]=F[j+1]+P-t1, F[j+1]+=t1;
            }
            for(int i=4; i<n; i<<=1) for(int j=0; j<n; j+=i<<1){
                Z *W=w+i;
                ull *F0=F+j, *F1=F+j+i;
                for(int k=j; k<j+i; k+=4, W+=4, F0+=4, F1+=4){
                    int t0=(W->x)**F1%P;
                    int t1=(W+1)->x**(F1+1)%P;
                    int t2=(W+2)->x**(F1+2)%P;
                    int t3=(W+3)->x**(F1+3)%P;
                    *F1=*F0+P-t0, *F0+=t0;
                    *(F1+1)=*(F0+1)+P-t1, *(F0+1)+=t1;
                    *(F1+2)=*(F0+2)+P-t2, *(F0+2)+=t2;
                    *(F1+3)=*(F0+3)+P-t3, *(F0+3)+=t3;
                }
            }
        }
        for(int i=0; i<n; ++i) f[i]=F[i]%P;
    }
    void IDFT(vector<Z> &f, int n){
        f.resize(n), reverse(f.begin()+1, f.end()), DFT(f, n);
        Z I=1;
        for(int i=1; i<n; i<<=1) I*=(P+1)/2;
        for(int i=0; i<n; ++i) f[i]*=I;
    }
    vector<Z> operator +(const vector<Z> &f, const vector<Z> &g){
        vector<Z> ans=f;
        ans.resize(max(f.size(), g.size()));
        for(int i=0; i<(int)g.size(); ++i) ans[i]+=g[i];
        return ans;
    }
    vector<Z> operator *(const vector<Z> &f, const vector<Z> &g){
        static vector<Z> F, G;
        F=f, G=g;
        int p=Get(f.size()+g.size()-2);
        DFT(F, p), DFT(G, p);
        for(int i=0; i<p; ++i) F[i]*=G[i];
        IDFT(F, p);
        return F.resize(f.size()+g.size()-1), F;
    }
    vector<Z> operator *(const vector<Z> &f, Z g){
        vector<Z> ans=f;
        for(Z &i:ans) i*=g;
        return ans;
    }
}
using namespace Poly;
int T,n,b[N],c[N];
Z fac[N],ifac[N];
vector<Z>work(int l, int r){
    if(l==r){
		return{c[l],b[l]};
    } 
	int mid=(l+r)>>1;
    return work(l,mid)*work(mid+1,r);
}
void init(int n){
	fac[0]=1;
    for(int i=1;i<=n;++i){
		fac[i]=fac[i-1]*i;
    } 
	ifac[n]=Pow(fac[n]);
    for(int i=n;i;--i){
    	ifac[i-1]=ifac[i]*i;
    }
}
vector<Z>a,d,res;
int main(){
    Init();
    init(N-1);
    scanf("%d",&T);
    while(T--){
        scanf("%d",&n);
        a.resize(n+1);
        for(int i=0;i<=n;++i){
    	    scanf("%d",&a[i].x);
			a[i]*=fac[i];
        } 
		for(int i=2;i<=n;++i){
			scanf("%d",&b[i]);
		}
		for(int i=2;i<=n;++i){
			scanf("%d",&c[i]);
		}
        d=work(2,n);
        reverse(d.begin(),d.end());
        res=d*a;
        for(int i=0;i<=n;++i){
			printf("%d%c",(res[i+n-1]*ifac[i]).x," \n"[i==n]);
		} 
    }
    return 0;
}

一个短很多的板子，感谢cls的板子…

namespace NTT//优化过的ntt，使用时记得初始化。998244353,1004535809,469762049,985661441,167772161. g = 3.950009857,g=7
{
    const int p = 998244353,g = 3;
    int w[maxn<<2],inv[maxn<<2],r[maxn<<2],last;
    int mod(int x){return x >= p ? x - p : x;}
    ll qp(ll base,ll n)
    {
        base %= p;
        ll res = 1;
        while(n){
            if(n&1) (res *= base) %= p;
            (base *= base) %= p;
            n >>= 1;
        }
        return res;
    }

    void init()
    {
        int lim = maxn << 1;//最长数组的两倍
        inv[1] = 1;
        for(int i=2;i<=lim;i++) inv[i] = mod(p - 1ll * (p / i) * inv[p%i] % p);
        for(int i=1;i<lim;i<<=1)
        {
            int wn = qp(g,(p - 1) / (i<<1));
            for(int j=0,ww=1;j<i;j++,ww=1ll*ww*wn%p) w[i+j] = ww;
        }
    }

    void ntt(vector<int> &f,int n,int op)
    {
        if(last!=n)
        {
            for(int i=1;i<n;i++) r[i] = (r[i>>1]>>1)|((i&1)?(n>>1):0);
            last=n;
        }
        for(int i=1;i<n;i++) if(i<r[i])swap(f[i],f[r[i]]);
        for(int i=1;i<n;i<<=1)
        for(int j=0;j<n;j+=i<<1)
            for(int k=0;k<i;k++)
            {
                int x=f[j+k],y=1ll*f[i+j+k]*w[i+k]%p;
                f[j+k]=mod(x+y);f[i+j+k]=mod(x-y+p);
            }
        if(op==-1)
        {
            reverse(&f[1],&f[n]);
            for(int i=0;i<n;i++)f[i]=1ll*f[i]*inv[n]%p;
        }
    }
}
using NTT::ntt;