07-29 03:31 阅读 226

关注

FFT & NTT(入门)

FFT & NTT(入门)

写此篇只为加深印象。

1.FFT(快速傅里叶变换)

参考

对于求解F(n+m)=g(n)∗h(m)F(n+m)=g(n)∗h(m)的系数该怎么办?

首先要了解如何能表示一个多项式?

一系数表示法:这是我们通常在用的方法,比如一个二次函数f(x)=x2+2x+1f(x)=x2+2x+1
二点值表示法:对于一个k次多项式,我们可以用k+1个点确定这个多项式,比如f(x)=x2+2x+1f(x)=x2+2x+1,我们可以用(0,1)(−1,0)(−2.1)(0,1)(−1,0)(−2.1)这三个点确定这个二次多项式

因此我们有两种方法

一:直接将对应系数相乘
二:选择n+m+1n+m+1个点,将这些点值全都求出来,最后确定这个多项式
复杂度都为O(n2)O(n2)

FFT就是第二种方法的优化

FFT的核心是什么?
(引用神仙学长chasedeath的话)
1.单位根
2.多项式与点值式的转化
3.分治思想

(1).首先介绍单位根

定义:数学上，n次单位根是n次幂为1的复数。它们位于复平面的单位圆上，构成正n边形的顶点，其中一个顶点是1。(by baidu)
令ωkn=cos2kπn+isin2kπnωnk=cos⁡2kπn+isin⁡2kπn

性质一:ω0n,ω1n...,ωn−1nωn0,ωn1...,ωnn−1互不相同

性质二:ωkn=ω2k2nωnk=ω2n2k(代入化简一下即可)

性质三:ωk+n/2n=−ωknωnk+n/2=−ωnk(代入化简一下即可)

性质四:(ωkn)2=ω2kn(ωnk)2=ωn2k(复数相乘,辐角相加)

性质五:∑n−1i=0ωkin=0&k=/0∑i=0n−1ωnki=0&k=⧸0(等比数列求和)

(2).快速傅里叶变换(将系数转化为点值)

我们可以将单位根的0到n-1次幂代入其中,求出答案
我们设g(x)=∑n−1i=0aixig(x)=∑i=0n−1aixi将这些点代入,如果硬算复杂度为O(n2)O(n2)

分治开始登场

我们按照下标的奇偶将系数分为两类(不要问我为什么,我也不知道)
g1(x)=a1+a3x1+a5x2+...+a(n−1)xn/2−1g1(x)=a1+a3x1+a5x2+...+a(n−1)xn/2−1(默认n为2的幂)
g2(x)=a0+a2x1+a4x2+...+a(n−2)xn/2−1g2(x)=a0+a2x1+a4x2+...+a(n−2)xn/2−1
g(x)=x∗g1(x2)+g2(x2)g(x)=x∗g1(x2)+g2(x2)
这时我们将ωknωnk和ωk+n/2nωnk+n/2代入式子得:
g(ωkn)=ωkng1((ωkn)2)+g2((ωkn)2)=ωkng1(ω2kn)+g2(ω2kn)g(ωnk)=ωnkg1((ωnk)2)+g2((ωnk)2)=ωnkg1(ωn2k)+g2(ωn2k)`
g(ωk+n/2n)=ωk+n/2ng1((ωk+n/2n)2)+g2((ωk+n/2n)2)=ωk+n/2ng1(ω2k+nn)+g2(ω2k+nn)=−ωkng1(ω2kn)+g2(ω2kn)g(ωnk+n/2)=ωnk+n/2g1((ωnk+n/2)2)+g2((ωnk+n/2)2)=ωnk+n/2g1(ωn2k+n)+g2(ωn2k+n)=−ωnkg1(ωn2k)+g2(ωn2k)
这时我们有一个惊人的发现只有符号发生了变化(oh,大师我悟了)
于是可以做到递归求解,时间就变成T(n)=T(n/2)+T(n/2)+O(n)T(n)=T(n/2)+T(n/2)+O(n)

(3)快速傅里叶逆变换(将点值重新转化为系数)

将(y0,y1,y2,...,yn−1)(y0,y1,y2,...,yn−1)当作系数,做一次傅里叶变换
即ak=1n∑n−1i=0yi(ω−kn)iak=1n∑i=0n−1yi(ωn−k)i

证明:

我们有式子:yk=∑n−1i=0ai(ωkn)iyk=∑i=0n−1ai(ωnk)i
将其代入得:ak=1n∑n−1i=0∑n−1j=0ajωi(j−k)nak=1n∑i=0n−1∑j=0n−1ajωni(j−k)
当j=kj=k时,右边Ans=akAns=ak
当j=k̸j=k̸时,右边对于每个jj,Ans=aj∗∑n−1s=0ω(j−k)n=0Ans=aj∗∑s=0n−1ωn(j−k)=0(由性质5易得)

(4)优化FFT

优化一:蝴蝶效应(B格很高)

就是复数的乘法据说比较慢,我们不要计算两次.

优化二:迭代反转(这个就很有用)

众所周知,递归跑地巨慢无比,我们可以用循环实现,然后一直向上合并.

发现位置跟这个数的二进制有关,于是将它们最终位置固定住即可

代码

#include<bits/stdc++.h>using namespace std;bool f1;#define Max(a,b) ((a)<(b))&&((a)=(b))#define Min(a,b) ((a)>(b))&&((a)=(b))#define LL long long#define db doubleinline LL rd() {
	LL res=0;bool f=0;	char ch;	while(ch=getchar(),ch<48||ch>57)if(ch=='-')f=true;	do res=(res<<1)+(res<<3)+(ch^48);	while(ch=getchar(),ch>47&&ch<58);	return f?-res:res;
}const int M=4e6+5;const db Pi=acos(1);struct com {
	db x,y;//x+iy
	com operator +(const com &_) const {return (com)<%x+_.x,y+_.y%>;}
	com operator -(const com &_) const {return (com)<%x-_.x,y-_.y%>;}
	com operator *(const com &_) const {return (com)<%x*_.x-y*_.y,y*_.x+x*_.y%>;}
}A[M],B[M];int n,m,Mx=1,cnt;int pos[M];void FFT(com *a,int opt) {	for(int i=0; i<Mx; ++i)if(i<pos[i])swap(a[i],a[pos[i]]);	for(int len=1; len<Mx; len<<=1) {
		com now=(com)<%cos(Pi/len),opt*sin(Pi/len)%>;		for(int L=0; L<Mx; L+=len<<1) {
			com sw=(com)<%1,0%>;			for(int R=L; R<L+len; ++R,sw=sw*now) {
				com xx=a[R],yy=a[R+len]*sw;
				a[R]=xx+yy;
				a[R+len]=xx-yy;
			}
		}
	}
}bool f2;int main(){//	printf("%lf\n",(&f2-&f1)/1024.0/1024.0);//	freopen("a.in","r",stdin);//	freopen("a.ans","w",stdout);
	n=rd(),m=rd();	for(int i=0; i<=n; ++i)A[i]=(com)<%rd(),0%>;	for(int i=0; i<=m; ++i)B[i]=(com)<%rd(),0%>;	while(Mx<=n+m)Mx<<=1,cnt++;	for(int i=0; i<Mx; ++i)pos[i]=(pos[i>>1]>>1)|((i&1)<<cnt-1);	// 在原序列中 i 与 i/2 的关系是 ： i可以看做是i/2的二进制上的每一位左移一位得来 (by attack)
    // 那么在反转后的数组中就需要右移一位，同时特殊处理一下复数 
	FFT(A,1),FFT(B,1);	for(int i=0; i<Mx; ++i)A[i]=A[i]*B[i];
	FFT(A,-1);	for(int i=0; i<=n+m; ++i)printf("%d ",(int)(A[i].x/Mx+0.5));	return 0;
}

2.NTT(快速数论变换)

参考
由于FFT是利用复数计算的,导致精度可能会出现问题,于是就出现了NTT.
NTT是用原根来代替单位根,其它与FFT没什么区别.

注意:模数需为质数

(1)原根

定义:ar≡1(modP)ar≡1(modP)若最小的r=ϕ(m)r=ϕ(m)则称aa为模PP的原根

令P=qn+1P=qn+1(n为2的幂次),定义g为原根,则ωn=gqωn=gq,集合为(1,gq,g2q,...,g(n−1)q)(1,gq,g2q,...,g(n−1)q)
既然被用来替代单位根,那么其就具备单位根应该有的性质(性质对应上面的FFT)

性质一:根据原根的定义可知,gϕ(P)≡1(modP)gϕ(P)≡1(modP), ϕ(P)ϕ(P)为满足条件中最小的数,因此这些数两两不同

性质二:ω2n=gq/2(P=q/2∗2n+1)ω2n=gq/2(P=q/2∗2n+1)因此得证

性质三:ωk+n/2nωkn=ωn/2nωnk+n/2ωnk=ωnn/2,(ωn/2n)2≡1(modP)(ωnn/2)2≡1(modP)因为ωn/2n!=1ωnn/2!=1所以ωn/2n=−1ωnn/2=−1得证

性质四:易证

性质五:同理用等比数列求和即可

(2)找原根

那么我们该怎么寻找原根,是不是所有数都有原根?
首先寻找一个质数的原根,我们可以用枚举法来做

有个重要的小结论:若ar≡1(modP)ar≡1(modP)最小的r一定是P−1P−1的约数,因为由费马小定理可得aP−1≡1(modP)aP−1≡1(modP)恒成立(也可利用反证法)

10000000071000000007的原根是55
998244353998244353的原根是33
73400337340033的原根是33
57671695767169的原根是33
其它需要的话用此程序(跑得比较慢,随便扒下来的)

#include<bits/stdc++.h>using namespace std;#define LL long longusing namespace std;int mod, phi;int main() {    scanf("%d", &mod), phi=mod-1;    for(int i = 2, p = 1; i <= mod; ++i, p = 1) {    	for(int a = i; a != 1; (a *= i) %= mod, ++p) ;    	if(p == phi) {printf("%d\n", i);return 0;}
    }
}

代码

#include<bits/stdc++.h>using namespace std;bool f1;#define Max(a,b) ((a)<(b))&&((a)=(b))#define Min(a,b) ((a)>(b))&&((a)=(b))#define LL long longinline LL rd() {
	LL res=0;bool f=0;	char ch;	while(ch=getchar(),ch<48||ch>57)if(ch=='-')f=true;	do res=(res<<1)+(res<<3)+(ch^48);	while(ch=getchar(),ch>47&&ch<58);	return f?-res:res;
}const int P=998244353,G=3,G_i=332748118;const int M=1<<21;int n,m;int A[M],B[M];int Mx=1,cnt,now,sw=1;int fpow(int a,int cnt) {	int res=1;	for(int i=cnt; i; i>>=1,a=1LL*a*a%P)if(i&1)res=1LL*res*a%P;	return res;
}int pos[M];void NTT(int *a,int f) {	for(int i=0; i<Mx; ++i)if(i<pos[i])swap(a[pos[i]],a[i]);	for(int len=1; len<Mx; len<<=1) {
		now=fpow(f?G:G_i,(P-1)/(len<<1));		for(int L=0; L<Mx; sw=1,L+=len<<1) {			for(int R=L; R<L+len; ++R,sw=1LL*sw*now%P) {				int xx=a[R],yy=1LL*a[R+len]*sw%P;
				a[R]=xx+yy;(a[R]>=P)&&(a[R]-=P);
				a[R+len]=xx-yy;(a[R+len]<0)&&(a[R+len]+=P);
			}
		}
	}
}bool f2;int main(){//	printf("%lf\n",(&f2-&f1)/1024.0/1024.0);//	freopen("a.in","r",stdin);//	freopen("a.ans","w",stdout);
	n=rd(),m=rd();	for(int i=0; i<=n; ++i)A[i]=rd();	for(int i=0; i<=m; ++i)B[i]=rd();	while(Mx<=n+m)Mx<<=1,cnt++;	for(int i=0; i<Mx; ++i)pos[i]=(pos[i>>1]>>1)|((i&1)<<cnt-1);
	NTT(A,1),NTT(B,1);	for(int i=0; i<Mx; ++i)A[i]=1LL*A[i]*B[i]%P;
	NTT(A,0);int inv_Mx=fpow(Mx,P-2);	for(int i=0; i<=n+m; ++i)printf("%d ",1LL*A[i]*inv_Mx%P);	return 0;
}

来源https://www.cnblogs.com/syhyyy/p/15071226.html