NTT快速数论变换

NTT快速数论变换

some 多项式变换

DFT:离散傅里叶变换 $\longrightarrow $ \(O(n^2)\)计算多项式乘法
FFT:快速傅里叶变换 $\longrightarrow $ \(O(n\log n)\)计算多项式乘法
FNTT/NTT:快速傅里叶变换的优化版 $\longrightarrow $ 优化常数及误差
FWT:快速沃尔什变换 $\longrightarrow $ 利用类似FFT的东西解决一类卷积问题
MTT:毛爷爷的FFT $\longrightarrow $ 非常nb/任意模数
FMT 快速莫比乌斯变化

定义

数论变换(\(NTT\))是通过将离散傅立叶变换化为 \(F={\mathbb {Z}/p}\),整数模质数 \(p\)

这是一个 有限域,只要 \(n\) 可除 \(p-1\),就存在本原 \(n\) 次方根,所以对于正整数 \(\xi\),我们有 \(p=\xi n+1\)

具体来说,对于质数 \(p=qn+1, (n=2^m)\),原根 \(g\) 满足 \(g^{qn} \equiv 1 \pmod p\), 将 \(g_n=g^q\pmod p\) 看做 \(\omega_n\) 的等价,则其满足相似的性质,比如 \( g_n^n \equiv 1 \pmod p, g_n^{n/2} \equiv -1 \pmod p\)

因为这里涉及到数论变化,所以 \(N\)(为了区分 FFT 中的 \(n\),我们把这里的 \(n\) 称为 \(N\))可以比 \(FFT\) 中的 \(n\) 大,但是只要把 \(\frac{qN}{n}\) 看做这里的 \(q\) 就行了,能够避免大小问题。

就是 \(g^{qn}\) 的等价 \(e^{2\pi n}\)

迭代到长度 \(l\)\(g_l = g^{\frac{p-1}{l}}\),或者 \(\omega_n = g_l = g_N^{\frac{N}{l}} = g_N^{\frac{p-1}{l}}\)

实现

快速数论变换(\(FNTT\))是数论变换(\(NTT\))增加分治操作之后的快速算法。
快速数论变换使用的分治办法,与快速傅里叶变换使用的分治办法完全一致。

\(FFT\)中,我们需要用到复数,复数虽然很神奇,但是它也有自己的局限性——需要用\(double\)类型计算,精度太低
那有没有什么东西能够代替复数且解决精度问题呢?这个东西,叫原根

原根为什么能代替单位根进行运算

因为它具有和单位根相同的性质
在FFT中,我们用到了单位根的四条性质,而原根也满足这四条性质
这样我们最终可以得到一个结论

\[\omega_n\equiv g^{\frac{p-1}{n}}\pmod p \]

然后把FFT中的\(\omega_n\)都替换掉就好了

code

il void ntt(int q[],int lim,int k){
    for(ri int i=0;i<lim;++i){
        if(i<r[i]) swap(q[i],q[r[i]]);
    }
    for(ri int mid=1,wn;mid<lim;mid<<=1){
        wn=qpow(k?G:I,(P-1)/(mid<<1));
        for(ri int R=mid<<1,j=0,w=1,x,y;j<lim;j+=R,w=1){
            for(ri int t=0;t<mid;++t,w=w*wn%P){
                x=q[j+t],y=w*q[j+mid+t]%P;
                q[j+t]=(x+y+P)%P,q[j+mid+t]=(x-y+P)%P;
            }
        }
    }
    return;
}
code
#include<bits/stdc++.h>
#define il inline
#define cs const
#define ri register
#define int long long
using namespace std;

namespace Q{
    il int rd(){
        ri int x=0;ri bool f=0;ri char c=getchar();
        while(!isdigit(c)) f|=(c==45),c=getchar();
        while(isdigit(c)) x=x*10+(c^48),c=getchar();
        return f?-x:x;
    }
    il void wt(int x){
        if(x<0) x=-x,putchar(45);
        if(x>=10) wt(x/10);
        return putchar(x%10|48),void();
    }
} using namespace Q;

cs int N=2097152,P=998244353,G=3,I=332748118;

namespace NTT{
    int l,r[N],f[N],g[N];
    il int qpow(int a,int b,int p=P){
        ri int as=1;
        while(b>0){
            if(b&1) as=as*a%p;
            a=a*a%p,b>>=1;
        }
        return as;
    }
    il int inv(int a,int p=P){
        return qpow(a,p-2,p);
    }
    il void ntt(int q[],int lim,int k){
        for(ri int i=0;i<lim;++i){
            if(i<r[i]) swap(q[i],q[r[i]]);
        }
        for(ri int mid=1,wn;mid<lim;mid<<=1){
            wn=qpow(k?G:I,(P-1)/(mid<<1));
            for(ri int R=mid<<1,j=0,w=1,x,y;j<lim;j+=R,w=1){
                for(ri int t=0;t<mid;++t,w=w*wn%P){
                    x=q[j+t],y=w*q[j+mid+t]%P;
                    q[j+t]=(x+y+P)%P,q[j+mid+t]=(x-y+P)%P;
                }
            }
        }
        return;
    }
} using namespace NTT;

signed main(){
    int n=rd(),m=rd(),lim=1;
    while(lim<=n+m) lim<<=1,l++;
    for(ri int i=0;i<=n;++i) f[i]=rd();
    for(ri int i=0;i<=m;++i) g[i]=rd();
    for(ri int i=0;i<lim;++i){
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    }
    ntt(f,lim,1),ntt(g,lim,1);
    for(ri int i=0;i<lim;++i) f[i]=f[i]*g[i]%P;
    ntt(f,lim,0);
    for(ri int i=0,s=inv(lim);i<=n+m;++i){
        wt(f[i]*s%P),putchar(32);
    }
    return 0;
}

edit

posted @ 2023-02-15 22:03  雨夜风月  阅读(152)  评论(0编辑  收藏  举报