模板:多项式乘法(FFTNTT)
文章目錄
- 前言
- 系數表示法和點值表示法
- 單位根
- 離散傅立葉變換(DFT)
- 位逆序置換(蝴蝶變換)
- 離散傅立葉逆變換(IDFT)
- 代碼
- NTT
- 代碼
所謂快速傅立葉變換,就是傅立葉發明的一種快速的變換
(逃)
前言
FFT,用于在nlognnlognnlogn的復雜度內求兩個多項式相乘。
算是多項式的入門吧。
十一被rabbit的數學打蒙了所以開始碰這些東西。
不知道以后會不會直接鴿掉。
(update on 2022.1.4:三個月后終于回來填坑了)。
很妙的算法。
代碼寫成迭代,精煉之后很優美。
系數表示法和點值表示法
系數表示法就是我們常用的函數的表示形式。
對于一個n次的函數,點值表示法則是給出了n+1個互異的點,通過這n+1個互異的點就可以唯一確定的描述出這個函數(可以想想高斯消元)。
這東西有啥用?
考慮如果兩個多項式都表示成了點值形式,且選取的點的橫坐標相同,我把它們的縱坐標對應乘起來,就能得到它們乘積的多項式的點值表示
乘起來的復雜度降到了O(n)O(n)O(n)
但是這個點值表示法求起來暴力似乎還是n方啊…
所以快速傅立葉變換其實解決的就是系數表示法和點值表示法快速互相轉換的問題
單位根
規定:
若復數ω\omegaω滿足ωn=1\omega^n=1ωn=1,就稱ω\omegaω是n次單位根
對于一個n項的多項式,我們讓它的點值表示取的橫坐標為ωnk(0≤k<n)\omega_n^k(0\leq k <n)ωnk?(0≤k<n)
然后就會有很多很妙的性質:
都可以結合幾何意義較為顯然的證明。
離散傅立葉變換(DFT)
那么回到問題。
首先,我們要把多項式填零直到 n=2kn=2^kn=2k 的長度。
接下來,我們要求:
f(x)=a0+a1x+a2x2+a3x3+...+an?1xn?1f(x)=a_0+a_1x+a_2x^2+a_3x^3+...+a_{n-1}x^{n-1}f(x)=a0?+a1?x+a2?x2+a3?x3+...+an?1?xn?1
讓我們把這個函數分奇偶拆成兩個函數:
A(x)=a0+a2x+a4x2+a6x3+...+an?1xn2A(x)=a_0+a_2x+a_4x^2+a_6x^3+...+a_{n-1}x^{\frac{n}{2}}A(x)=a0?+a2?x+a4?x2+a6?x3+...+an?1?x2n?
B(x)=a1+a3x+a5x2+a7x3+...+an?2xn2B(x)=a_1+a_3x+a_5x^2+a_7x^3+...+a_{n-2}x^{\frac{n}{2}}B(x)=a1?+a3?x+a5?x2+a7?x3+...+an?2?x2n?
那么就有:
f(x)=A(x2)+xB(x2)f(x)=A(x^2)+xB(x^2)f(x)=A(x2)+xB(x2)
我們把 x=ωnkx=\omega_n^kx=ωnk? 帶入。
分情況討論:
(設:0≤k<n20\le k<\dfrac{n}{2}0≤k<2n?)
當指數 <n2< \dfrac{n}{2}<2n? 時:
f(ωnk)=A((ωnk)2)+ωnkB((ωnk)2)f(\omega_n^k)=A((\omega_n^k)^2)+\omega_n^kB((\omega_n^k)^2)f(ωnk?)=A((ωnk?)2)+ωnk?B((ωnk?)2)
=A(ωn2k)+ωnkB(ωn2k)=A(\omega_n^{2k})+\omega_n^kB(\omega_n^{2k})=A(ωn2k?)+ωnk?B(ωn2k?)
=A(ωn2k)+ωnkB(ωn2k)=A(\omega_{\frac{n}{2}}^{k})+\omega_n^kB(\omega_{\frac{n}{2}}^{k})=A(ω2n?k?)+ωnk?B(ω2n?k?)
類似的,當指數 ≥n2\ge\dfrac{n}{2}≥2n? 時:
f(ωnk+n2)=A(ωn2k+n)+ωnk+n2B(ωn2k+n)f(\omega_n^{k+\frac{n}{2}})=A(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}B(\omega_n^{2k+n})f(ωnk+2n??)=A(ωn2k+n?)+ωnk+2n??B(ωn2k+n?)
=A(ωn2k)?ωnkB(ωn2k)=A(\omega_n^{2k})-\omega_n^{k}B(\omega_n^{2k})=A(ωn2k?)?ωnk?B(ωn2k?)
=A(ωn2k)?ωnkB(ωn2k)=A(\omega_{\frac{n}{2}}^{k})-\omega_n^{k}B(\omega_{\frac{n}{2}}^{k})=A(ω2n?k?)?ωnk?B(ω2n?k?)
這樣,我們就可以在 O(len)O(len)O(len) 的復雜度內合并兩個長度為 lenlenlen 的多項式,就可以利用分治把復雜度降到 O(nlog?n)O(n\log n)O(nlogn)。
位逆序置換(蝴蝶變換)
這個也很妙
有點線性求逆元那味了
通過觀察發現,系數 i 調換位置之后會去往它的二進制表示反過來的位置(記為rir_iri?)
而上面那個東西可以通過:
r[i]=r[i>>1]>>1+[i&1]?(n>>1)r[i]=r[i>>1]>>1+[i\&1]*(n>>1) r[i]=r[i>>1]>>1+[i&1]?(n>>1)
來線性的求
這個還是挺好理解的畫一畫就行了
離散傅立葉逆變換(IDFT)
我們變完之后當然還要變回來。
我會 n3n^3n3 高斯消元!
設之前DFT得到的點值表示 bk=f(ωnk)=∑i=0n?1ai(ωnk)ib_k=f(\omega_n^k)=\sum_{i=0}^{n-1}a_i(\omega_{n}^k)^ibk?=f(ωnk?)=∑i=0n?1?ai?(ωnk?)i
構造一個函數:
g(x)=b0+b1x+b2x2+b3x3+...+bn?1xn?1g(x)=b_0+b_1x+b_2x^2+b_3x^3+...+b_{n-1}x_{n-1}g(x)=b0?+b1?x+b2?x2+b3?x3+...+bn?1?xn?1?
然后我們用 ωn0,ωn?1,ωn?2,...,ωn?n\omega_n^{0},\omega_n^{-1},\omega_n^{-2},...,\omega_n^{-n}ωn0?,ωn?1?,ωn?2?,...,ωn?n? 作為橫坐標帶入。(不難發現剛才 DFT 需要的性質依然成立,所以可以用同樣的方法求解)
那么我們在第 kkk 位得到的新的點值就是:
g(ωn?k)=∑i=0n?1bi(ωn?k)ig(\omega_n^{-k})=\sum_{i=0}^{n-1}b_i(\omega_n^{-k})^ig(ωn?k?)=i=0∑n?1?bi?(ωn?k?)i
=∑i=0n?1∑j=0n?1aj(ωni)j(ωn?k)i=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_{n}^i)^j(\omega_n^{-k})^i=i=0∑n?1?j=0∑n?1?aj?(ωni?)j(ωn?k?)i
=∑j=0j?1aj∑i=0n?1(ωnj?k)i=\sum_{j=0}^{j-1}a_j\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i=j=0∑j?1?aj?i=0∑n?1?(ωnj?k?)i
考慮 ∑i=0n?1(ωnj?k)i\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i∑i=0n?1?(ωnj?k?)i,這一項,當 j=kj=kj=k 時,顯然等于0;當 j≠kj\ne kj?=k 時,根據等比公式,有:
∑i=0n?1(ωnj?k)i=(wnj?k)n?1wnj?k?1\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i=\frac{(w_n^{j-k})^n-1}{w_n^{j-k}-1}i=0∑n?1?(ωnj?k?)i=wnj?k??1(wnj?k?)n?1?
=(wnn)j?k?1wnj?k?1=1j?k?1wnj?k?1=0=\frac{(w_n^{n})^{j-k}-1}{w_n^{j-k}-1}=\frac{1^{j-k}-1}{w_n^{j-k}-1}=0=wnj?k??1(wnn?)j?k?1?=wnj?k??11j?k?1?=0
所以原來的式子只有當 j=kj=kj=k 的時候有值,等于 nakna_knak?。
所以我們直接 NTT 完除 nnn 即可。
代碼
#include<bits/stdc++.h> using namespace std; #define ll long long #define il inline const int N=5e6+100; const int M=150; const int mod=998244353; const double pi=acos(-1.0); inline ll read(){ll x=0,f=1;char c=getchar();while(!isdigit(c)){if(c=='-') f=-1;c=getchar();}while(isdigit(c)){x=x*10+c-'0';c=getchar();}return x*f; } int n,m,len,k; struct node{double x,y;node(double a=0,double b=0){x=a;y=b;} }A[N],B[N]; il node operator * (node a,node b){return (node){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x}; } il node operator + (node a,node b){return (node){a.x+b.x,a.y+b.y}; } il node operator - (node a,node b){return (node){a.x-b.x,a.y-b.y}; } int r[N]; il void fft(node *x,int flag){for(int i=0;i<=len;i++){if(i<r[i]) swap(x[i],x[r[i]]);}for(int l=1;l<len;l<<=1){node o(cos(pi/l),flag*sin(pi/l));for(int st=0;st<len;st+=l<<1){node t(1,0);for(int j=0;j<l;j++,t=t*o){node u=x[st+j],v=t*x[st+j+l];x[st+j]=u+v;x[st+j+l]=u-v;}}}return; } int main(){n=read();m=read();for(int i=0;i<=n;i++) A[i].x=read();for(int i=0;i<=m;i++) B[i].x=read();len=n+m;while((1<<k)<=len) k++;len=1<<k;for(int i=1;i<=len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));fft(A,1);fft(B,1);for(int i=0;i<=len;i++) A[i]=A[i]*B[i];fft(A,-1);for(int i=0;i<=n+m;i++) printf("%d ",(int)(A[i].x/(len)+0.5));return 0; } /**/NTT
NTT說簡單也簡單說難也難。
簡單是因為板子幾乎和FFT一模一樣。
難是因為那個群的概念完全看不懂qwq。
然后我就選擇調過了證明部分。
背下來背下來!——宋隊
總的來說,寫一寫對背板子有幫助的我可以理解的。
對于一個質數P,若其存在原根(基本全是3)且可以表示為:
P=k?2n+1P=k*2^n+1P=k?2n+1
它就是一個ntt模數。
ggg為什么叫原根?
因為它有一些關鍵性質:
1.gk?n=1(modP)1. g^{k*n}=1(mod P)1.gk?n=1(modP)
2.gi(modP)g^i(mod P)gi(modP) 在 0≤i≤p0\leq i\leq p0≤i≤p 的范圍內取值兩兩不同。
設:gkg^kgk 為單位根,記為gng_ngn?。
那么它就有很多熟悉的性質:
沒錯,推FFT時 ω\omegaω 需要有的性質它全有!
所以就可以直接上FFT的板子啦!
然后煩人的復數、浮點運算和精度問題全都拜拜了。
針不戳。
宋隊誠不欺我。
逆變換還是倒數,從FFT的求共軛改成求逆元。
就像喝水一樣。
代碼
il void ntt(ll *x,int flag){for(int i=0;i<lim;i++){if(i<r[i]) swap(x[i],x[r[i]]);}for(register int l=1;l<lim;l<<=1){register ll g=ksm(3,(mod-1)/(l<<1));if(flag==-1) g=ksm(g,mod-2);for(register int st=0;st<lim;st+=l<<1){register ll now=1;for(register int i=0;i<l;i++,now=now*g%mod){ll u=x[st+i],v=x[st+i+l]*now%mod;x[st+i]=(u+v)%mod;x[st+i+l]=(u-v+mod)%mod;}}}if(flag==-1){register ll ni=ksm(lim,mod-2);for(register int i=0;i<lim;i++) x[i]=x[i]*ni%mod;} }總結
以上是生活随笔為你收集整理的模板:多项式乘法(FFTNTT)的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: LOL日女技能介绍
- 下一篇: 逆水寒岁月神偷最后奖励选什么好?奖励一览