[多项式算法](Part 2)NTT 快速数论变换学习笔记

最新推荐文章于 2024-09-04 16:36:43 发布

weixin_30337157

最新推荐文章于 2024-09-04 16:36:43 发布

阅读量174

点赞数

原文链接：http://www.cnblogs.com/LanrTabe/p/11305607.html

版权

其他多项式算法传送门：

[多项式算法](Part 1)FFT 快速傅里叶变换学习笔记

[多项式算法](Part 3)MTT 任意模数FFT/NTT 学习笔记

[多项式算法](Part 4)FWT 快速沃尔什变换学习笔记

[多项式算法](Part 5)分治FFT 学习笔记

\(2.Medium-NTT(FNT)\)

定义

NTT\((Number\ Theoretic\ Transforms)\)

(也称为\(Fast\ Number-Theoretic\ Transform\)，简称\(FNT\))

中文名称：快速数论变换

~~(Not True Transforms)~~

\(Q:FFT\)懂了，\(NTT\)又有什么用呢？\(FFT\)已经够用了啊？？

\(A:\)别急，相对于\(FFT\)来说，\(NTT\)具有更强的针对性，\(NTT\)可以在取模意义下进行多项式乘法计算，从而避免\(FFT\)的\(double\)失精的情况，但\(NTT\)对模数有特殊要求，这在下面会提到。

那么\(NTT\)是怎么实现的呢？

其实\(NTT\)的实现方法与\(FFT\)几乎无异，只需把复数运算换成取模运算即可。

现在来思考为什么\(FFT\)的点值表示要用单位根呢？

因为单位根有如下性质供\(FFT\)利用以加速：

\(1.\)

\[\omega_n^n=1\]

这一点在\(DFT\)时需要用到

\(2.\)

所有单位根互不相同，这样才能算出正确答案（例如\(n\)个\(n\)元方程组只有线性无关才有唯一解）。

\(3.\)

\[\begin{cases} \omega_{2n}^{2k}=\omega_n^k(k<\frac n2)\\ \omega_n^{k+\frac n2}=-\omega_n^k(k\ge\frac n2) \end{cases}\]

显然这样分治才能继续进行

\(4.\)

\[\sum_{k=0}^{n-1}\omega_n^{k(j-i)}= \begin{cases} 0(i\not=j)\\ n(i=j) \end{cases}\]

这一点\(IDFT\)时有需要。

所以说接下来我们需要找到合适的数满足上面几条性质来代替单位根。

原根

设有质数\(p\)，则\(p\)的原根\(g\)满足\(g^k\mod p(1\le k<p-1)\)互不相同。

那么若质数\(p=q*n+1(n=2^x)\)，则根据费马小定理有\(a^{p-1}\equiv 1(\mod p)\)

显然，\(g_0\equiv g_{p-1}\equiv 1(\mod p)\)

若设\(\omega_n=g^q\)，则可以得到\(n\)个不相同的数：\(\omega_n^k(0\le k<n)\)

满足性质2

并且\(\omega_n^n\equiv g^{p-1}\equiv 1(\mod p)\)

满足性质1

那么就可以得到：

\(\because p=q*n+1=\frac q2*2n+1,\omega_n=g^q\)

\(\therefore \omega_{2n}=g^{\frac q2}\)

\(\omega_{2n}^{2k}=g^{\frac q2*2k}=g^{qk}=\omega_n^k\)

且

\(\omega_n^{k+\frac n2}\)

\(=\omega_n^k*\omega_n^{\frac n2}\)

\(\because (\omega_n^{\frac n2})^2=\omega_n^n=1\)

\(\therefore \omega_n^{\frac n2}=\pm1\)

\(\because \omega\)互不相同

\(\therefore \omega_n^{\frac n2}=-1\)

\(\therefore \omega_n^{k+\frac n2}=-\omega_n^k\)

满足性质3

最后：对于\(\sum_{k=0}^{n-1}\omega_n^{k(j-i)}\)

若\(i=j\)，则：

\(\sum_{k=0}^{n-1}\omega_n^{k(j-i)}\)

\(=n*\omega_n^0\)

\(=n\)

若\(i\not=j\)，则有：

\(\sum_{k=0}^{n-1}\omega_n^{k(j-i)}\)

\(=\frac{(q^n-1)a_0}{q-1}\)(等比数列求和，此处\(q\)为公比\((\omega_n^{j-i})\)，\(a_0\)为首项\((\omega_n^0)\))

\(\because q^n-1=\omega_n^{n(j-i)}-1=(\omega_n^n)^{j-i}-1=0\)

\(\therefore \sum_{k=0}^{n-1}\omega_n^{k(j-i)}=0\)

综上所述，满足性质4

\(Q:\)什么？原根怎么求？

\(A:\)~~我怎么知道，百度啊~~

一般情况下只需要记住\(998244353\)的原根是\(3\)就好\((998244353=119*2^{23}+1)\)，也可以查表：[Miskcoo's Space]FFT用到的各种素数

于是，照着上面说的，\(NTT\)的框架就出来了：

把复数运算换成取模即可。

\(Q:n\)不满\(2^{23}\)怎么办？补满太慢。

\(A:\)把\(\omega_n\)换成\(g^{\frac{p-1}{n}}\)即可。

P3803 【模板】多项式乘法（FFT）

\(Q:\)为什么是\(FFT\)模板？

\(A:\)~~找不到NTT的~~ 因为这题答案不会超过\(998244353\)，那么用来取模就并不会产生影响。

代码：

#include <cstdio>
#include <cctype>
#include <algorithm>

char In[1000005],*p1=In,*p2=In;
#define Getchar (p1==p2&&(p2=(p1=In)+fread(In,1,1000000,stdin),p1==p2)?EOF:*p1++)
inline int Getint()
{
    register int x=0,c;
    while(!isdigit(c=Getchar));
    for(;isdigit(c);c=Getchar)x=x*10+(c^48);
    return x;
}

const int p=998244353,g=3;

int Pow(int a,int b)
{
    if(b<0)return Pow(Pow(a,p-2),-b);//逆元
    int Res=1;
    for(;b;b>>=1)
    {
        if(b&1)Res=Res*1LL*a%p;
        a=a*1LL*a%p;
    }
    return Res%p;
}

int n,m,Maxl;
int a[2100005],b[2100005];

void NTT(int Pol[],int op)//op=1为DFT,-1为IDFT
{
    for(int i=0,j=0;i<n;++i)
    {
        if(i>j)std::swap(Pol[i],Pol[j]);
        for(int l=n>>1;(j^=l)<l;l>>=1);
    }
    for(register int i=2;i<=n;i<<=1)
    {
        int m=i>>1,Rs=Pow(g,op*(p-1)/i);
        for(register int j=0;j<n;j+=i)
        {
            int Root=1;
            for(register int k=0;k<m;++k)
            {
                int Tmp=Root*1LL*Pol[j+k+m]%p;
                Root=Root*1LL*Rs%p;
                //int Root=Pow(Omega,op*(p-1)/n*k)
                //不预处理单位根，使用秦九韶算法加速
                Pol[j+k+m]=Pol[j+k]-Tmp;
                Pol[j+k]=Pol[j+k]+Tmp;
                Pol[j+k+m]<0?Pol[j+k+m]+=p:0;//减少取模
                Pol[j+k]>=p?Pol[j+k]-=p:0;
            }
        }
    }
}

int main()
{
    n=Getint(),m=Getint();
    for(register int i=0;i<=n;++i)a[i]=Getint();
    for(register int i=0;i<=m;++i)b[i]=Getint();
    for(Maxl=n+m,n=2;n<=Maxl;n<<=1);
    NTT(a,1),NTT(b,1);
    for(int i=0;i<n;++i)a[i]=a[i]*1LL*b[i]%p;
    NTT(a,-1);
    int Invn=Pow(n,p-2);
    for(int i=0;i<=Maxl;++i)printf("%d%c",int(a[i]*1LL*Invn%p),i==n-1?'\n':' ');
    //答案同样/n
    return 0;
}

相信有了\(FFT\)的基础，\(NTT\)应该很简单。

参考资料：（\(STO\ Dalao\)）

从多项式乘法到快速傅里叶变换-Miskcoo

从傅里叶变换到数论变换-Menci

转载于:https://www.cnblogs.com/LanrTabe/p/11305607.html