NTT）

前言

众所周知，这两个东西都是用来算多项式乘法的。

对于这种常人思维难以理解的东西，就少些理解，多背板子吧！

~~因此只总结一下思路和代码，什么概念和推式子就靠巨佬们吧~~

推荐自为风月马前卒巨佬的概念和定理都非常到位的总结

推荐ppl巨佬的简明易懂的总结

FFT

多项式乘法的蹊径——点值表示法

一般我们把两个长度为 $n$ 的多项式乘起来，就类似于做竖式乘法，一位一位地乘再加到对应位上，是 $O(n^2)$ 的

如何优化？直接看是没有思路的，只好另辟蹊径了。

多项式除了我们常用的系数表示法 $y=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}$ 以外，还可以用点值表示法。

所谓点值表示法，就是相当于用函数图像上 $n$ 个点的坐标 $(x_i,y_i)$ 表示一个 $n$ 次多项式，显然这个表示是唯一的。

我们可以把系数表示转化成点值表示。点值表示下的多项式怎么相乘呢？就是选相同的 $x_i$ ，把对应的 $y_i$ 相乘。

当然，两个长度为 $n$ 的多项式相乘得到的是长度为 $2n-1$ 的多项式，需要 $2n-1$ 个点值才能唯一表示，因此一开始两个多项式也要选 $2n-1$ 个点表示。

举一个例子

$(x+1)(x+2)$ → → → → $x^2+3x+2$
↓ $\qquad\qquad\qquad\qquad\qquad\qquad$ ↑
↓（点值） $\qquad\qquad\qquad$ （系数）↑
↓ $\qquad\qquad\qquad\qquad\qquad\qquad$ ↑
$(0,1)(1,2)(2,3)$ （相乘） $\qquad\quad$ ↑
$(0,2)(1,3)(2,4)$ → → → → $(0,2)(1,6)(2,12)$
可是，随意选 $O(n)$ 个点，计算 $y$ 是 $O(n)$ 的，总时间还是 $O(n^2)$ 的，把它还原成系数表达式更不轻松。

所以，如果选的点很普通，这只是一条蹊径，并不是一条捷径。

神奇的单位根

介绍一个神奇的东西—— $n$ 次单位根（记为 $\omega$ ）。

它是 $n$ 个复数的集合，每一个的 $n$ 次方都等于 $1$ 。其中的一个是 $e^{{2\pi i}\over n}$ ，记为 $\omega_n$ 。

普及一下欧拉公式， $e^{\theta i}=\cos\theta+(\sin\theta)i$ ， $\theta$ 就是这个复数向量的旋转角。显然满足 $(\omega_n)^n=e^{2\pi i}=1$ ，那么它的 $k$ 次方 $(\omega_n^k)^n=e^{2k\pi i}$ 也等于 $1$ 。

于是可以看出，满足 $n$ 次方等于 $1$ 的复数的取值只会有 $n$ 个，为 $\omega_n^k(k\in[0,n-1])$ ，因为会有 $\omega_n^{n+k}=\omega_n^k$ 。

这 $n$ 个复数向量在单位圆上呈放射状。下面是算导的图片：

介绍消去引理 $ω^{dk}_{dn}=ω^k_n$ ，证明很容易的。

DFT&IDFT

单位根有什么用呢？

看看我们把 $\omega_n^k(k\in[0,n-1])$ 分别带入多项式求点值会发生什么就知道了。这个过程叫DFT。

假设 $n$ 为偶数，那么我们可以把它的奇偶项分开，用两个多项式表示它

$A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n−2}x^{{\frac n 2}−1}$

$A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n−1}x^{{\frac n 2}−1}$

那么 $A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$ 。

注意，下面的变化用到了 $ω^{2k}_n=ω^k_{\frac n 2}$ ， $ω^n_n=1$ 和 $ω^{\frac n 2}_n=-1$ 。

首先带入单位根

$A(ω^k_n)=A^{[0]}(ω^{2k}_n)+ω^k_nA^{[1]}(ω^{2k}_n)\\\qquad=A^{[0]}(ω^k_{\frac n 2})+ω^k_nA^{[1]}(ω^k_{\frac n 2})(k<\frac n 2)$

那 $k\geq\frac n 2$ 时又会发生什么呢？把它变成 $\frac n 2+k$

$A(ω^{\frac n 2+k}_n)=A^{[0]}(ω^{n+2k}_n)+ω^{\frac n 2+k}_nA^{[1]}(ω^{n+2k}_n)\\\qquad\quad=A^{[0]}(ω^{2k}_n)+ω^{\frac n 2}_nω^k_nA^{[1]}(ω^{2k}_n)\\\qquad\quad =A^{[0]}(ω^k_{\frac n 2})-ω^k_nA^{[1]}(ω^k_{\frac n 2})$

也就是说，如果我们要求一个长度为 $n$ 的多项式取 $\omega_n^k(k\in[0,n-1])$ 的 $n$ 个点值，我们只需要求出两个长度为 $n/2$ 的多项式取 $\omega_{\frac n 2}^k(k\in[0,\frac n 2-1])$ 的 $\frac n 2$ 个点值，再通过上述两个式子合并结果。这完全就是个递归过程，很容易写一个函数来计算。

可以由算法写出DFT的复杂度 $T(n)=2T(\frac n 2)+O(n)=O(n\log n)$

当然，别忘了还原成系数表示，这个过程叫做IDFT。

蒟蒻觉得这里理解太麻烦了，因此不再证明IDFT的过程，想了解的参见其它的总结。

只需要记住，IDFT的 $\omega_n$ 是 $e^{-\frac{2\pi i}n}$ ，最后的结果除以 $n$ ，其它过程同DFT，可以写在一个函数里。具体看下面的代码：

void FFT(R complex<double>*a,R int n,R int op){//op=1为DFT，-1为IDFT
    if(!n)return;//为了方便，n的意义与上面不一样，这里的n是a0、a1的长度
    complex<double>a0[n],a1[n];
    for(R int k=0;k<n;++k)
        a0[k]=a[k<<1],a1[k]=a[k<<1|1];//奇偶项分离
    FFT(a0,n>>1,op);FFT(a1,n>>1,op);//递归处理
    R complex<double>wn(cos(PI/n),sin(PI/n)*op),w(1,0);//单位根
    for(R int k=0;k<n;++k,w*=wn)//k从到n/2
        a[k]=a0[k]+w*a1[k],a[k+n]=a0[k]-w*a1[k];
}

递归版过程

引入例题：洛谷P3803 【模板】多项式乘法（FFT）

由于系数很小，我们不必担心精度的问题了（~~当然float是使不得的~~）

递归版代码：

#include<cstdio>
#include<cmath>
#include<complex>
#include<iostream>
#define R register
#define G c=getchar()
using namespace std;
const int N=4.2e6;
const double PI=acos(-1);//自定义π
complex<double>f[N],g[N];
inline int in(){
    R char G;
    while(c<'-')G;
    return c&15;
}
void FFT(R complex<double>*a,R int n,R int op){//同上
    if(!n)return;
    complex<double>a0[n],a1[n];
    for(R int i=0;i<n;++i)
        a0[i]=a[i<<1],a1[i]=a[i<<1|1];
    FFT(a0,n>>1,op);FFT(a1,n>>1,op);
    R complex<double>W(cos(PI/n),sin(PI/n)*op),w(1,0);
    for(R int i=0;i<n;++i,w*=W)
        a[i]=a0[i]+w*a1[i],a[i+n]=a0[i]-w*a1[i];
}
int main(){
    R int n,m,i;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i]=in();
    for(i=0;i<=m;++i)g[i]=in();
    for(m+=n,n=1;n<=m;n<<=1);//把长度补到2的幂，不必担心高次项的系数，因为默认为0
    FFT(f,n>>1,1);FFT(g,n>>1,1);//DFT
    for(i=0;i<n;++i)f[i]*=g[i];//点值直接乘
    FFT(f,n>>1,-1);//IDFT
    for(i=0;i<=m;++i)printf("%.0f ",fabs(f[i].real())/n);//注意结果除以n，小心“-0”
    puts("");
    return 0;
}

好记又好写的递归版结果怎样呢？

~~Fast Fast TLE!~~只有77分。

最主要的原因在于，空间调用太多了。

蝴蝶操作和Rader排序

针对效率太低的问题，我们继续观察FFT实现过程中的更多特点。

观察这一句代码a[k]=a0[k]+w*a1[k],a[k+n]=a0[k]-w*a1[k];，在操作过程中，取出了a0[k]和a1[k]的值，通过计算又求出了a[k]和a[k+n]的值。我们把这样的一次运算叫做“蝴蝶操作”。

这样的操作有什么特点呢？我们试着将a0和a1合并成一个数组a，每一次蝴蝶操作后，取出了a[k]和a[k+n]的值，又求出了a[k]和a[k+n]的值。最后，整个数组都完成了求值，而并没有用到两个数组！

以 $n=8$ 为例，看看递归过程的结构

其实，我们完全可以递推求解！假设 $a$ 数组已经变成了第四层的样子，先对a0和a4、a2和a6、a1和a5、a3和a7蝴蝶操作，变成第三层；再对a0和a2、a4和a6、a1和a3、a5和a7蝴蝶操作，变成第二层；最后对a0和a1、a2和a3、a4和a5、a6和a7蝴蝶操作，变成第一层，FFT就完成了，而空间复杂度仅为 $O(n)$ 。这个过程可以用循环来控制。

剩下的问题就是把初始的数组变成最后一层的样子了。先别急着写一个递归函数暴力把位置换过去。来观察一下最后序列的编号的二进制表示000,100,010,110,001,101,011,111，是不是与原来000,001,010,011,100,101,110,111相比，每个位置上的二进制位都反过来了？这样的变化叫做Rader排序。

我们可以 $O(n)$ 将Rader排序的映射关系求出。设 $i$ Rader排序后的数为 $r_i$ ，我们可以通过 $r_{\frac i 2}$ 递推求出，具体方法可以看看代码&动动脑筋qwq

#include<cmath>
#include<cstdio>
#define R register
#define I inline
using namespace std;
const int N=4.2e6;
const double PI=acos(-1);
int n,r[N];
struct C{//手写complex，比STL快一点点
	double r,i;
	I C(){r=i=0;}
	I C(R double x,R double y){r=x;i=y;}
	I C operator+(R C&x){return C(r+x.r,i+x.i);}
	I C operator-(R C&x){return C(r-x.r,i-x.i);}
	I C operator*(R C&x){return C(r*x.r-i*x.i,r*x.i+i*x.r);}
	I void operator+=(R C&x){r+=x.r;i+=x.i;}
	I void operator*=(R C&x){R double t=r;r=r*x.r-i*x.i;i=t*x.i+i*x.r;}
}f[N],g[N];
I int in(){
	R char c=getchar();
	while(c<'-')c=getchar();
	return c&15;
}
I void FFT(R C*a,R int op){
	R C W,w,t,*a0,*a1;
	R int i,j,k;
	for(i=0;i<n;++i)//根据映射关系交换元素
		if(i<r[i])t=a[i],a[i]=a[r[i]],a[r[i]]=t;
	for(i=1;i<n;i<<=1)//控制层数
		for(W=C(cos(PI/i),sin(PI/i)*op),j=0;j<n;j+=i<<1)//控制一层中的子问题个数
			for(w=C(1,0),a1=i+(a0=a+j),k=0;k<i;++k,++a0,++a1,w*=W)
				t=*a1*w,*a1=*a0-t,*a0+=t;//蝴蝶操作
}
int main(){
    R int m,i,l=0;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i].r=in();
    for(i=0;i<=m;++i)g[i].r=in();
    for(m+=n,n=1;n<=m;n<<=1,++l);
    for(i=0;i<n;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//递推求r
    FFT(f,1);FFT(g,1);
    for(i=0;i<n;++i)f[i]*=g[i];
    FFT(f,-1);
    for(i=0;i<=m;++i)printf("%.0f ",fabs(f[i].r)/n);
    puts("");
    return 0;
}

NTT

原根

FFT的全过程都是基于小数运算，一旦数据较大就会丢失精度。我们如何在整数范围内寻找和单位根 $\omega_n^n=1$ 性质一样的数呢？

显然，如果不在模意义下这样的式子是不成立的。

引入原根，设为 $g$ ，如果 $g^n=1(\mod p)$ 中 $n$ 的最小正整数解为欧拉函数 $\phi(p)$ ，则称 $g$ 为 $p$ 的一个原根。

为了理解NTT，了解定义就够了，当然更多内容可以参考YL的总结

用 $\omega_n$ 替换 $g$ ，根据定义显然也满足消去引理 $ω^{dk}_{dn}=ω^k_n$ 。那么放到多项式乘法里没问题。

一个问题在于，递推过程中如何根据求出不同子问题的 $ω$ ？显然因为 $g^{p-1}=\omega_n^n=1$ ，有 $\omega_n=g^{\frac{p-1}n}$

当然，多项式的长度都被补成了 $2$ 的幂，这就要求 $\phi(p)$ 的质因子中含有大量的 $2$ ；另外，IDFT中的单位根要对原根求逆，最后除以 $n$ 也需要求逆，为了方便还要求 $p$ 是质数。这样，能被用作NTT的模数的数非常少，常见的就是 $998244353(998244352=2^{23}×7×17)$ ， $3$ 是它的一个原根。

上面那题的NTT代码

#include<cstdio>
#define I inline
#define RG register
#define RL RG LL
#define R RG int
typedef long long LL;
const int N=4.2e6,YL=998244353;//一起来%YL
int n,f[N],g[N],r[N];
I int in(){
    RG char c=getchar();
    while(c<'-')c=getchar();
    return c&15;
}
I int qpow(RL b,R k){//快速幂
    RL a=1;
    for(;k;k>>=1,b=b*b%YL)
        if(k&1)a=a*b%YL;
    return a;
}
I void NTT(R*a,R op){
    R i,j,k,d,wn,w,t,*a0,*a1;
    for(i=0;i<n;++i)
        if(i<r[i])t=a[i],a[i]=a[r[i]],a[r[i]]=t;
    for(i=1;i<n;i<<=1){
        wn=qpow(3,(YL-1)/(d=i<<1));//原根变换
        if(op&2)wn=qpow(wn,YL-2);//注意要求逆
        for(j=0;j<n;j+=d)
            for(w=1,a1=(a0=a+j)+i,k=0;k<i;++k,++a0,++a1,w=(LL)w*wn%YL)
                t=(LL)*a1*w%YL,*a1=(*a0-t+YL)%YL,*a0=(*a0+t)%YL;
    }
}
int main(){
    R m,i,l=0;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i]=in();
    for(i=0;i<=m;++i)g[i]=in();
    for(m+=n,n=1;n<=m;n<<=1,++l);
    for(i=0;i<n;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    NTT(f,1);NTT(g,1);
    for(i=0;i<n;++i)f[i]=(LL)f[i]*g[i]%YL;
    NTT(f,-1);
    R inv=qpow(n,YL-2);//同样注意要求逆
    for(i=0;i<=m;++i)printf("%lld ",(LL)f[i]*inv%YL);
    puts("");
    return 0;
}

posted @ 2018-07-17 15:17 Flash_Hu 阅读(851) 评论(3) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

公告

点击按钮来显示日期。

欢迎各路大佬和我QQ或邮件交流，也欢迎查错，我会尽量及时回复，感谢！

本作品由胡锦添采用知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议进行许可。

昵称： Flash_Hu
园龄： 7年2个月
粉丝： 309
关注： 80

+加关注

2025年3月

日

一

二

三

四

五

六

FlashHu

——蝌

FFT/NTT总结+洛谷P3803 【模板】多项式乘法（FFT）（FFT/NTT）

前言

FFT

多项式乘法的蹊径——点值表示法

神奇的单位根

DFT&IDFT

递归版过程

蝴蝶操作和Rader排序

NTT

原根

公告

搜索

积分与排名

随笔分类 (432)

随笔档案 (135)

被同级巨佬全方位吊打

长江后浪推前浪的巨佬学弟学妹们

一走出校门又被外校巨佬吊打

自学成才，爆肝裕民

阅读排行榜

评论排行榜

推荐排行榜

最新评论