捲積方法 OAO 本篇的重點應該是放在理解與使用與競程相關的捲積，其中會利用到FFT或NTT加速∗ * ∗
( a ∗ b ) x = ∑ i + j = x a i b j
(a * b) _ x = \sum _ {i+j=x} a_i b_j
( a ∗ b ) x = i + j = x ∑ a i b j 實際上就等價於我們常見的多項式乘法O ( n 2 ) \mathcal{O}(n^2) O ( n 2 )
1
2
3
4
5
for ( int i = 0 ; i < A . size (); i ++ ) {
for ( int j = 0 ; j < B . size (); j ++ ) {
res [ i + j ] += A [ i ] * B [ j ];
}
}
DFT 先假設我們有兩個多項式
A ( x ) = ∑ a i x i , B ( x ) = ∑ b i x i C ( x ) = A ( x ) B ( x )
\begin{gathered}
A(x) = \sum a_i x^i, B(x) = \sum b_i x^i\\
C(x) = A(x)B(x)
\end{gathered}
A ( x ) = ∑ a i x i , B ( x ) = ∑ b i x i C ( x ) = A ( x ) B ( x ) C C C A ( x ) A(x) A ( x ) B ( x ) B(x) B ( x ) n n n C C C n = deg C + 1 n = \deg C + 1 n = deg C + 1
[ 1 x 0 x 0 2 ⋯ x 0 n − 1 1 x 1 x 1 2 ⋯ x 1 n − 1 ⋮ ⋮ ⋮ ⋱ ⋮ 1 x n − 1 x n − 1 2 ⋯ x n − 1 n − 1 ] [ a 0 a 1 ⋮ a n − 1 ] = [ A ( x 0 ) A ( x 1 ) ⋮ A ( x n − 1 ) ]
\begin{bmatrix}
1 & x_0 & x_0^2 & \cdots & x_0^{n-1}\newline
1 & x_1 & x_1^2 & \cdots & x_1^{n-1}\newline
\vdots & \vdots & \vdots & \ddots & \vdots\newline
1 & x _ {n-1} & x _ {n-1}^2 & \cdots & x _ {n-1}^{n-1}
\end{bmatrix}
\begin{bmatrix}
a_0\newline a_1\newline \vdots\newline a _ {n-1}
\end{bmatrix} =
\begin{bmatrix}
A(x_0)\newline A(x_1)\newline \vdots\newline A(x _ {n-1})
\end{bmatrix}
1 1 ⋮ 1 x 0 x 1 ⋮ x n − 1 x 0 2 x 1 2 ⋮ x n − 1 2 ⋯ ⋯ ⋱ ⋯ x 0 n − 1 x 1 n − 1 ⋮ x n − 1 n − 1 a 0 a 1 ⋮ a n − 1 = A ( x 0 ) A ( x 1 ) ⋮ A ( x n − 1 ) [ 1 x 0 x 0 2 ⋯ x 0 n − 1 1 x 1 x 1 2 ⋯ x 1 n − 1 ⋮ ⋮ ⋮ ⋱ ⋮ 1 x n − 1 x n − 1 2 ⋯ x n − 1 n − 1 ] [ b 0 b 1 ⋮ b n − 1 ] = [ B ( x 0 ) B ( x 1 ) ⋮ B ( x n − 1 ) ]
\begin{bmatrix}
1 & x_0 & x_0^2 & \cdots & x_0^{n-1}\newline
1 & x_1 & x_1^2 & \cdots & x_1^{n-1}\newline
\vdots & \vdots & \vdots & \ddots & \vdots\newline
1 & x _ {n-1} & x _ {n-1}^2 & \cdots & x _ {n-1}^{n-1}
\end{bmatrix}
\begin{bmatrix}
b_0\newline b_1\newline \vdots\newline b _ {n-1}
\end{bmatrix} =
\begin{bmatrix}
B(x_0)\newline B(x_1)\newline \vdots\newline B(x _ {n-1})
\end{bmatrix}
1 1 ⋮ 1 x 0 x 1 ⋮ x n − 1 x 0 2 x 1 2 ⋮ x n − 1 2 ⋯ ⋯ ⋱ ⋯ x 0 n − 1 x 1 n − 1 ⋮ x n − 1 n − 1 b 0 b 1 ⋮ b n − 1 = B ( x 0 ) B ( x 1 ) ⋮ B ( x n − 1 ) [ 1 x 0 x 0 2 ⋯ x 0 n − 1 1 x 1 x 1 2 ⋯ x 1 n − 1 ⋮ ⋮ ⋮ ⋱ ⋮ 1 x n − 1 x n − 1 2 ⋯ x n − 1 n − 1 ] − 1 [ A ( x 0 ) B ( x 0 ) A ( x 1 ) B ( x 1 ) ⋮ A ( x n − 1 ) B ( x n − 1 ) ] = [ c 0 c 1 ⋮ c n − 1 ]
\begin{bmatrix}
1 & x_0 & x_0^2 & \cdots & x_0^{n-1}\newline
1 & x_1 & x_1^2 & \cdots & x_1^{n-1}\newline
\vdots & \vdots & \vdots & \ddots & \vdots\newline
1 & x _ {n-1} & x _ {n-1}^2 & \cdots & x _ {n-1}^{n-1}
\end{bmatrix}
^{-1}
\begin{bmatrix}
A(x_0)B(x_0)\newline A(x_1)B(x_1)\newline \vdots\newline A(x _ {n-1})B(x _ {n-1})
\end{bmatrix} =
\begin{bmatrix}
c_0\newline c_1\newline \vdots\newline c _ {n-1}
\end{bmatrix}
1 1 ⋮ 1 x 0 x 1 ⋮ x n − 1 x 0 2 x 1 2 ⋮ x n − 1 2 ⋯ ⋯ ⋱ ⋯ x 0 n − 1 x 1 n − 1 ⋮ x n − 1 n − 1 − 1 A ( x 0 ) B ( x 0 ) A ( x 1 ) B ( x 1 ) ⋮ A ( x n − 1 ) B ( x n − 1 ) = c 0 c 1 ⋮ c n − 1 最後一步直接高斯消去或用拉格朗日/牛頓插值法可以做到O ( n 2 ) \mathcal{O}(n^2) O ( n 2 )
上述步驟稱為DFT(對序列的版本叫離散傅立葉變換，與使用積分的連續傅立葉變換相對)和IDFTx x x
Root of Unity 首先先來介紹單位根ω \omega ω
ω n = 1 ∀ 0 ≤ i < j < n , ω i ≠ ω j
\begin{gathered}
\omega ^ n = 1\\
\forall 0 \leq i < j < n, \omega^i \neq \omega^j
\end{gathered}
ω n = 1 ∀0 ≤ i < j < n , ω i = ω j 的數
複習一下歐拉公式
e i x = cos x + i sin x e^{ix} = \cos x + i\sin x e i x = cos x + i sin x ω n = e − 2 π i n \omega_n = e^{-\frac{2\pi i}{n}} ω n = e − n 2 π i n n n i i i ω i = 1 \omega^i = 1 ω i = 1 ord ( ω n ) = n \operatorname{ord}(\omega_n) = n ord ( ω n ) = n
推薦觀賞3B1B系列https://youtu.be/v0YEaeIClKY https://youtu.be/mvmuCPvRoWQ
引理們 Lemma a. ω d n d k = ( e 2 π i d n ) d k = ( e 2 π i n ) k = ω n k
\omega _ {dn}^{dk} = (e^{\frac{2\pi i}{dn}})^{dk} = (e^{\frac{2\pi i}{n}})^k = \omega_n^k
ω d n d k = ( e d n 2 π i ) d k = ( e n 2 π i ) k = ω n k Lemma b. ω n n 2 = ω 2 = e i π = − 1
\omega_n^{\frac{n}{2}} = \omega_2 = e^{i\pi} = -1
ω n 2 n = ω 2 = e iπ = − 1 Cooley-Tukey FFT algorithm 先假設n n n i i i x i x_i x i ω n i \omega_n^i ω n i i ∈ [ 0 , n − 1 ] i \in [0, n-1] i ∈ [ 0 , n − 1 ]
y i = ∑ j = 0 n − 1 a j ( ω n i ) j
y_i = \sum _ {j=0}^{n-1} a_j (\omega_n^i)^j
y i = j = 0 ∑ n − 1 a j ( ω n i ) j
y i = ∑ j = 0 n − 1 a j ( ω n i ) j = ∑ j = 0 n 2 − 1 a 2 j ( ω n i ) 2 j + ∑ j = 0 n 2 − 1 a 2 j + 1 ( ω n i ) 2 j + 1 = ∑ j = 0 n 2 − 1 a 2 j ( ω n 2 i ) j + ω n i ∑ j = 0 n 2 − 1 a 2 j + 1 ( ω n 2 i ) j = F e v e n ( i ) + ω n i F o d d ( i )
\begin{align}
y_i = \sum _ {j=0}^{n-1} a_j (\omega_n^i)^j
&= \sum _ {j=0}^{\frac{n}{2}-1} a _ {2j} (\omega_n^i)^{2j} + \sum _ {j=0}^{\frac{n}{2}-1} a _ {2j+1} (\omega_n^i)^{2j+1}\newline
&= \sum _ {j=0}^{\frac{n}{2}-1} a _ {2j} (\omega _ {\frac{n}{2}}^i)^j + \omega_n^i \sum _ {j=0}^{\frac{n}{2}-1} a _ {2j+1} (\omega _ {\frac{n}{2}}^i)^j\newline
&= F _ {even}(i) + \omega_n^i F _ {odd}(i)
\end{align}
y i = j = 0 ∑ n − 1 a j ( ω n i ) j = j = 0 ∑ 2 n − 1 a 2 j ( ω n i ) 2 j + j = 0 ∑ 2 n − 1 a 2 j + 1 ( ω n i ) 2 j + 1 = j = 0 ∑ 2 n − 1 a 2 j ( ω 2 n i ) j + ω n i j = 0 ∑ 2 n − 1 a 2 j + 1 ( ω 2 n i ) j = F e v e n ( i ) + ω n i F o dd ( i ) 其中F e v e n , F o d d F _ {even}, F _ {odd} F e v e n , F o dd n 2 \frac{n}{2} 2 n i ∈ [ 0 , n 2 − 1 ] i \in [0, \frac{n}{2}-1] i ∈ [ 0 , 2 n − 1 ] F e v e n F _ {even} F e v e n F o d d F _ {odd} F o dd n 2 \frac{n}{2} 2 n
for 0 ≤ i < n 2 , { y i = F e v e n ( i ) + ω n i F o d d ( i ) y i + n 2 = F e v e n ( i ) − ω n i F o d d ( i )
\text{for } 0 \leq i < \frac{n}{2},
\left\{\begin{matrix}
y _ i &= F _ {even}(i) + \omega_n^i F _ {odd}(i) \newline
y _ {i+\frac{n}{2}} &= F _ {even}(i) - \omega_n^i F _ {odd}(i)
\end{matrix}\right.
for 0 ≤ i < 2 n , { y i y i + 2 n = F e v e n ( i ) + ω n i F o dd ( i ) = F e v e n ( i ) − ω n i F o dd ( i ) 時間複雜度有T ( n ) = 2 T ( n / 2 ) + O ( n ) T(n) = 2T(n/2) + \mathcal{O}(n) T ( n ) = 2 T ( n /2 ) + O ( n ) T ( n ) = O ( n log n ) T(n) = \mathcal{O}(n\log n) T ( n ) = O ( n log n )
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
const double PI = acos ( - 1 );
typedef complex < double > cd ;
vector < cd > FFT ( const vector < cd > & F ) { // assume F.size() == 2^k
if ( F . size () == 1 ) return F ; // base case (important)
vector < cd > rec [ 2 ], ans ;
for ( int i = 0 ; i < F . size (); i ++ ) rec [ i & 1 ]. push_back ( F [ i ]);
rec [ 0 ] = FFT ( rec [ 0 ]);
rec [ 1 ] = FFT ( rec [ 1 ]);
double theta = - 2 * PI / F . size ();
cd now = 1 , omega ( cos ( theta ), sin ( theta ));
ans . resize ( F . size ());
for ( int i = 0 ; i < F . size () / 2 ; i ++ ) {
ans [ i ] = rec [ 0 ] + now * rec [ 1 ];
ans [ i + F . size () / 2 ] = rec [ 0 ] - now * rec [ 1 ];
}
return ans ;
}
Inverse-FFT 那麼要怎麼做IFFT(傅立葉變換的逆變換)，也就是把點值表示轉換回係數呢？
[ 1 1 1 ⋯ 1 1 ω ω 2 ⋯ ω n − 1 ⋮ ⋮ ⋮ ⋱ ⋮ 1 ω n − 1 ( ω n − 1 ) 2 ⋯ ( ω n − 1 ) n − 1 ] [ c 0 c 1 ⋮ c n − 1 ] = [ C ( 1 ) C ( ω ) ⋮ C ( ω n − 1 ) ]
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1\newline
1 & \omega & \omega^2 & \cdots & \omega^{n-1}\newline
\vdots & \vdots & \vdots & \ddots & \vdots\newline
1 & \omega^{n-1} & (\omega^{n-1})^2 & \cdots & (\omega^{n-1})^{n-1}
\end{bmatrix}
\begin{bmatrix}
c_0\newline c_1\newline \vdots\newline c _ {n-1}
\end{bmatrix} =
\begin{bmatrix}
C(1)\newline C(\omega)\newline \vdots\newline C(\omega^{n-1})
\end{bmatrix}
1 1 ⋮ 1 1 ω ⋮ ω n − 1 1 ω 2 ⋮ ( ω n − 1 ) 2 ⋯ ⋯ ⋱ ⋯ 1 ω n − 1 ⋮ ( ω n − 1 ) n − 1 c 0 c 1 ⋮ c n − 1 = C ( 1 ) C ( ω ) ⋮ C ( ω n − 1 ) V = [ ω i j ] V = [\omega^{ij}] V = [ ω ij ] V ′ = [ 1 n ω − i j ] V' = [\frac{1}{n}\omega^{-ij}] V ′ = [ n 1 ω − ij ]
說明:
[ V ∗ V ′ ] i , j = ∑ k = 0 n − 1 V i , k V k , j ′ = 1 n ∑ k = 0 n − 1 ω k ( i − j )
[V * V'] _ {i,j} = \sum _ {k=0}^{n-1} V _ {i,k} V' _ {k,j} = \frac{1}{n}\sum _ {k=0}^{n-1} \omega^{k(i-j)}
[ V ∗ V ′ ] i , j = k = 0 ∑ n − 1 V i , k V k , j ′ = n 1 k = 0 ∑ n − 1 ω k ( i − j ) i = j i=j i = j i ≠ j i \neq j i = j 0 0 0
可以發現我們只需要把FFT的ω \omega ω n n n
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
const double PI = acos ( - 1 );
typedef complex < double > cd ;
vector < cd > FFT ( const vector < cd > & F , bool inv ) { // assume F.size() == 2^k
if ( F . size () == 1 ) return F ; // base case (important)
vector < cd > rec [ 2 ];
for ( int i = 0 ; i < F . size (); i ++ ) rec [ i & 1 ]. push_back ( F [ i ]);
rec [ 0 ] = FFT ( rec [ 0 ], inv );
rec [ 1 ] = FFT ( rec [ 1 ], inv );
double theta = ( inv ? 1 : - 1 ) * 2 * PI / F . size ();
cd now = 1 , omega ( cos ( theta ), sin ( theta ));
vector < cd > ans ( F . size ());
for ( int i = 0 ; i < F . size () / 2 ; i ++ ) {
ans [ i ] = rec [ 0 ][ i ] + now * rec [ 1 ][ i ];
ans [ i + F . size () / 2 ] = rec [ 0 ][ i ] - now * rec [ 1 ][ i ];
now *= omega ;
}
if ( inv ) for ( int i = 0 ; i < ans . size (); i ++ ) ans [ i ] /= 2 ;
return ans ;
}
Convolution 有了FFT和IFFT兩個工具，我們要做捲積就很簡單了
確定兩個多項式相乘的次數，並且選擇一個足夠大的n = 2 k n = 2^k n = 2 k 利用Cooley-Tukey演算法求出A , B A,B A , B A ^ , B ^ \hat A, \hat B A ^ , B ^ 將A ^ , B ^ \hat A, \hat B A ^ , B ^ C ^ \hat C C ^ 再利用Cooley-Tukey演算法求出C C C 1
2
3
4
5
6
7
8
9
10
11
12
13
vector < cd > A { 1 , 3 , 4 };
vector < cd > B { 1 , 2 , 5 };
signed main () {
int n = 1 << __lg ( A . size () + B . size ()) + 1 ;
A . resize ( n );
B . resize ( n );
A = FFT ( A , 0 );
B = FFT ( B , 0 );
vector < cd > C ( n );
for ( int i = 0 ; i < n ; i ++ ) C [ i ] = A [ i ] * B [ i ];
C = FFT ( C , 1 );
for ( int i = 0 ; i < n ; i ++ ) cout << C [ i ]. real () << " \n " [ i == n - 1 ];
}
Iterative Version 迭代的版本不但簡單執行時間又快，值得記一下
觀察遞迴的情況
100010 1 → 1 100010 111110 1 → 1 111110 101010 0 → 0 101010
\begin{gathered}
100010\textbf{1} \rightarrow \textbf{1}100010\\
111110\textbf{1} \rightarrow \textbf{1}111110\\
101010\textbf{0} \rightarrow \textbf{0}101010
\end{gathered}
100010 1 → 1 100010 111110 1 → 1 111110 101010 0 → 0 101010 k = log 2 n k = \log_2{n} k = log 2 n i i i j j j j j j i i i k k k
那要拿哪些合併呢？其實每個相鄰的兩塊的相同位置對應的就是F e v e n F _ {even} F e v e n F o d d F _ {odd} F o dd y i y_i y i
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
typedef complex < double > cd ;
void FFT ( cd F [], int n , bool inv ) { // in-place FFT, also assume n = 2^k
for ( int i = 0 , j = 0 ; i < n ; i ++ ) {
if ( i < j ) swap ( F [ i ], F [ j ]);
// magic! (maintain j to be the bit reverse of i)
for ( int k = n >> 1 ; ( j ^= k ) < k ; k >>= 1 );
}
for ( int step = 1 ; step < n ; step <<= 1 ) {
double theta = ( inv ? 1 : - 1 ) * PI / step ;
cd omega ( cos ( theta ), sin ( theta ));
for ( int i = 0 ; i < n ; i += step * 2 ) {
cd now ( 1 , 0 );
for ( int j = 0 ; j < step ; j ++ ) {
cd a = F [ i + j ];
cd b = F [ i + j + step ] * now ;
F [ i + j ] = a + b ;
F [ i + j + step ] = a - b ;
now *= omega ;
}
}
}
if ( inv ) for ( int i = 0 ; i < n ; i ++ ) F [ i ] /= n ;
}
NTT 注意到我們可以實行分治的關鍵就是存在一個ω \omega ω
ω n = 1 ∧ ∀ 0 ≤ i < j < n , ω i ≠ ω j
\omega ^ n = 1 \wedge
\forall 0 \leq i < j < n, \omega^i \neq \omega^j
ω n = 1 ∧ ∀0 ≤ i < j < n , ω i = ω j 現在我們想要在模一個質數p p p
費馬小定理表明
∀ ( a , p ) = 1 , a φ ( p ) ≡ 1 ( m o d p )
\forall (a,p) = 1, a^{\varphi(p)} \equiv 1 \pmod p
∀ ( a , p ) = 1 , a φ ( p ) ≡ 1 ( mod p ) 如果有原根g g g
∀ 0 ≤ i < j < φ ( p ) , g i ≢ g j ( m o d p )
\forall 0 \leq i < j < \varphi(p), g^i \not \equiv g^j \pmod p
∀0 ≤ i < j < φ ( p ) , g i ≡ g j ( mod p ) 那麼ω \omega ω ω n ≡ g φ ( p ) n \omega_n \equiv g ^ {\frac{\varphi(p)}{n}} ω n ≡ g n φ ( p ) ω n \omega_n ω n
這樣做必須滿足n ∣ φ ( p ) n | \varphi(p) n ∣ φ ( p ) n n n φ ( p ) = p − 1 = t ⋅ 2 k \varphi(p) = p-1 = t \cdot 2^k φ ( p ) = p − 1 = t ⋅ 2 k t t t p p p n n n 2 k 2^k 2 k p − 1 p-1 p − 1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
const int64_t MOD = 998244353 , G = 3 ;
int64_t modpow ( int64_t e , int64_t p , int64_t m ) {
int64_t r = 1 ;
for (; p ; p >>= 1 ) {
if ( p & 1 ) r = r * e % m ;
e = e * e % m ;
}
return r ;
}
void NTT ( int64_t F [], int n , bool inv ) { // assume n = 2^k!
for ( int i = 0 , j = 0 ; i < n ; i ++ ) {
if ( i < j ) swap ( F [ i ], F [ j ]);
for ( int k = n >> 1 ; ( j ^= k ) < k ; k >>= 1 );
}
for ( int step = 1 ; step < n ; step <<= 1 ) {
//may preprocess to boost
int64_t omega = modpow ( G , ( MOD - 1 ) / ( step * 2 ), MOD );
if ( inv ) omega = modpow ( omega , MOD - 2 , MOD );
for ( int i = 0 ; i < n ; i += step * 2 ) {
int64_t now = 1 ;
for ( int j = 0 ; j < step ; j ++ ) {
cd a = F [ i + j ];
cd b = F [ i + j + step ] * now % MOD ;
// reduce the use of % operator
F [ i + j ] = ( a + b < MOD ? a + b : a + b - MOD );
F [ i + j + step ] = ( a - b < 0 ? a - b + MOD : a - b );
now = now * omega % MOD ;
}
}
}
if ( inv ) {
int64_t invn = modpow ( n , MOD - 2 , MOD );
for ( int i = 0 ; i < n ; i ++ ) F [ i ] = F [ i ] * invn % MOD ;
}
}
→NTT模數表←
中國剩餘? 一個模數合不合適取決於最後答案的大小c c c n n n n c 2 nc^2 n c 2
如果真實的答案不是指數或階乘那種直接爆炸的數值
End FFT與NTT的利用其實滿少的，大部分不是大數乘法就是生成函數，以後有時間再放一篇講好了
