多项式与生成函数学习笔记

发表于 更新于 本文字数： 1.5k 阅读时长 ≈ 6 分钟

基本概念

定义

记一个关于 xx 的式子：

f(x)=∑i=0nai×xif(x)=\sum\limits_{i=0}^na_i\times x^i

为一个 nn 次多项式。其中 aia_i 为常数，nn 为该多项式的度，也称次数，记作 deg⁡f\deg f。

显然，f(x)f(x) 可以看作一个关于 xx 的 nn 次函数 y=f(x)y=f(x)。

根据我们的初中常识，n+1n+1 个点的坐标可以唯一确定一个 nn 次函数的解析式，进而可以确定这个多项式，这种方式被称为多项式的点值表示法.

对应的，写成 f(x)=∑i=0nai×xif(x)=\sum\limits_{i=0}^na_i\times x^i 形式被称为多项式的系数表示法.

多项式加减法

设 nn 次多项式 A(x)=∑i=0nai×xiA(x)=\sum_{i=0}^n a_i\times x^i，B(x)=∑i=0nbi×xiB(x)=\sum_{i=0}^nb_i\times x^i，则：

f(x)=A(x)±B(x)=∑i=0n(ai±bi)xif(x)=A(x)\pm B(x)=\sum\limits_{i=0}^n(a_i\pm b_i)x^i

多项式乘法

f(x)=A(x)×B(x)=∑i=0n∑j=0nai×bj×xi+jf(x)=A(x)\times B(x)= \sum\limits_{i=0}^n\sum\limits_{j=0}^n a^i\times b^j\times x^{i+j}

使用点值表示法时，直接将其横坐标对应的纵坐标进行加减乘即可，不过做乘法时，需要 A(x)A(x) 和 B(x)B(x) 各给出 2n2n 组点值。

快速傅里叶变换（FFT）

前置知识：复数、单位根

FFT 能够支持在 O(nlog⁡n)O(n\log n) 的时间内计算两个度为 nn 的多项式的乘法。

我们知道系数表示法需要 O(n2)O(n^2) 的时间复杂度，而点值表示法只需要 O(n)O(n) ，因为只需要将对应点值的纵坐标相乘就可以了。因此我们考虑将一个多项式从系数表示法改为点值表示法，相乘后再改回系数表示法。如果这两个操作的时间复杂度能够低于 O(n2)O(n^2) ，那么就可以得到一个比暴力更优的做法。

而求出 nn 次多项式在每个 nn 次单位根下的点值的过程叫做离散傅里叶变换（DFT），将这些点值重新插值为系数表示法的过程叫做离散傅里叶逆变换（IDFT）。

以下设进行变换的多项式为 (n−1)(n - 1) 次多项式 A(x)=∑i=0n−1ai×xiA(x) = \sum_{i = 0}^{n - 1} a_i \times x^i。其中 nn 是 22 的整数幂。如果不是的话，则向 A(x)A(x) 的更高次数位 nn 补充 an=0a_n = 0 ，令其成为 nn 次多项式，一直进行直到其次数加一的值是 22 的整数幂，取 nn 等于其次数，m=n2m = \frac{n}{2}。

DFT

考虑求出一个长度为 nn 的数列 {bi}\{b_i\}，这个数列的第 kk 项为 A(x)A(x) 在 nn 次单位根的 kk 次幂处的点值。

因此有

bk = ∑i=0n−1ai×ωnib_k~=~\sum_{i = 0}^{n - 1} a_i \times \omega_{n}^i

这个过程是 O(n2)O(n^2) 的，可以使用 FFT 来优化。

FFT

我们对 A(x)A(x) 按照系数的编号的奇偶性进行分类：

A(x) = ∑i=0n−1aixi = ∑i=0m−1a2i×x2i+∑i=0m−1a2i+1×x2i+1= ∑i=0m−1a2ixi2+x∑i=0m−1a2i+1x2i = ∑i=0m−1a2i(x2)i+x∑i=0m−1a2i+1(x2)i\begin{aligned} A(x)~&=~\sum_{i = 0}^{n - 1} a_i x^i~\\ &=~\sum_{i = 0}^{m-1} a_{2i} \times x^{2i} + \sum_{i = 0}^{m-1} a_{2i + 1} \times x^{2^i + 1}\\ &=~\sum_{i = 0}^{m - 1} a_{2i} x^{i2} + x\sum_{i = 0}^{m - 1} a_{2i + 1} x^{2i}~\\ &=~\sum_{i = 0}^{m - 1} a_{2i} (x^{2})^{i} + x\sum_{i = 0}^{m - 1} a_{2i + 1} (x^2)^{i} \end{aligned}

设 A0(x)A_0(x) 是一个 (m−1)(m - 1) 次多项式，满足

A0 = ∑i=0m−1a2ixiA_0~=~\sum_{i = 0}^{m - 1} a_{2i} x^i

设 A1(x)A_1(x) 是一个 (m−1)(m - 1) 次多项式，满足

A1 = ∑i=0m−1a2i+1xiA_1~=~\sum_{i = 0}^{m - 1} a_{2i + 1}x^i

联立以上三式，可以得到

A(x) = A0(x2)+x×A1(x2)A(x)~=~A_0(x^2) + x \times A_1(x^2)

接着，参考单位根的性质：

(ωnk)2=ωmkωnk+m=−ωnk(\omega_n^k)^2=\omega_m^k\\ \omega_n^{k+m}=-\omega_n^k

我们可以分开考虑小于 mm 的点值和大于 mm 的点值。

对于 0≤k<m0\leq k<m，有：

A(ωnk) = A0((ωnk)2)+wnkA1((ωnk)2)= A0(ωmk)+wnkA1(ωmk)\begin{aligned} A(\omega_n^k)~&=~A_0((\omega_n^k)^2) + w_n^kA_1((\omega_n^k)^2)\\ &=~A_0(\omega_m^k) + w_n^kA_1(\omega_m^k) \end{aligned}

再考虑大于 mm 次的点值：

A(ωnm+k) = A0((ωnm+k)2)+ωnm+kA1((ωnm+k)2)= A0((wnk)2)+−ωnkA1((ωnk)2)= A0(ωmk)−wnkA1(ωmk)\begin{aligned} A(\omega_n^{m + k})~&=~A_0((\omega_n^{m + k})^2) + \omega_n^{m + k} A_1((\omega_n^{m + k})^2)\\ &=~A_0((w_n^k)^2)+-\omega_n^k A_1((\omega_n^k)^2)\\ &=~A_0(\omega_m^k) - w_n^kA_1(\omega_m^k) \end{aligned}

可以发现，只要求出 A0A_0 和 A1A_1 小于 mm 次的点值，就可以线性求出 AA 在整个 nn 次幂处的点值，而 A0A_0 和 A1A_1 在小于 mm 次的点值是可以递归求解的。时间复杂度 O(nlog⁡n)O(n\log n)。

IDFT

至此我们得到了 bk = ∑i=0n−1ai×ωnib_k~=~\sum_{i = 0}^{n - 1} a_i \times \omega_{n}^i，而现在希望还原系数数列 {ai}\{a_i\}。

这里直接给出结论，推导和证明请参考其它文章：

ak = 1n∑i=0n−1biωn−kia_k~=~\frac{1}{n} \sum_{i = 0}^{n - 1} b_i \omega_n^{-ki}

仍然要尝试优化，那就是 IFFT，即快速傅里叶逆变换。

IFFT

由于 wn−kw_n^{-k} 可以看做 nn 次本原单位根每次逆时针旋转本原单位根幅角的弧度，因此 ωn−k\omega_n^{-k} 和 ωnk\omega_n^k 是一一对应的。具体的，wn−k=wnk+nw_n^{-k} = w_n^{k + n}。因此我们只需要使用 FFT 的方法，求出 B(x)B(x) 在 ωn\omega_n 各个幂次下的值，然后数组反过来，即令 ak = 1n∑i=0nB(wnn−k)a_k~=~\frac{1}{n} \sum_{i = 0}^n B(w_n^{n - k}) 即可。

时间复杂度 O(nlog⁡n)O(n\log 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
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=6e6+5;
const double pi=acos(-1);
int n,m,sm,k;
complex<double> f[maxn],g[maxn];
int tax[maxn];
inline void makerev(int x)
{
	int d=x>>1,num=0;
	tax[num++]=0,tax[num++]=d;
	for(int w=2;w<=x;w<<=1)
	{
		d>>=1;
		for(int p=0;p<w;p++) tax[num++]=tax[p]|d;
	}
}
inline void fft(complex<double> *a,int x)
{
	for(int i=1;i<x;i++) 
        if(tax[i]>i) swap(a[i],a[tax[i]]);
	for(int len=2,mid=1;len<=x;mid=len,len<<=1)
	{
		complex<double> add(cos(pi/mid),sin(pi/mid)),init(1.0,0.0);
		for(int l=0,r=len-1;r<=x;l+=len,r+=len)
		{
			auto w=init;
			for(int p=l;p<l+mid;p++)
			{
				auto x=a[p]+w*a[p+mid];
                auto y=a[p]-w*a[p+mid];
				a[p]=x,a[p+mid]=y;
				w*=add;
			}
		}
	}
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(0),cout.tie(0);
	cin>>n>>m;
	for(int i=0;i<=n;i++) 
	{
		int tmp;
		cin>>tmp;
		f[i]=tmp;
	}
	for(int i=0;i<=m;i++)
	{
		int tmp;
		cin>>tmp;
		g[i]=tmp;
	}
	sm=n+m,k=1;
	while(k<=sm) k<<=1;
	makerev(k);
	fft(f,k);
	fft(g,k);
	for(int i=0;i<k;i++) f[i]*=g[i];
	fft(f,k);
	reverse(f+1,f+k);
	for(int i=0;i<=sm;i++)  
        cout<<int(f[i].real()/k+0.5)<<" ";
	return 0;
}

关于如何递归转迭代，先挖个坑，以后再补。

未完待续……