FFT & NTT(入门)
写此篇只为加深印象。
1.FFT(快速傅里叶变换)
对于求解\(F(n+m)=g(n)*h(m)\)的系数该怎么办?
首先要了解如何能表示一个多项式?
一系数表示法:这是我们通常在用的方法,比如一个二次函数\(f(x)=x^2+2x+1\)
二点值表示法:对于一个k次多项式,我们可以用k+1个点确定这个多项式,比如\(f(x)=x^2+2x+1\),我们可以用\((0,1) (-1,0) (-2.1)\)这三个点确定这个二次多项式
因此我们有两种方法
一:直接将对应系数相乘
二:选择\(n+m+1\)个点,将这些点值全都求出来,最后确定这个多项式
复杂度都为\(O(n^2)\)
FFT就是第二种方法的优化
FFT的核心是什么?
(引用神仙学长chasedeath的话)
1.单位根
2.多项式与点值式的转化
3.分治思想
(1).首先介绍单位根
定义:数学上，n次单位根是n次幂为1的复数。它们位于复平面的单位圆上，构成正n边形的顶点，其中一个顶点是1。(by baidu)
令\(\omega_{n}^{k}=\cos \frac{2k\pi}{n}+i\sin \frac{2k\pi}{n}\)
性质一:\(\omega_{n}^{0},\omega_{n}^{1}...,\omega_{n}^{n-1}\)互不相同
性质二:\(\omega_{n}^{k}=\omega_{2n}^{2k}\)(代入化简一下即可)
性质三:\(\omega_{n}^{k+n/2}=-\omega_{n}^{k}\)(代入化简一下即可)
性质四:\({(\omega_{n}^{k})}^2=\omega_{n}^{2k}\)(复数相乘,辐角相加)
性质五:\(\sum_{i=0}^{n-1}\omega_{n}^{ki}=0 \& k=\not 0\)(等比数列求和)
(2).快速傅里叶变换(将系数转化为点值)
我们可以将单位根的0到n-1次幂代入其中,求出答案
我们设\(g(x)=\sum_{i=0}^{n-1}a_ix^{i}\)将这些点代入,如果硬算复杂度为\(O(n^2)\)
分治开始登场
我们按照下标的奇偶将系数分为两类(不要问我为什么,我也不知道)
\(g_1(x)=a_1+a_3x^1+a_5x^2+...+a_{(n-1)}x^{n/2-1}\)(默认n为2的幂)
\(g_2(x)=a_0+a_2x^1+a_4x^2+...+a_{(n-2)}x^{n/2-1}\)
\(g(x)=x*g_1(x^2)+g_2(x^2)\)
这时我们将\(\omega_{n}^{k}\)和\(\omega_{n}^{k+n/2}\)代入式子得:
\(g(\omega_{n}^{k})=\omega_{n}^{k}g_1({(\omega_{n}^{k})}^2)+g_2({(\omega_{n}^{k})}^2)
=\omega_{n}^{k}g_1(\omega_{n}^{2k})+g_2(\omega_{n}^{2k})\)`
\(g(\omega_{n}^{k+n/2})=\omega_{n}^{k+n/2}g_1({(\omega_{n}^{k+n/2})}^2)+g_2({(\omega_{n}^{k+n/2})}^2)
=\omega_{n}^{k+n/2}g_1(\omega_{n}^{2k+n})+g_2(\omega_{n}^{2k+n})
=-\omega_{n}^{k}g_1(\omega_{n}^{2k})+g_2(\omega_{n}^{2k})
\)
这时我们有一个惊人的发现只有符号发生了变化(oh,大师我悟了)
于是可以做到递归求解,时间就变成\(T(n)=T(n/2)+T(n/2)+O(n)\)
(3)快速傅里叶逆变换(将点值重新转化为系数)
将\((y_0,y_1,y_2,...,y_n-1)\)当作系数,做一次傅里叶变换
即\(a_k=\frac{1}{n}\sum_{i=0}^{n-1}y_i(\omega_{n}^{-k})^i\)
证明:
我们有式子:\(y_k=\sum_{i=0}^{n-1}a_i(\omega_{n}^{k})^i\)
将其代入得:\(a_k=\frac{1}{n}\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\omega_{n}^{i(j-k)}\)
当\(j=k\)时,右边\(Ans=a_k\)
当\(j=\not k\)时,右边对于每个\(j\),\(Ans=a_j*\sum_{s=0}^{n-1}\omega_{n}^{(j-k)}=0\)(由性质5易得)
(4)优化FFT
优化一:蝴蝶效应(B格很高)
就是复数的乘法据说比较慢,我们不要计算两次.
优化二:迭代反转(这个就很有用)
众所周知,递归跑地巨慢无比,我们可以用循环实现,然后一直向上合并.
| 层数 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|---|
| 1 | \(a_0\) | \(a_1\) | \(a_2\) | \(a_3\) | \(a_4\) | \(a_5\) | \(a_6\) | \(a_7\) |
| 2 | \(a_0\) | \(a_2\) | \(a_4\) | \(a_6\) | \(a_1\) | \(a_3\) | \(a_5\) | \(a_7\) |
| 3 | \(a_0\) | \(a_4\) | \(a_2\) | \(a_6\) | \(a_1\) | \(a_5\) | \(a_3\) | \(a_7\) |
| 4 | \(a_0\) | \(a_4\) | \(a_2\) | \(a_6\) | \(a_1\) | \(a_5\) | \(a_3\) | \(a_7\) |
| 发现位置跟这个数的二进制有关,于是将它们最终位置固定住即可(本质上是一个二进制翻转) |
代码
#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 double
inline 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)原根
定义:\(a^r\equiv{1}\pmod{P}\)若最小的\(r=\phi(m)\)则称\(a\)为模\(P\)的原根
令\(P=qn+1\)(n为2的幂次),定义g为原根,则\(\omega_{n}=g^q\),集合为\((1,g^q,g^{2q},...,g^{(n-1)q})\)
既然被用来替代单位根,那么其就具备单位根应该有的性质(性质对应上面的FFT)
性质一:根据原根的定义可知,\(g^{\phi(P)}\equiv{1}\pmod{P}\), \(\phi(P)\)为满足条件中最小的数,因此这些数两两不同
性质二:\(\omega_{2n}=g^{q/2}(P=q/2*2n+1)\)因此得证
性质三:\(\frac{\omega_{n}^{k+n/2}}{\omega_{n}^{k}}=\omega_{n}^{n/2}\),\((\omega_{n}^{n/2})^2\equiv{1} \pmod{P}\)因为\(\omega_{n}^{n/2}!=1\)所以\(\omega_{n}^{n/2}=-1\)得证
性质四:易证
性质五:同理用等比数列求和即可
(2)找原根
那么我们该怎么寻找原根,是不是所有数都有原根?
首先寻找一个质数的原根,我们可以用枚举法来做
有个重要的小结论:若\(a^r\equiv{1}\pmod{P}\)最小的r一定是\(P-1\)的约数,因为由费马小定理可得\(a^{P-1}\equiv{1}\pmod{P}\)恒成立(也可利用反证法)
\(1000000007\)的原根是\(5\)
\(998244353\)的原根是\(3\)
\(7340033\)的原根是\(3\)
\(5767169\)的原根是\(3\)
其它需要的话用此程序(跑得比较慢,随便扒下来的)
#include<bits/stdc++.h>
using namespace std;
#define LL long long
using 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;}
}
}
update
不得已自己写了一个.(因为每次都要翻题解看原根是多少)
LL fpow(LL a,int cnt) {
LL res=1;
for(int i=cnt; i; i>>=1,a=a*a%P)if(i&1)res=res*a%P;
return res;
}
namespace Get_G {
int stk[M],top;
bool Check(int now) {
for(int i=1; i<=top; ++i) {
int y=(P-1)/stk[i];
if(fpow(now,y)==1)return false;
}
return true;
}
void Solve() {
int p=sqrt(P-1),x=P-1;
for(int i=2; i<=p&&i<=x; ++i) {
if(x%i)continue;
stk[++top]=i;
while(x%i==0)x/=i;
}if(x>1)stk[++top]=x;
for(int i=2; i<P; ++i) {
if(Check(i)) {
printf("%d",i);
break;
}
}
}
}
代码
#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
inline 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;
}
浙公网安备 33010602011771号