【演算法講15:多項式卷積】
- 模板題
- 多項式的表示法
- 多項式的乘法與卷積
- 多項式乘法
- ? \lceil ?快速傅里葉變換 ? \rfloor ? F F T FFT FFT
- 單位根
- 離散傅里葉變換 D F T DFT DFT
- D F T DFT DFT 的詳細討論
- 離散傅里葉逆變換 I D F T IDFT IDFT
- 減小常數 (蝴蝶變換)
- ? \lceil ?快速數論變換 ? \rfloor ? (NTT)
- 原根
- 同樣的性質
- 參考與查閱
模板題
- 【模板】多項式乘法(FFT) | 洛谷 P3803
求兩個多項式的乘積
多項式的表示法
- 系數表示:
n
n
n 個系數確定一個
n
?
1
n-1
n?1 次多項式
( a 0 , a 1 , a 2 , ? ? , a n ? 1 ) (a_0,a_1,a_2,\cdots,a_{n-1}) (a0?,a1?,a2?,?,an?1?) - 點值表示:
n
n
n 個點確定一個
n
?
1
n-1
n?1 刺多項式
( x 0 , x 1 , x 1 , ? ? , x n ? 1 ) ( y 0 , y 1 , y 1 , ? ? , y n ? 1 ) (x_0,x_1,x_1,\cdots,x_{n-1})\\ (y_0,y_1,y_1,\cdots,y_{n-1})\\ (x0?,x1?,x1?,?,xn?1?)(y0?,y1?,y1?,?,yn?1?)
其中滿足
y i = ∑ j = 0 n ? 1 a j x i j y_i=\sum_{j=0}^{n-1}a_jx_i^j yi?=j=0∑n?1?aj?xij?
多項式的乘法與卷積
多項式乘法
- 有兩個多項式 A ( x ) , B ( x ) A(x),B(x) A(x),B(x),兩個多項式相乘結果為 C ( x ) C(x) C(x),
- 新的多項式系數可以表示為:
c k = ∑ i + j = k a i b j = ∑ i = 0 k a i b k ? i c_k=\sum_{i+j=k}a_ib_j=\sum_{i=0}^ka_ib_{k-i} ck?=i+j=k∑?ai?bj?=i=0∑k?ai?bk?i? - 多項式乘法得到的結果又可以叫做這兩個多項式的卷積,可以用 F F T FFT FFT 進行快速計算,
? \lceil ?快速傅里葉變換 ? \rfloor ? F F T FFT FFT
- 如果暴力求解兩個多項式的卷積,時間復雜度為
O
(
n
m
)
O(nm)
O(nm)
利用 F F T FFT FFT 可以把時間復雜度降到 O ( n log ? n ) O(n\log n) O(nlogn)
具體分為兩個程序: - 把 A ( x ) 、 B ( x ) A(x)、B(x) A(x)、B(x) 的系數表示轉化成點值表示,然后用 O ( n ) O(n) O(n) 的方法求得 C ( x ) C(x) C(x) 的點值表示,然后還原成 C ( x ) C(x) C(x) 的系數表示,
- 那么第一步,我們想要找一個快速的方法把 A ( x ) 、 B ( x ) A(x)、B(x) A(x)、B(x) 的點值表示求出來,
單位根
- F F T FFT FFT 的方法是找到了一組特殊的點代到多項式當中去,可以快速找到一組點的數值,那組特殊的點就是復數中的單位根,
- 一個復數可以在復平面上用一個向量來表示,在單位元上的向量都可以成為單位根,
我們把一個單位元分成 n n n 份,我們就得到了 n n n 個單位根,其中的第 k k k 個可以即為 ω n k \omega_n^k ωnk?,并且 ω n 0 = 1 \omega_n^0=1 ωn0?=1 - 歐拉公式:
e
i
θ
=
c
o
s
θ
+
i
?
sin
?
θ
e^{i\theta}=cos\theta+i\cdot\sin\theta
eiθ=cosθ+i?sinθ
于是,
ω n k = e 2 π k n i \omega_n^k=e^{2\pi\frac{k}{n}i} ωnk?=e2πnk?i - 性質:
ω n k = ω 2 n 2 k \omega_n^k=\omega_{2n}^{2k} ωnk?=ω2n2k?
ω n k + 2 n = ? ω n k \omega_n^{k+\frac{2}{n}}=-\omega_n^k ωnk+n2??=?ωnk?
ω n k ? ω n l = ω n k + l \omega_n^k\cdot\omega_n^l=\omega_n^{k+l} ωnk??ωnl?=ωnk+l?
離散傅里葉變換 D F T DFT DFT
- 考慮 n = 2 b n=2^b n=2b 項的多項式 A ( x ) A(x) A(x) 系數表示為 ( a 0 , a 1 , ? ? , a n ? 1 ) (a_0,a_1,\cdots,a_{n-1}) (a0?,a1?,?,an?1?)
- 我們現在想匯入一組單位根 ( ω n 0 , ω n 1 , ? ? , ω n n ? 1 ) (\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}) (ωn0?,ωn1?,?,ωnn?1?),快速算出 ( A ( ω n 0 ) , A ( ω n 1 ) , ? ? , A ( ω n n ? 1 ) ) (A(\omega_n^0),A(\omega_n^1),\cdots,A(\omega_n^{n-1})) (A(ωn0?),A(ωn1?),?,A(ωnn?1?)),這個程序成為 D F T DFT DFT
- 首先,我們把多項式按奇偶性分類:
A ( x ) = ( a 0 + a 2 x 2 + ? + a n ? 2 x n ? 2 ) + x ( a 1 + a 3 x 2 + ? + a n ? 1 x n ? 2 ) A(x)=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2}) A(x)=(a0?+a2?x2+?+an?2?xn?2)+x(a1?+a3?x2+?+an?1?xn?2)
現在,我們設兩個多項式:
A 0 ( x ) = a 0 + a 2 x 1 + ? + a n ? 2 x n ? 2 2 A 1 ( x ) = a 1 + a 3 x 1 + ? + a n ? 1 x n ? 2 2 A0(x)=a_0+a_2x^1+\cdots+a_{n-2}x^{\frac{n-2}{2}}\\ A1(x)=a_1+a_3x^1+\cdots+a_{n-1}x^{\frac{n-2}{2}} A0(x)=a0?+a2?x1+?+an?2?x2n?2?A1(x)=a1?+a3?x1+?+an?1?x2n?2?
于是,原來的多項式就變為:
A ( x ) = A 0 ( x 2 ) + x A 1 ( x 2 ) A(x)=A0(x^2)+xA1(x^2) A(x)=A0(x2)+xA1(x2)
D F T DFT DFT 的詳細討論
- 當
0
≤
k
≤
n
2
?
1
0\le k\le \frac{n}{2}-1
0≤k≤2n??1 時,
A ( ω n k ) = A 0 ( ( ω n k ) 2 ) + ω n k ? A 1 ( ( ω n k ) 2 ) A(\omega_n^k)=A0((\omega_n^k)^2)+\omega_n^k\cdot A1((\omega_n^k)^2) A(ωnk?)=A0((ωnk?)2)+ωnk??A1((ωnk?)2)
化簡得到:
A ( ω n k ) = A 0 ( ω n 2 k ) + ω n k ? A 1 ( ω n 2 k ) A(\omega_n^k)=A0(\omega_{\frac{n}{2}}^k)+\omega_n^k\cdot A1(\omega_{\frac{n}{2}}^k) A(ωnk?)=A0(ω2n?k?)+ωnk??A1(ω2n?k?)
可以看到,我們分為了兩個子問題,可以遞回求解, - 那么
n
2
≤
k
+
n
2
≤
n
?
1
\frac{n}{2}\le k+\frac{n}{2}\le n-1
2n?≤k+2n?≤n?1,我們帶入得到:
A ( ω n k + n 2 ) = A 0 ( ( ω n k + n 2 ) 2 ) + ω n k + n 2 ? A 1 ( ( ω n k + n 2 ) 2 ) A(\omega_n^{k+\frac{n}{2}})=A0((\omega_n^{k+\frac{n}{2}})^2)+\omega_n^{k+\frac{n}{2}}\cdot A1((\omega_n^{k+\frac{n}{2}})^2) A(ωnk+2n??)=A0((ωnk+2n??)2)+ωnk+2n???A1((ωnk+2n??)2)
化簡得到:
A ( ω n k + n 2 ) = A 0 ( ω n 2 k ) ? ω n k ? A 1 ( ω n 2 k ) A(\omega_n^{k+\frac{n}{2}})=A0(\omega_{\frac{n}{2}}^k)-\omega_n^k\cdot A1(\omega_{\frac{n}{2}}^k) A(ωnk+2n??)=A0(ω2n?k?)?ωnk??A1(ω2n?k?) - 時間復雜度: T ( n ) = 2 T ( n 2 ) + O ( n ) = O ( n log ? n ) T(n)=2T(\frac{n}{2})+O(n)=O(n\log n) T(n)=2T(2n?)+O(n)=O(nlogn)
離散傅里葉逆變換 I D F T IDFT IDFT
- 將 D F T DFT DFT 得到的點值轉化為系數,這一程序就叫做 I D F T IDFT IDFT
- 假設我們得到了
A
(
x
)
A(x)
A(x) 的點值,設為
(
d
0
,
d
1
,
?
?
,
d
n
?
1
)
(d_0,d_1,\cdots,d_{n-1})
(d0?,d1?,?,dn?1?),我們要獲得其多項式
F
(
x
)
F(x)
F(x)
設 c ( k ) c(k) c(k) 為 x = ω n ? k x=\omega_n^{-k} x=ωn?k? 時的點值,
c k = ∑ i = 0 n ? 1 d i ? ( ω n ? k ) i = ∑ i = 0 n ? 1 ∑ j = 0 n ? 1 a j ? ( ω n i ) j ? ( ω n ? k ) i = ∑ j = 0 n ? 1 a j ∑ i = 0 n ? 1 ( ω n j ? k ) i = ∑ j = 0 n ? 1 a j ? n ? [ j = k ] = n ? a k c_k=\sum_{i=0}^{n-1}d_i\cdot (\omega_n^{-k})^i=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\cdot (\omega_n^i)^j\cdot (\omega_n^{-k})^i=\sum_{j=0}^{n-1}a_j \sum_{i=0}^{n-1}(\omega_n^{j-k})^i=\sum_{j=0}^{n-1}a_j\cdot n\cdot [j=k]=n\cdot a_k ck?=i=0∑n?1?di??(ωn?k?)i=i=0∑n?1?j=0∑n?1?aj??(ωni?)j?(ωn?k?)i=j=0∑n?1?aj?i=0∑n?1?(ωnj?k?)i=j=0∑n?1?aj??n?[j=k]=n?ak?
減小常數 (蝴蝶變換)
- 我們一遍遍遞回下去,在最后一層發現一個規律:原序列和置換后的序列的對應二進制位剛好是翻轉的,比如 1 → 8 1\rightarrow 8 1→8, 14 → 7 14\rightarrow 7 14→7
- 代碼:
2.77 s 2.77s 2.77s
const int MAX = 1e7+50;
const double PI = acos(-1);
struct comp{
double x,y;
comp(double xx = 0,double yy = 0){x = xx,y = yy;}
};
comp operator +(comp a,comp b){return comp(a.x+b.x , a.y+b.y);}
comp operator -(comp a,comp b){return comp(a.x-b.x , a.y-b.y);}
comp operator *(comp a,comp b){return comp(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);}
int lim = 1;
int l,r[MAX];
void FFT(comp *A,int type){
for(int i = 0;i < lim;++i)
if(i < r[i])swap(A[i],A[r[i]]);
for(int mid = 1;mid < lim;mid<<=1){
comp omg(cos(PI/mid),type * sin(PI/mid));
for(int R = mid<<1,j = 0;j < lim;j += R){
comp w(1,0);
for(int k = 0;k < mid;++k,w = w * omg){
comp x = A[j+k],y = w * A[j + mid + k];
A[j + k] = x + y;
A[j + mid + k] = x - y;
}
}
}
}
comp aa[MAX],bb[MAX];
int main()
{
int N,M;N = read();M = read();
for(int i = 0;i <= N;++i)aa[i].x = read();
for(int i = 0;i <= M;++i)bb[i].x = read();
while(lim <= N + M)lim<<=1,l++;
for(int i = 0;i < lim;++i)
r[i] = (r[i>>1]>>1) | ((i&1)<<(l-1));
FFT(aa,1);
FFT(bb,1);
for(int i = 0;i <= lim;++i)aa[i] = aa[i] * bb[i];
FFT(aa,-1);
for(int i = 0;i <= N + M;++i)
printf("%d ",(int)(aa[i].x/lim+0.5));
return 0;
}
? \lceil ?快速數論變換 ? \rfloor ? (NTT)
- 求的是多項式乘法的系數取模之后的結果,這里模數應該要是一個質數,
- 我們希望帶入一組冪次以便于更好取模運算,
原根
- 如果
g
i
?
m
o
d
?
p
(
1
≤
i
≤
p
?
1
)
g^i\bmod p (1\le i\le p-1)
gimodp(1≤i≤p?1) 的值互不相同,那么
g
g
g 稱為
p
p
p 的原根,
并非所有的數都有原根,且如果一個數存在原根,其個數也不唯一, - 事實上,原根存在的重要條件是:
m
=
2
,
4
,
p
a
,
2
p
a
(
a
>
0
)
m=2,4,p^a,2p^a(a>0)
m=2,4,pa,2pa(a>0) 且
a
a
a 為奇素數,
如果這個數有原根,那么一定有 φ ( φ ( m ) ) \varphi(\varphi(m)) φ(φ(m)) 個原根, - 原根沒有快速求法,只能通過暴力列舉判斷,或者查表,
常見的 N T T NTT NTT 模數是 998244353 998244353 998244353,其原根為 3 3 3
同樣的性質
- 加入模數為
p
p
p,我們發現構造一下原根,得到與單位根有同樣的性質的數,即:
ω n i = g i ( p ? 1 ) n ? m o d ? p \omega_n^i=g^{\frac{i(p-1)}{n}}\bmod p ωni?=gni(p?1)?modp
直接將 F F T FFT FFT 中的原根替換為 N T T NTT NTT 中的原根, - 代碼
1.64 s 1.64s 1.64s
const int MAX = 1e7+50;
ll qpow(ll a,ll n){/* */ll res = 1LL;while(n){if(n&1)res=res*a%MOD;a=a*a%MOD;n>>=1;}return res;}
const int P = 998244353,G = 3,Gi = 332748118;
int lim = 1;
int l,r[MAX];
void NTT(ll *A,int type){
for(int i = 0;i < lim;++i)
if(i < r[i])swap(A[i],A[r[i]]);
for(int mid = 1;mid < lim;mid<<=1){
ll omg = qpow(type == 1 ? G : Gi , (P - 1) / (mid << 1));
for(int j = 0;j < lim;j += (mid << 1)){
ll w = 1;
for(int k = 0;k < mid;++k,w = w * omg % P){
int x = A[j+k],y = w * A[j + mid + k] % P;
A[j + k] = (x + y) % P;
A[j + mid + k] = (x - y + P) % P;
}
}
}
}
ll aa[MAX],bb[MAX];
int main()
{
int N,M;N = read();M = read();
for(int i = 0;i <= N;++i)aa[i] = (read() + P) % P;
for(int i = 0;i <= M;++i)bb[i] = (read() + P) % P;
while(lim <= N + M)lim<<=1,l++;
for(int i = 0;i < lim;++i)
r[i] = (r[i>>1]>>1) | ((i&1)<<(l-1));
NTT(aa,1);
NTT(bb,1);
for(int i = 0;i < lim;++i)aa[i] = (aa[i] * bb[i]) % P;
NTT(aa,-1);
ll iv = qpow(lim,P - 2);
for(int i = 0;i <= N + M;++i)
printf("%d ",(aa[i] * iv) % P);
return 0;
}
參考與查閱
- 【SDUACM-暑期專題div1】FFT、NTT
- 洛谷改題的 d a l a o dalao dalao 們的題解
轉載請註明出處,本文鏈接:https://www.uj5u.com/qita/261091.html
標籤:其他
上一篇:用python玩微信跳一跳
下一篇:關于聯邦濾波器
