「多项式」学习笔记

多项式乘法

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

\[c_n = \sum\limits_{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=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}}{k})+\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}}{k})-\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)=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})\)

\[\begin{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 \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\)

综上

\[z_k = \sum\limits_{i=0}^{n-1}a_j(\sum\limits_{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 (\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;
}
posted @ 2019-07-11 14:23  DennyQi  阅读(519)  评论(0编辑  收藏  举报