博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
[多项式算法](Part 2)NTT 快速数论变换 学习笔记
阅读量:5109 次
发布时间:2019-06-13

本文共 3387 字,大约阅读时间需要 11 分钟。

其他多项式算法传送门:

\(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)\),也可以查表:


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

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

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

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

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

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

代码:

#include 
#include
#include
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
j)std::swap(Pol[i],Pol[j]); for(int l=n>>1;(j^=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
=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

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

参考资料:(\(STO\ Dalao\)


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

你可能感兴趣的文章
最近看NCZ的JS高级程序设计整理的一些代码
查看>>
浙江省第十二届省赛 Beauty of Array(思维题)
查看>>
NOIP2013 提高组 Day1
查看>>
UVA 1602 Lattice Animals
查看>>
bzoj千题计划219:bzoj1568: [JSOI2008]Blue Mary开公司
查看>>
[笔记]STM32使用非8M晶振时如何修改代码
查看>>
个人对vue生命周期的理解
查看>>
cocos2dx 3.x simpleAudioEngine 长音效被众多短音效打断问题
查看>>
Section 1.2 dualpal
查看>>
存储(硬件方面的一些基本术语)
查看>>
Dithering-视觉的奇特现象
查看>>
观察者模式
查看>>
转】MyEclipse使用总结——MyEclipse文件查找技巧
查看>>
Weka中数据挖掘与机器学习系列之基本概念(三)
查看>>
Java-文件上传和下载
查看>>
Memory and Trident(CodeForces 712B)
查看>>
Win磁盘MBR转换为GUID
查看>>
大家在做.NET B/S项目的时候多用什么设技术啊?
查看>>
投资策略 ——摘自凤凰网
查看>>
Java SE和Java EE应用的性能调优
查看>>