「多项式」学习笔记
多项式乘法
已知两个多项式的系数表达式，求其乘积的系数表达式。
FFT
系数表达式逐项相乘，复杂度\(O(n^2)\)，而点值表达式相乘复杂度为\(O(n)\)。因此我们要快速地将两个多项式转化为点值表达式，完成点值表达式的乘法，然后转为系数表达式得到结果。表达式转换的过程分别是离散傅里叶变换(\(DFT\))和离散傅里叶逆变换(\(IDFT\))。这两个过程统称为快速傅里叶变换(\(FFT\))。
离散傅里叶变换(DFT)
目标：快速求得系数表达式对应的点值表达式。
设有多项式
其中\(n\)为2的幂，不足则补0。我们按照其指数的奇偶将其分为两部分，变成两个新多项式。
则有\(A(x)=A_1(x^2)+xA_2(x^2)\)
单位根
在复平面内以原点为圆心，半径为1的圆称为单位圆。以\((1,0)\)为起点，将单位圆\(n\)等分，这些等分点是\(n\)个单位根。记第一个单位根为\(\omega_n\)，根据复数相乘模长相乘，幅角相加的原则，第\(k\)个单位根记作\(\omega_n^k\)。
根据三角函数推得单位根的表示方法：\(\omega_n^k=cos\frac{2k\pi}{n}+isin\frac{2k\pi}{n}\)
定理一（相消定理，折半引理）：\(\omega_{d*n}^{d*k}=\omega_n^k,\omega_n^k=\omega_{\frac{n}{2}}^{\frac{k}{2}}\)
定理二：\(\omega_n^k=-\omega_n^{k+\frac{n}{2}}\)
我们选择将单位根带入这两个子多项式求点值。
设\(k<\dfrac{n}{2}\)，将任意\(\omega_n^k\)代入\(A(x)\)可得$$A(\omega_nk)=A_1(\omega_{\frac{n}{2}})+\omega_nkA_2(\omega_{\frac{n}{2}}k)$$
将任意\(\omega_n^{k+\frac{n}{2}}\)代入\(A(x)\)可得$$A(\omega_n{k+\frac{n}{2}})=A_1(\omega_{\frac{n}{2}})-\omega_nkA_2(\omega_{\frac{n}{2}}k)$$
所以只要求出\(A_1(\omega_{\frac{n}{2}}^{k})\)和\(A_2(\omega_{\frac{n}{2}}^k)\)，就可以求得\(A(\omega_n^k)\)和\(A(\omega_n^{k+\frac{n}{2}})\)。也就是说，\(k\)只要取\([0,\frac{n}{2})\)就可以得出另一半\([\frac{n}{2},n)\)。因此就得到了整个区间\([0,n)\)，也就求出了\(A(x)\)的点值表达式。
递归求解即可，复杂度\(O(n \log n)\)
离散傅里叶逆变换(IDFT)
目标：将点值表达式转化回系数表达式
设点值表达式乘法结果为\((y_0,y_1,...,y_{n-1})\)。设它对应的系数表达式（答案）为\((a_0,a_1,...,a_{n-1})\)
我们构造一个以\((y_0,y_1,...,y_{n-1})\)为系数的多项式
此时多项式\(B(x)\)是系数表示的。我们依次将单位根的倒数（\(\omega_n^0,\omega_n^{-1},...,\omega_n^{1-n}\)）代入得到一个点值表达式\((z_0,z_1,...,z_{n-1})\)。
而因为后半部分\(\sum_{j=0}^{n-1}(\omega_n^{j-k})^i\)可由等比数列求和公式得：
注意，要使此式有意义，需满足分母不为0。而当\(j=k\)时分母为0。因此分类讨论：
当\(j \neq k\)时，\(\sum_{j=0}^{n-1}(\omega_n^{j-k})^i=0\)。
当\(j = k\)时，\(\sum_{j=0}^{n-1}(\omega_n^{j-k})^i=n\)
综上
所以
因此我们只需要把我们得到的点值表达式看做是系数表达式，再做一次\(DFT\)（其中单位根以倒数形式代入）就能够得到我们所要的结果了。最后还要除以\(n\)。
迭代实现
刚才这样是通过递归实现的，常数较大。
如果能事先确定叶节点的值，就可以直接向上合并了。之前系数是按照奇偶分左右的。第\(i\)层的分配看二进制的倒数第\(i\)位……因此系数的位置就是其初始位置的二进制倒序。
当前这一轮我们需要利用两个点值\(A_1(\omega_{\frac{n}{2}}^{k})\)和\(A_2(\omega_{\frac{n}{2}}^k)\)，我们利用他们造出了\(A(\omega_n^k)\)并推出了\(A(\omega_n^{k+\frac{n}{2}})\)，恰好可以把它们保存到原来的位置。这就是所谓的“蝴蝶操作”。
NTT
\(FFT\)需要用复数运算，存在精度问题。\(NTT\)（快速数论变换）是模意义下的\(FFT\)，可以避免精度问题，还快了不少。而其实现原理就是用模数的原根来代替单位根。
设模数\(P\)的原根为\(g\)，根据费马小定理我们知道\(g^{P-1} \equiv 1 (\text{mod} P)\)。\(n\)个原根是\(g^{\frac{P-1}{n}*k}\)。这就要求\(n\)是\(P-1\)的因数，而\(n\)又是\(2\)的幂。\(998244353\)就是这样一个满足条件的模数，它的原根是\(3\)。对于其他的模数，这里有{一张表}。
/*DennyQi 2019*/
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <queue>
using namespace std;
const int N = 4000010;
const int P = 998244353;
const int INF = 0x3f3f3f3f;
inline int Max(const int& a, const int& b){ return a>b?a:b; }
inline int Min(const int& a, const int& b){ return a<b?a:b; }
inline int mul(const int& a, const int& b){ return 1ll*a*b%P; }
inline int add(const int& a, const int& b){ return (a+b>=P)?a+b-P:a+b; }
inline int sub(const int& a, const int& b){ return (a-b<0)?a-b+P:a-b; }
inline int read(){
int x = 0, w = 0; char c = getchar();
while(c^'-' && (c<'0' || c>'9')) c = getchar();
if(c=='-') w = 1, c = getchar();
while(c>='0' && c<='9') x = x*10+c-'0', c = getchar(); return w?-x:x;
}
int n,m,lim=1,len,invn,a[N],b[N],r[N];
inline int qpow(int x, int y){
int res = 1;
while(y){
if(y & 1) res = mul(res,x);
y >>= 1, x = mul(x,x);
}
return res;
}
inline void NTT(int* a, int Tp){
for(int i = 1; i < lim; ++i) if(i < r[i]) swap(a[i],a[r[i]]);
for(int i = 1; i < lim; i <<= 1){
int w0 = Tp ? qpow(qpow(3,(P-1)/i/2),P-2) : qpow(3,(P-1)/i/2);
for(int j = 0; j < lim; j += (i<<1)){
int w = 1;
for(int k = j; k < j+i; ++k){
const int tmp = mul(w,a[k+i]);
a[k+i] = sub(a[k],tmp);
a[k] = add(a[k],tmp);
w = mul(w,w0);
}
}
}
}
int main(){
n = read(), m = read();
for(int i = 0; i <= n; ++i) a[i] = read();
for(int i = 0; i <= m; ++i) b[i] = read();
len = n+m+1;
while(lim <= len) lim <<= 1;
for(int i = 1; i < lim; ++i) r[i] = (r[i>>1]>>1)|((i&1)?(lim>>1):0);
NTT(a,0);
NTT(b,0);
for(int i = 0; i < lim; ++i) a[i] = mul(a[i],b[i]);
NTT(a,1);
invn = qpow(lim,P-2);
for(int i = 0; i < len; ++i) a[i] = mul(a[i],invn);
for(int i = 0; i < len; ++i) printf("%d ",a[i]);
printf("\n");
return 0;
}
