FFT
问题引入
示例问题
给定一个 \(n\) 次多项式 \(F(x)\) ，和一个 \(m\) 次多项式 \(G(x)\) 。
请求出 \(F(x)\) 和 \(G(x)\) 的乘积。
这里如果我们暴力去做就是 \(\mathcal O(n m)\) 的，但是貌似这里已经没有更好的办法了。
但是，如果题目提供的是 "两个多项式上相同 \(x\) 下的多个点" （这里被称作 "点表示法" ） 那么对于在两个多项式上的两个点 \((x, F(x)), (x, G(x))\) ，他们在乘积多项式上的点就直接是 \((x, F(x)\cdot G(x))\) ，只需要 \(\mathcal O(n+m+1)\) 就可以求解。
于是现在的问题就变成了需要把原来的多项式转化为 "点表示法" 然后合并之后快速切换回去。
Step 1： 考虑如何转化过去，我们发现此时就算我们随机找 \(n + m + 1\) 个点也需要 \(\mathcal O\left((n+m+1)^2\right)\) ，但是由于选择任意不同点都是一样的，所以可以考虑选择一些特殊点进行计算。
Step 2： 考虑如何转化回去，首先有一个结论，对于一个 \(n\) 次多项式，需要至少 \(n+1\) 个点就能确定它，既然我们通过 Step 1 已经求出这 \(n+1\) 个点了，那么一定可以做出来。而最直观的办法就是使用高斯消元，但是时间复杂度 \(\mathcal O(n^3)\) 。
所以现在两个步骤都是不能直接使用前面的办法解决的，这时候就需要使用 \(\texttt{FFT}\) 。
复数
关于 Step 1 我们不是需要找到 \(n+1\) 个特殊点吗，所以这里我们可以使用性质良好的 复数 作为特殊点。
当然对于这个多项式，他当然支持 复数 ， 并且 Step 2 的那个结论仍然存在。
复数运算
加法 \((a + b i) + (c+di) = (a+c) + (b+d)i\) ，在几何意义中就是和向量加法相同。
乘法 \((a+bi)\times (c+di) = (ac - bd) + (ad + bc)i\) ，在几何意义中就是复平面上角度相加，长度相乘。
单位根
我们在复平面上画一个单位圆，然后均分这个园到 \(n\) 份，然后每两份之间的交界点就是一个单位根，使用 \(w_n^k\) ，图示如下：
此时有以下性质：
\(w_n^k = \left( cos(\frac{k \cdot 2\pi}{n}), sin(\frac{k \cdot 2\pi}{n}) \right)\)
\(w_n^k = w_{2n}^{2k}\)
\(w_n^k = w_n^{k+n}\)
\(w_n^{k + n / 2} = -w_n^{k}\)
\(w_n^{a+b} = w_n^{a} * w_n^{b}\)
FFT
现在终于到真正的 FFT 了，这里我们令一个多项式 \(f(x) = \sum_{i=0}^{n-1} a_i x^i\) 。
我们考虑求出对于 \(f(w_n)，k \in [0, n)\) ，那么这里我们考虑分成两个子问题：
那么不难发现: \(f(x) = f_1(x^2) + x * f_2(x^2)\) 。
那么把特殊点 \(w_n\) 带入：
\[
f(w_n^k) = f_1((w_n^k)^2) + w_n^k \times f_2((w_n^k)^2)
\]
此时分为两种情况：
\(0\le k < \frac{n}{2}\) : \(f(w_n^k) = f_1(w_n^{2k}) + w_n^k \times f_2(w_n^{2k}) = f_1(w_{n/2}^k) + w_n^k \times f_2(w_{n/2}^{k})\)
\(\frac{n}{2} \le \frac{n}{2} + k < n\) : \(f(w_n^{k+\frac{n}{2}}) = f_1(w_n^{2k+n}) + w_n^{k+\frac{n}{2}} \times f_2(w_n^{2k+n}) = f_1(w_{n/2}^k) - w_n^k \times f_2(w_{n/2}^{k})\)
所以我们发现这两个式子居然只有中间的符号不一样, 于是现在两个问题被分成了来嗯个规模相同的子问题。只需要递归下去。
然后当问题规模为 \(1\) 的时候直接结束，总共时间复杂度 \(\mathcal O(n \log n)\) 。
IFFT
这里我们已经得到了 \(c = (c_0, c_1, c_2, /dots, c_{n-1})\) ，然后我们考虑对于 \(f(x) = \sum_{i=1}^{n-1} c_i x^i\) 求解他的 FFT : \(d = (d_0, d_1, d_2, /dots, d_{n-1})\) 。
所以可以发现:
\[
d_k = f(W_n^{-k})
\]
\[
= \sum_{i=0}^{n-1} c_i(W_n^{-k})^i
\]
\[
= \sum_{i=0}^{n-1} c_i(W_n^{-k})^i
\]
\[
= \sum_{i=0}^{n-1} \sum_{j=0}^{n-1} a_j (W_n^{i})^j (W_n^{-k})^i
\]
\[
= \sum_{j=0}^{n-1} a_j \sum_{i=0}^{n-1} (W_n^{j-k})^i
\]
\[
= \sum_{j=0}^{n-1} a_j \frac{1 - (W_n^{j-k})^{n}}{1-W_n^{j-k}}
\]
\[
= \sum_{j=0}^{n-1} a_j \frac{1 - (W_n^n)^{j-k}}{1-W_n^{j-k}}
\]
所以只有 \(j = k\) 时，所以 \(d_i = n \cdot a_i\) ， 那么 \(a_i = \frac{d_i}n\) 。
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
36
37
38
39 void FFT ( int n , vector < com >& v , bool invert ) {
if ( n == 1 ) return ;
int mid = n >> 1 ;
vector < com > al ( mid + 1 ), ar ( mid + 1 );
for ( int i = 0 ; i < mid ; i ++ )
al [ i ] = v [ i << 1 ], ar [ i ] = v [ i << 1 | 1 ];
FFT ( mid , al , invert ), FFT ( mid , ar , invert );
com j ( 1 , 0 );
com step ( cos ( 2 * M_PI / n ), sin ( 2 * M_PI / n * ( invert ? -1 : 1 )));
for ( int i = 0 ; i < mid ; i ++ , j = j * step ) {
v [ i ] = al [ i ] + j * ar [ i ];
v [ i + mid ] = al [ i ] - j * ar [ i ];
}
}
/*~~~~~~~~~~~~~~~~~~~~ Boundary Line ~~~~~~~~~~~~~~~~~~~~*/
signed main () {
cin >> n >> m ; n ++ , m ++ ;
a . resize ( n + 1 ), b . resize ( m + 1 );
for ( int i = 0 ; i < n ; i ++ ) cin >> a [ i ];
for ( int i = 0 ; i < m ; i ++ ) cin >> b [ i ];
int mxn = n + m -1 , t = 1 ;
while ( t <= mxn ) t <<= 1 ;
a . resize ( t + 1 ), b . resize ( t + 1 );
c . resize ( t + 1 );
FFT ( t , a , 0 ), FFT ( t , b , 0 );
for ( int i = 0 ; i < t ; i ++ ) c [ i ] = a [ i ] * b [ i ];
FFT ( t , c , 1 );
for ( int i = 0 ; i < mxn ; i ++ ) cout << ( int ) round ( c [ i ]. real () / t ) << ' ' ;
return 0 ;
}
NTT
首先我们发现 FFT 有一个巨大的问题，其不能取膜，如果题目就是要求你取膜应当如何是好。
所以这里我们引入几个概念：
价: 在模 \(p\) 下，如果 \(a^x \equiv 1 \mod p\) ， 那么称 \(x\) 为 \(a\) 在 \(p\) 下的价，记作 \(|a|\) 。
原根 \(g\) : 对于 \(|g| = \varphi(p)\) ， 那么称 \(g\) 为 \(p\) 下的原根。对于 \(p\) 为质数的情况，其一定有原根，其原根满足 \(g^{p-1} \equiv 1 \mod p\) 。
然后我们现在实际上就是要选择一个单位根，满足前面的那些性质，这里我们选择单位根为 \(g^{\frac{p-1}{n}}\) ，那么此时 \(w_n^k = \left(g^{\frac{p-1}{n}}\right)^k\) 。
此时我们验证以下各个性质: 验证个球 。
然后代码:
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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51 emplate < int mod >
struct Pol {
typedef long long ll ;
static constexpr int G = 3 ;
int ksm ( int x , int y ) {
int ans = 1 ;
while ( y ) {
if ( y & 1 ) ans = ( ll ) ans * x % mod ;
x = ( ll ) x * x % mod , y >>= 1 ;
}
return ans ;
}
int rev [ N ];
void init ( int n ) {
int bit = 0 ;
while (( 1 << bit ) < n ) bit ++ ;
for ( int i = 0 ; i < n ; i ++ ) rev [ i ] = ( rev [ i >> 1 ] >> 1 ) + (( i & 1 ) << ( bit -1 ));
}
void NTT ( int n , vector < int >& v , bool invert ) {
static int invG = ksm ( G , mod -2 );
for ( int i = 0 ; i < n ; i ++ ) if ( i < rev [ i ]) swap ( v [ i ], v [ rev [ i ]]);
for ( int len = 1 ; len <= ( n >> 1 ); len <<= 1 ) {
int W = ksm ( invert ? invG : G , ( mod -1 ) / ( len << 1 ));
for ( int i = 0 ; i < n ; i += ( len << 1 )) {
int w = 1 ;
for ( int j = 0 ; j < len ; j ++ , w = 1L L * w * W % mod ) {
int a1 = v [ i + j ], a2 = 1L L * w * v [ i + j + len ] % mod ;
v [ i + j ] = ( a1 + a2 ) % mod , v [ i + j + len ] = (( a1 - a2 ) % mod + mod ) % mod ;
}
}
}
}
vector < int > main ( vector < int > a , vector < int > b ) {
int n = a . size (), m = b . size (), t = 1 ;
while ( t <= n + m -1 ) t <<= 1 ;
a . resize ( t ), b . resize ( t ), init ( t );
NTT ( t , a , 0 ), NTT ( t , b , 0 );
for ( int i = 0 ; i < t ; i ++ ) a [ i ] = ( ll ) a [ i ] * b [ i ] % mod ;
NTT ( t , a , 1 );
int invT = ksm ( t , mod -2 );
for ( int i = 0 ; i < n + m -1 ; i ++ ) a [ i ] = (( ll ) a [ i ] * invT % mod + mod ) % mod ;
return a ;
}
};
多项式其他运算
多项式求逆元
这里我们另 \(F(x)\) 是原多项式， \(G(x)\) 为新的多项式，直接推式子：
\[
F(x) \times G(x) \equiv 1 \pmod{x^n}
\]
后面就不写同余了。
\[
F(x) - \frac{1}{G(x)} = 0
\]
因此我们直接上牛顿迭代法：
\[
G(x) = G_0(x) - \frac{F(x) - \frac{1}{G_0(x)}}{\frac{1}{G_0^2(x)}}
\]
（这里 \(G_0\) 表示牛顿迭代法上一个，在 \(\bmod~ {x^{\lceil \frac{n}{2} \rceil}}\) 下，这是因为在多项式中每牛顿迭代法一次精度翻倍）
那么:
\[
G(x) = 2 G_0(x) - F(x) G_0^2(x)
\]
\[
= G_0(x)\times\left( 2 - F(x) G_0(x)\right)
\]
然后递归求解就可以了，代码：
1
2
3
4
5
6
7
8
9
10
11
12 void inv ( int n , const Poly & a , Poly & res ) {
if ( n == 1 ) return res . clean ( 1 ), res [ 0 ] = ksm ( a [ 0 ], mod -2 ), void ();
inv (( n + 1 ) >> 1 , a , res ), res . resize ( n );
static Poly b ; b . clean ();
int t = 1 ; while ( t <= ( n << 1 )) t <<= 1 ;
init ( t ), b = a , b . resize ( n );
NTT ( t , b , 0 ), NTT ( t , res , 0 );
for ( int i = 0 ; i < t ; i ++ ) res [ i ] = 1L L * norm ( 2 - 1L L * b [ i ] * res [ i ]) * res [ i ] % mod ;
NTT ( t , res , 1 ), res . resize ( n );
}
多项式对数
\(F, G\) 和上面相同，同样退式子：
\[G(x) = \ln F(x)\]
\[G'(x) = \frac{1}{F(x)} \cdot F'(x)\]
\[G'(x) = F'(x) \cdot F^{-1}(x)\]
然后直接使用 求导 ， 积分 ， 求逆元 就可以了。代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 // 多项式求积分
void Int ( const Poly & a , Poly & res ) {
res . resize ( a . size () + 1 );
for ( int i = 1 ; i <= a . size (); i ++ ) res [ i ] = 1L L * a [ i -1 ] * ksm ( i , mod -2 ) % mod ;
res [ 0 ] = 0 ;
}
// 多项式求导
void Der ( const Poly & a , Poly & res ) {
res . resize ( a . size () -1 );
for ( int i = 1 ; i < a . size (); i ++ ) res [ i -1 ] = 1L L * i * a [ i ] % mod ;
}
// 多项式 ln
void Ln ( const Poly & x , Poly & res ) {
static Poly a , b ; a . clean (), b . clean ();
Der ( x , a ), inv ( x . size (), x , b );
mul ( a , b , a ), Int ( a , res );
res . resize ( x . size ());
}
多项式指数
同样直接推式子：
\[G(x) = \exp(F(x))\]
\[F(x) = \ln G(x)\]
\[\ln G(X) - F(x) = 0\]
然后直接使用牛顿迭代法：
\[
G(x) = G_0(x) - \frac{\ln G_0(x) - F(x)}{\frac{1}{G_0(x)}}
\]
\[
G(x) = G_0(x) (1 - \ln G(x) + F(x))
\]
真好，前面写的 \(\ln\) 有可以用上了阿
代码：
1
2
3
4
5
6
7
8
9
10
11
12 void Exp ( int n , const Poly & x , Poly & res ) {
if ( n == 1 ) return res . clean ( 1 ), res [ 0 ] = 1 , void ();
Exp (( n + 1 ) >> 1 , x , res ), res . resize ( n );
static Poly lng , f ; Ln ( res , lng );
int t = 1 ; while ( t <= ( n << 1 )) t <<= 1 ;
init ( t ), f = x , f . resize ( n );
NTT ( t , f , 0 ), NTT ( t , lng , 0 ), NTT ( t , res , 0 );
for ( int i = 0 ; i < t ; i ++ ) res [ i ] = 1L L * norm ( 1 - lng [ i ] + f [ i ]) * res [ i ] % mod ;
NTT ( t , res , 1 ), res . resize ( n );
}
公用模板
参考代码
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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115 struct Poly {
#define norm(x) (((x)%mod+mod)%mod)
typedef long long ll ;
static constexpr int mod = 998244353 , G = 3 ;
static constexpr int N = 5e5 + 5 ;
vector < int > a ;
inline void resize ( int n ) { a . resize ( n ); }
inline int operator [] ( int x ) const { return a [ x ]; }
inline int & operator [] ( int x ) { return a [ x ]; }
inline int size () { return a . size (); }
inline int size () const { return a . size (); }
inline void clean () { a . clear (); }
inline void clean ( int m ) { a . clear (), a . resize ( m ); }
static int rev [ N ];
static inline void init ( int n ) {
int bit = 1 ;
while (( 1 << bit ) < n ) bit ++ ;
for ( int i = 1 ; i < n ; i ++ ) rev [ i ] = ( rev [ i >> 1 ] >> 1 ) | (( i & 1 ) << ( bit -1 ));
}
static int ksm ( int x , int y ) {
int ans = 1 ;
while ( y ) {
if ( y & 1 ) ans = 1L L * ans * x % mod ;
x = 1L L * x * x % mod , y >>= 1 ;
}
return ans ;
}
static void NTT ( int n , Poly & v , bool invert ) {
v . resize ( n );
for ( int i = 0 ; i < n ; i ++ ) if ( i < rev [ i ]) swap ( v [ i ], v [ rev [ i ]]);
static int invG = ksm ( G , mod -2 );
for ( int len = 1 ; len <= ( n >> 1 ); len <<= 1 ) {
int W = ksm ( invert ? invG : G , ( mod -1 ) / ( len << 1 ));
for ( int i = 0 ; i < n ; i += ( len << 1 )) {
int w = 1 ;
for ( int j = 0 ; j < len ; j ++ , w = 1L L * w * W % mod ) {
int a0 = v [ i + j ], a1 = 1L L * w * v [ i + j + len ] % mod ;
v [ i + j ] = ( a0 + a1 ) % mod , v [ i + j + len ] = norm ( a0 - a1 );
}
}
}
if ( invert ) {
int invN = ksm ( n , mod -2 ); // 这里不能加 static
for ( int i = 0 ; i < n ; i ++ ) v [ i ] = 1L L * v [ i ] * invN % mod ;
}
}
static void mul ( const Poly & va , const Poly & vb , Poly & res ) {
static Poly a , b ; a = va , b = vb ;
int t = 1 ; while ( t <= a . size () + b . size () - 1 ) t <<= 1 ;
init ( t ), res . clean (), res . resize ( t );
NTT ( t , a , 0 ), NTT ( t , b , 0 );
for ( int i = 0 ; i < t ; i ++ ) res [ i ] = 1L L * a [ i ] * b [ i ] % mod ;
NTT ( t , res , 1 ), res . resize ( a . size () + b . size () - 1 );
}
// 多项式逆元
static void inv ( int n , const Poly & a , Poly & res ) {
if ( n == 1 ) return res . clean ( 1 ), res [ 0 ] = ksm ( a [ 0 ], mod -2 ), void ();
inv (( n + 1 ) >> 1 , a , res ), res . resize ( n );
static Poly b ; b . clean ();
int t = 1 ; while ( t <= ( n << 1 )) t <<= 1 ;
init ( t ), b = a , b . resize ( n );
NTT ( t , b , 0 ), NTT ( t , res , 0 );
for ( int i = 0 ; i < t ; i ++ ) res [ i ] = 1L L * norm ( 2 - 1L L * b [ i ] * res [ i ]) * res [ i ] % mod ;
NTT ( t , res , 1 ), res . resize ( n );
}
// 多项式求积分
static void Int ( const Poly & a , Poly & res ) {
res . resize ( a . size () + 1 );
for ( int i = 1 ; i <= a . size (); i ++ ) res [ i ] = 1L L * a [ i -1 ] * ksm ( i , mod -2 ) % mod ;
res [ 0 ] = 0 ;
}
// 多项式求导
static void Der ( const Poly & a , Poly & res ) {
res . resize ( a . size () -1 );
for ( int i = 1 ; i < a . size (); i ++ ) res [ i -1 ] = 1L L * i * a [ i ] % mod ;
}
// 多项式 ln
static void Ln ( const Poly & x , Poly & res ) {
static Poly a , b ; a . clean (), b . clean ();
Der ( x , a ), inv ( x . size (), x , b );
mul ( a , b , a ), Int ( a , res );
}
// 多项式 exp
static void Exp ( int n , const Poly & x , Poly & res ) {
if ( n == 1 ) return res . clean ( 1 ), res [ 0 ] = 1 , void ();
Exp (( n + 1 ) >> 1 , x , res ), res . resize ( n );
static Poly lng , f ; Ln ( res , lng );
int t = 1 ; while ( t <= ( n << 1 )) t <<= 1 ;
init ( t ), f = x , f . resize ( n );
NTT ( t , f , 0 ), NTT ( t , lng , 0 ), NTT ( t , res , 0 );
for ( int i = 0 ; i < t ; i ++ ) res [ i ] = 1L L * norm ( 1 - lng [ i ] + f [ i ]) * res [ i ] % mod ;
NTT ( t , res , 1 ), res . resize ( n );
}
};
int Poly :: rev [ Poly :: N ];
