摘要 操纵 有限项/无限项 的多项式是 OI 数学中，尤其是生成函数中的重要内容。
以 快速傅里叶变换 为基石的多项式算法赋予了算法竞赛选手直接操纵生成函数的能力。
摘自 OI-Wiki
多项式的核心操作为多项式乘法，以此为延伸，可以进行一系列运算。例如：多项式求逆，多项式开根，多项式除法，多项式对数函数，多项式指数函数，多项式牛顿迭代，多项式幂函数等。
下面逐一剖析这些操作。
约定 多项式 f ( x ) f(x) f ( x ) f ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n x n f(x)=a_{0}+a_{1}x+a_{2}x^{2}+\cdots+a_{n}x^{n} f ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n x n { a i } \{a_i\} { a i } 系数表示法 。
y → f ( x ) y\rightarrow f(x) y → f ( x ) y y y f ( x ) f(x) f ( x )
由拉格朗日插值定理得到：使用 n + 1 n+1 n + 1 ( x i , y i ) (x_i,y_i) ( x i , y i ) x i x_i x i x i → f ( x ) = y i x_i\rightarrow f(x)=y_i x i → f ( x ) = y i n n n n + 1 n+1 n + 1 点值表示法 。
多项式乘法 朴素算法 多项式的乘法操作也被称为加法卷积操作。对给定多项式 f ( x ) f(x) f ( x ) g ( x ) g(x) g ( x )
f ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n x n g ( x ) = b 0 + b 1 x + b 2 x 2 + ⋯ + b m x m f(x)=a_{0}+a_{1}x+a_{2}x^{2}+\cdots+a_{n}x^{n}\\ g(x)=b_{0}+b_{1}x+b_{2}x^{2}+\cdots+b_{m}x^{m}\\ f ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n x n g ( x ) = b 0 + b 1 x + b 2 x 2 + ⋯ + b m x m
其乘法结果为
∑ i = 0 n ∑ j = 0 m a i b j x i + j \sum_{i=0}^{n}\sum_{j=0}^{m}a_{i}b_{j}x^{i+j} i = 0 ∑ n j = 0 ∑ m a i b j x i + j
有时也写作
∑ i + j = k a i b j x k \sum_{i+j=k}a_ib_jx^{k} i + j = k ∑ a i b j x k
这也是为什么其被称为加法卷积。
直接根据定义计算复杂度是 O ( n m ) O(nm) O ( nm )
单位根 如果多项式环上存在 2 n 2^{n} 2 n O ( n 2 n ) O(n2^{n}) O ( n 2 n ) 2 n 2^{n} 2 n
典型的两个例子是快速傅里叶变换和快速数论变换。复数域存在任意 2 n 2^{n} 2 n
定义 n n n ω n i , 0 ≤ i < n \omega_{n}^{i},0\le i<n ω n i , 0 ≤ i < n n n n
ω n a ω n b = ω n ( a + b ) m o d n \omega_{n}^{a}\omega_{n}^{b}=\omega_{n}^{(a+b)\bmod n} ω n a ω n b = ω n ( a + b ) mod n
( ω n a ) b = ω n a b m o d n (\omega_{n}^{a})^{b}=\omega_{n}^{ab\bmod n} ( ω n a ) b = ω n ab mod n
对于 2 ∣ n 2\vert n 2∣ n
ω n 2 i = ω n 2 i \omega_{n}^{2i}=\omega_{\frac{n}{2}}^{i} ω n 2 i = ω 2 n i
ω n i + n 2 = − ω n i \omega_{n}^{i+\frac{n}{2}}=-\omega_{n}^{i} ω n i + 2 n = − ω n i
下面称 ω n i \omega_{n}^{i} ω n i i i i n n n i i i
单位根优化多项式乘法 如果将 n n n n n n O ( n ) O(n) O ( n )
也就是说，使用单位根进行多项式乘法有三步：
将单位根代入多项式。
将点值相乘。
用点值乘积求系数表示。
而单位根的性质使得，对于 2 2 2 n n n O ( n log n ) O(n\log n) O ( n log n ) n n n n − 1 n-1 n − 1
第一步 对于多项式 f ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n − 1 x n − 1 f(x)=a_0+a_1x+a_2x^{2}+\cdots+a_{n-1}x^{n-1} f ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n − 1 x n − 1 n n n h ( x ) = a 0 + a 2 x + a 4 x 2 + ⋯ a n − 2 x n − 2 2 h(x)=a_0+a_2x+a_4x^{2}+\cdots a_{n-2}x^{\frac{n-2}{2}} h ( x ) = a 0 + a 2 x + a 4 x 2 + ⋯ a n − 2 x 2 n − 2 g ( x ) = a 1 + a 3 x + a 5 x 2 + ⋯ a n − 1 x n − 2 2 g(x)=a_1+a_3x+a_5x^{2}+\cdots a_{n-1}x^{\frac{n-2}{2}} g ( x ) = a 1 + a 3 x + a 5 x 2 + ⋯ a n − 1 x 2 n − 2
根据定义可以自然得到：
f ( x ) = h ( x 2 ) = x g ( x 2 ) \begin{equation} f(x)=h(x^2)=xg(x^2) \end{equation} f ( x ) = h ( x 2 ) = xg ( x 2 )
换句话说，我们提取 f ( x ) f(x) f ( x ) n 2 \frac{n}{2} 2 n
求 ω n i → f ( x ) \omega_{n}^{i}\rightarrow f(x) ω n i → f ( x ) i < n 2 i<\frac{n}{2} i < 2 n
ω n i → f ( x ) = ( ω n i ) 2 → h ( x ) + ω n i [ ( ω n i ) 2 → g ( x ) ] \omega_{n}^{i}\rightarrow f(x)=(\omega_{n}^{i})^2\rightarrow h(x)+\omega_{n}^{i}[(\omega_{n}^{i})^2\rightarrow g(x)] ω n i → f ( x ) = ( ω n i ) 2 → h ( x ) + ω n i [( ω n i ) 2 → g ( x )]
根据单位根性质进行变换，得到 ( ω n i ) 2 = ω n 2 i = ω n 2 i (\omega_{n}^{i})^2=\omega_{n}^{2i}=\omega_{\frac{n}{2}}^{i} ( ω n i ) 2 = ω n 2 i = ω 2 n i
ω n i → f ( x ) = ω n 2 i → h ( x ) + ω n i ( ω n 2 i → g ( x ) ) \begin{equation} \omega_{n}^{i}\rightarrow f(x)=\omega_{\frac{n}{2}}^{i}\rightarrow h(x)+\omega_{n}^{i}(\omega_{\frac{n}{2}}^{i}\rightarrow g(x)) \end{equation} ω n i → f ( x ) = ω 2 n i → h ( x ) + ω n i ( ω 2 n i → g ( x ))
求 ω n i + n 2 → f ( x ) \omega_{n}^{i+\frac{n}{2}}\rightarrow f(x) ω n i + 2 n → f ( x ) i < n 2 i<\frac{n}{2} i < 2 n
ω n i + n 2 → f ( x ) = ( ω n i + n 2 ) 2 → h ( x ) + ω n i + n 2 [ ( ω n i + n 2 ) 2 → g ( x ) ] \omega_{n}^{i+\frac{n}{2}}\rightarrow f(x)=(\omega_{n}^{i+\frac{n}{2}})^2\rightarrow h(x)+\omega_{n}^{i+\frac{n}{2}}[(\omega_{n}^{i+\frac{n}{2}})^2\rightarrow g(x)] ω n i + 2 n → f ( x ) = ( ω n i + 2 n ) 2 → h ( x ) + ω n i + 2 n [( ω n i + 2 n ) 2 → g ( x )]
根据单位根性质进行变换，得到 ( ω n i + n 2 ) 2 = ω n 2 i + n = ω n 2 i (\omega_{n}^{i+\frac{n}{2}})^2=\omega_{n}^{2i+n}=\omega_{\frac{n}{2}}^{i} ( ω n i + 2 n ) 2 = ω n 2 i + n = ω 2 n i ω n i + n 2 = − ω n i \omega_{n}^{i+\frac{n}{2}}=-\omega_{n}^{i} ω n i + 2 n = − ω n i
ω n i + n 2 → f ( x ) = ω n 2 i → h ( x ) − ω n i ( ω n 2 i → g ( x ) ) \begin{equation} \omega_{n}^{i+\frac{n}{2}}\rightarrow f(x)=\omega_{\frac{n}{2}}^{i}\rightarrow h(x)-\omega_{n}^{i}(\omega_{\frac{n}{2}}^{i}\rightarrow g(x)) \end{equation} ω n i + 2 n → f ( x ) = ω 2 n i → h ( x ) − ω n i ( ω 2 n i → g ( x ))
观察 (2) 式与 (3) 式，可以发现，两式只有一个符号的区别。也就是说，只要知道 ω n 2 i → h ( x ) \omega_{\frac{n}{2}}^{i}\rightarrow h(x) ω 2 n i → h ( x ) ω n 2 i → g ( x ) \omega_{\frac{n}{2}}^{i}\rightarrow g(x) ω 2 n i → g ( x ) O ( n ) O(n) O ( n ) ω n i → f ( x ) \omega_{n}^{i}\rightarrow f(x) ω n i → f ( x )
至此我们在 O ( n log n ) O(n\log n) O ( n log n )
第三步 记 f i = ω n i → f ( x ) f_i=\omega_{n}^{i}\rightarrow f(x) f i = ω n i → f ( x ) g ( x ) = f 0 + f 1 x + f 2 x 2 + ⋯ + f n − 1 x n − 1 g(x)=f_0+f_1x+f_2x^{2}+\cdots+f_{n-1}x^{n-1} g ( x ) = f 0 + f 1 x + f 2 x 2 + ⋯ + f n − 1 x n − 1 f ( x ) f(x) f ( x ) a i = ω n − i → g ( x ) n a_i=\dfrac{\omega_{n}^{-i}\rightarrow g(x)}{n} a i = n ω n − i → g ( x )
证明 由 f i f_i f i
f i = a 0 + a 1 ω n i + a 2 ω n 2 i + ⋯ ω n i ( n − 1 ) f_i=a_0+a_1\omega_{n}^{i}+a_2\omega_{n}^{2i}+\cdots\omega_{n}^{i(n-1)} f i = a 0 + a 1 ω n i + a 2 ω n 2 i + ⋯ ω n i ( n − 1 )
ω n − j → g ( x ) \omega_{n}^{-j}\rightarrow g(x) ω n − j → g ( x ) f i x i f_ix^{i} f i x i ω − i j f i \omega^{-ij}f_i ω − ij f i
展开 f i f_i f i
a 0 ω i ( 0 − j ) + a 1 ω n i ( 1 − j ) + a 2 ω n i ( 2 − j ) + ⋯ ω n i ( n − 1 − j ) a_0\omega^{i(0-j)}+a_1\omega_{n}^{i(1-j)}+a_2\omega_{n}^{i(2-j)}+\cdots\omega_{n}^{i(n-1-j)} a 0 ω i ( 0 − j ) + a 1 ω n i ( 1 − j ) + a 2 ω n i ( 2 − j ) + ⋯ ω n i ( n − 1 − j )
对于其中一项 a k a_k a k a k ω i ( k − j ) a_k\omega^{i(k-j)} a k ω i ( k − j )
当 k = j k=j k = j n a k na_{k} n a k
当 k ≠ j k\not=j k = j i ( k − j ) i(k-j) i ( k − j ) n n n [ 0 , n ) [0,n) [ 0 , n ) ∑ i = 0 n − 1 ω n i a k \sum_{i=0}^{n-1}\omega_{n}^{i}a_k ∑ i = 0 n − 1 ω n i a k a k ( 1 − ω n n ) 1 − ω n 1 = 0 \dfrac{a_k(1-\omega_n^{n})}{1-\omega_{n}^{1}}=0 1 − ω n 1 a k ( 1 − ω n n ) = 0
故对于 a k a_k a k ω n j , j = k \omega_{n}^{j},j=k ω n j , j = k n a k na_k n a k 0 0 0
所以再做一次类似第一步的操作即可还原系数。
快速傅里叶变换 对于 2 n 2^{n} 2 n f ( x ) f(x) f ( x ) g ( x ) g(x) g ( x ) O ( n 2 n ) O(n2^{n}) O ( n 2 n )
在快速傅里叶变换中，取 ω n i = cos ( i n π ) + sin ( i n π ) ı \omega_{n}^{i}=\cos(\frac{i}{n}\pi)+\sin(\frac{i}{n}\pi)\imath ω n i = cos ( n i π ) + sin ( n i π ) ı \imath
将第 i i i n n n
如此我们得到了快速傅里叶变换算法。
实现上，可以直接分治，但是常数较大，且不容易实现。推荐使用倍增法进行实现。
具体来说，先使用被称为 蝶形运算优化 的技巧来合并减小不必要的内存拷贝。观察系数 a i a_i a i 1 1 1 i i i 0 0 0 1 1 1 0 0 0 0 0 0
设 f[i] 表示 i 的二进制翻转后的结果。则可以递推得到，f[i]=f[i>>1]|((i&1)?(1<<s):0)，其中 s s s n n n
下面给出 洛谷模板题 的代码。
Sample Code(c++) c++ #include <bits/stdc++.h> using namespace std ; const int N = 1 << 21 ; const int mod = 998244353 , RT = 3 ; using ll = long long ; ll qpow ( ll a , ll b ) { ll res = 1 ; for (; b; b >>= 1 , a = a * a % mod) if (b & 1 ) res = res * a % mod; return res; } void ntt ( int * a , int n , bool inv ) { static int G[N], invG[N], rev[N]; static int lst_n = 0 ; if (n != lst_n) { lst_n = n; for ( int i = 0 ; i < n; ++ i) rev[i] = (rev[i >> 1 ] >> 1 ) | ((i & 1 ) ? (n >> 1 ) : 0 ); for ( int i = 1 ; i < n; i <<= 1 ) { int g1 = qpow (RT, (mod - 1 ) / (i << 1 )), ig1 = qpow (g1, mod - 2 ); int g = 1 , ig = 1 ; for ( int j = i; j < i + i; ++ j) { G[j] = g, invG[j] = ig; g = 1 ull * g * g1 % mod, ig = 1 ull * ig * ig1 % mod; } } } for ( int i = 0 ; i < n; ++ i) if (rev[i] < i) swap (a[rev[i]], a[i]); for ( int i = 1 , x; i < n; i <<= 1 ) { for ( int * j = a; j < a + n; j += (i << 1 )) { for ( int * k = j, * buf = (inv ? invG : G) + i; k < j + i; ++ k, ++ buf) { x = 1 ull * k[i] * * buf % mod; if ((k[i] = * k + mod - x) >= mod) k[i] -= mod; if (( * k += x) >= mod) * k -= mod; } } } if (inv) { int invn = qpow (n, mod - 2 ); for ( int i = 0 ; i < n; ++ i) a[i] = 1 ull * a[i] * invn % mod; } } void mul ( int * a , int * b , int l1 , int l2 ) { int n = 1 ; while (n <= l1 + l2) n <<= 1 ; ntt (a, n, false ), ntt (b, n, false ); for ( int i = 0 ; i < n; ++ i) a[i] = 1 ull * a[i] * b[i] % mod; ntt (a, n, true ); } int a[N], b[N]; int main () { cin. tie ( 0 )-> sync_with_stdio ( false ); int n, m; cin >> n >> m; for ( int i = 0 ; i <= n; ++ i) cin >> a[i]; for ( int i = 0 ; i <= m; ++ i) cin >> b[i]; mul (a, b, n, m); for ( int i = 0 ; i <= n + m; ++ i) cout << a[i] << ' ' ; cout << ' \n ' ; } 快速数论变换 由于快速傅里叶变换使用了复数，不可避免的引起了精度误差问题，并且实数运算是低效的，我们使用原根的一些性质来替代单位根，使得运算更为准确快捷。
参见：原根与快速数论变换