「多项式」学习笔记

多项式乘法

已知两个多项式的系数表达式，求其乘积的系数表达式。

[c_n = sumlimits_{i=0}^{n}a_ib_{n-i} ]

FFT

系数表达式逐项相乘，复杂度(O(n^2))，而点值表达式相乘复杂度为(O(n))。因此我们要快速地将两个多项式转化为点值表达式，完成点值表达式的乘法，然后转为系数表达式得到结果。表达式转换的过程分别是离散傅里叶变换((DFT))和离散傅里叶逆变换((IDFT))。这两个过程统称为快速傅里叶变换((FFT))。

离散傅里叶变换(DFT)

目标：快速求得系数表达式对应的点值表达式。

设有多项式

[A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1} ]

其中(n)为2的幂，不足则补0。我们按照其指数的奇偶将其分为两部分，变成两个新多项式。

[A_1(x)=a_0+a_2x+...+a_{n-2}x^{frac{n}{2}-1} ]

[A_2(x)=a_1+a_3x+...+a_{n-1}x^{frac{n}{2}-1} ]

则有(A(x)=A_1(x^2)+xA_2(x^2))

单位根
在复平面内以原点为圆心，半径为1的圆称为单位圆。以((1,0))为起点，将单位圆(n)等分，这些等分点是(n)个单位根。记第一个单位根为(omega_n)，根据复数相乘模长相乘，幅角相加的原则，第(k)个单位根记作(omega_n^k)。
根据三角函数推得单位根的表示方法：(omega_n^k=cosfrac{2kpi}{n}+isinfrac{2kpi}{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_n^k)=A_1(omega_{frac{n}{2}}^{k})+omega_n^kA_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}}^{k})-omega_n^kA_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)=y_0+y_1x+...+y_{n-1}x^{n-1} ]

此时多项式(B(x))是系数表示的。我们依次将单位根的倒数（(omega_n^0,omega_n^{-1},...,omega_n^{1-n})）代入得到一个点值表达式((z_0,z_1,...,z_{n-1}))。

[egin{aligned} z_k &= sum_{i=0}^{n-1}y_i(omega_n^{-k})^i\ &= sum_{i=0}^{n-1}(sum_{j=0}^{n-1}a_j*omega_n^{ij})(omega_n^{-k})^i\ &= sum_{i=0}^{n-1}sum_{j=0}^{n-1}a_j*(omega_n^{j-k})^i\ &= sum_{i=0}^{n-1}a_j(sum_{j=0}^{n-1}(omega_n^{j-k})^i)\ end{aligned} ]

而因为后半部分(sum_{j=0}^{n-1}(omega_n^{j-k})^i)可由等比数列求和公式得：

[sum_{j=0}^{n-1}(omega_n^{j-k})^i=dfrac{(omega_n^{j-k})^n-1}{omega_n^{j-k}-1} ]

注意，要使此式有意义，需满足分母不为0。而当(j=k)时分母为0。因此分类讨论：

当(j eq 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)

综上

[z_k = sumlimits_{i=0}^{n-1}a_j(sumlimits_{j=0}^{n-1}(omega_n^{j-k})^i) = na_k ]

所以

[a_k=dfrac{z_k}{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 ( ext{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("
");
	return 0;
}