Qinshihuang's Blog

tmc的小窝

0%

多项式与生成函数

发表于

咍咍。

前置知识

  • 复数。
  • 微积分。
  • 泰勒展开。
  • 牛顿迭代。
  • 多项式加法。

乘法

定义 $n$ 阶单位根为 $\mathbf C$ 内 $x^n=1$ 的根，记作 $\omega_1,\omega_2,\dots,\omega_n$，其中 $\omega_n=1$。

不难发现如果按照辐角主值排序那么 $\omega_2=\omega_1^2,\omega_3=\omega_1^3\dots$，同时还可以得到计算 $\omega$ 的公式：

原因很简单，$1=\omega_1^n=\angle2\pi$，因此 $\arg(\omega_1)=\dfrac{2\pi}{n}$，进而推得其他单位根。

以下记 $n$ 阶本源单位根 $\omega_n=\cos\dfrac{2\pi}{n}+\mathrm i\sin\dfrac{2\pi}{n}$。

注意到多项式乘法的定义说明了由点值表示到点值表示的乘法是 $O(n)$ 的。因此我们考虑快速实现一般多项式表示与点值表示的转换。

我们考虑

那么

我们不妨令 $b_k=f(\omega_n^k)$，并且构造一个新的多项式 $g(x)=\sum\limits_{k=0}^{n-1}b_kx^k$，那么将单位根的共轭复数带入得到：

这里用到了单位根

的性质。利用向量等知识证明不难。

（那个 $\sum(\text{bla})$ 就是表示对那些东西求和，因为我感觉写成矩阵形式更加直观。）

类似的，我们可以得到 $g(\omega_n^{-2})=na_2,\cdots,g(\omega_n^{-k})=na_k$。这样，将单位根作为点值表示有天然的优势：它天生支持逆变换。

而 FFT 算法正是利用了单位根的这一性质，又增加上了分治优化而允许我们在 $O(n2^n)$ 的时间复杂度内完成度数为 $2^n$ 的多项式的点值/一般多项式表示相互转化。

具体的，设多项式 $A(x)=\sum\limits_{k=0}^{n-1}a_kx^k$。我们将 $A(x)$ 分成两部分：

并记这两部分分别为 $A_1(x^2),A_2(x^2)$，也就是说 $A(x)=A_1(x^2)+xA_2(x^2)$。

现在让我们看看如果把 $\omega_n^i$ 带入会发生什么。

好了，$A$ 被分成了两个结构相同的子问题，这就可以使用分治/迭代实现，从而做到 $O(n\log n)$ 的时间复杂度。

核心代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
il void fft(cp *a,int inv) {
	for(int i=0;i<lim;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int mid=1;mid<lim;mid<<=1) {
		cp wn=cp(cos(pi/mid),inv*sin(pi/mid));
		for(int R=mid<<1,j=0;j<lim;j+=R) {	//迭代减常数
			cp w=cp(1,0);
			for(int k=0;k<mid;k++,w=w*wn) {
				cp x=a[j+k],y=w*a[j+mid+k];	//蝴蝶变换减常数。其实就是减少计算次数啦（
				a[j+k]=x+y;
				a[j+mid+k]=x-y;
			}
		}
	}
}

NTT

一般的多项式题里都会有取模。而 FFT 在取模下一无是处，因此要使用在模意义下的算法 NTT。

其实 NTT 就是把单位根 $\omega_n$ 换成了原根的幂 $g_n=g^{\frac{p-1}{n}}$。

我们需要使得这个新的替代品在模意义下与单位根的性质相同。

首先

这是原根的性质之一。

其次

而 $\left(g^{\frac{p-1}{2}}\right)^2\equiv g^{p-1}\equiv1\pmod p$，但 $g,g^2,\dots,g^{p-1}$ 模 $p$ 两两不同余，因此 $g_n^{\frac{n}{2}}\equiv-1\pmod p$。

最后

因此我们发现 $g_n$ 在模意义下与 $\omega_n$ 是等价的，可以直接替代。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
il void ntt(ll *A,ll inv) {
	for(int i=0;i<lim;i++) if(i<r[i]) swap(A[i],A[r[i]]);
	for(int mid=1;mid<lim;mid<<=1) {
		ll gn=qpow(inv==1?g:invg,(mod-1)/(mid<<1));
		for(int R=mid<<1,j=0;j<lim;j+=R) {
			ll G=1;
			for(int k=0;k<mid;k++,G=(G*gn)%mod) {
				ll x=A[j+k]%mod,y=G*A[j+mid+k]%mod;
				A[j+k]=(x+y)%mod;
				A[j+mid+k]=(x-y+mod)%mod;
			}
		}
	}
}

附：常用分治模数及其最小原根。

模数 分解 $g$
$65537$ $2^{16}+1$ $3$
$469762049$ $2^{26}\times7$ $3$
$998244353$ $2^{23}\times 7\times 17+1$ $3$
$1004535809$ $2^{21}\times479+1$ $3$

多项式求逆

定义：多项式 $f(x)$ 的逆为另一个多项式 $g(x)$ 满足

求满足条件的多项式 $f(x)g(x)\equiv 1\pmod{x^n}$，系数对 $998244353$ 取模。

我们考虑递归求解。

不妨设我们知道了一个多项式 $G(x)$ 使得 $f(x)G(x)\equiv1\pmod{x^{\frac{n}{2}}}$。

那么考虑做差

两边同时乘上 $f(x)$，根据定义有

化简得

可用 NTT 求解。递归单层时间复杂度是 $O(n\log 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
il void solve(ll len,ll *a,ll *b) {
	if(len==1) {
		b[0]=qpow(a[0],mod-2);
		return;
	}
	solve((len+1)>>1,a,b);
	lim=1,cnt=0;
	while(lim<(len<<1)) {
		lim<<=1;
		cnt++;
	}
	for(int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
	for(int i=0;i<len;i++) c[i]=a[i];
	for(int i=len;i<=lim;i++) c[i]=0;
	NTT(c,1);
	NTT(b,1);
	for(int i=0;i<lim;i++) b[i]=1ll*(2ll-b[i]*c[i]%mod+mod)%mod*b[i]%mod;
	NTT(b,-1);
	int invl=qpow(lim,mod-2);
	for(int i=0;i<len;i++) b[i]=1ll*b[i]*invl%mod;
	for(int i=len;i<lim;i++) b[i]=0;
}

多项式对数函数

在模 $x^n$ 意义下求一个多项式 $g(x)$ 使得

系数对 $9998244353$ 取模。

对这个式子做些变换：

直接求导、求逆、乘起来再积回去就可以了，时间复杂度仍然是 $O(n\log n)$。

核心代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
getInv(n,a,inva);
ll lim=1,cnt=0;
while(lim<(n<<1)) {
	lim<<=1;
	cnt++;
}
for(int i=0;i<lim;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
NTT(da,lim,1);
NTT(inva,lim,1);
for(int i=0;i<lim;i++) ans[i]=da[i]*inva[i]%mod;
NTT(ans,lim,-1);
ll invl=qpow(lim,mod-2);
for(int i=0;i<n;i++) ans[i]=ans[i]*invl%mod;
for(int i=n-1;i>0;i--) ans[i]=ans[i-1]*qpow(i,mod-2)%mod;
ans[0]=0;

多项式指数

给定多项式 $f(x)$，在模 $x^n$ 意义下求 $g(x)$ 使得

牛顿迭代

前置一个知识：什么是牛顿迭代。

假设你要求 $\sqrt n$，一个简单的方法是二分，这样需要计算 $\log n$ 次。

另一个思路就是使用牛顿迭代。这个问题的本质是求解 $f(x)=x^2-n$ 的零点。我们不妨假设一个初始值 $x_0$，并求出 $f(x)$ 在 $x=x_0$ 处的切线

与 $x$ 轴的坐标 $x_1$，最后将初始值设为 $x_1$ 继续迭代。这样做每次精度会翻倍。

应用在多项式也是如此。不妨设我们要求一个函数 $G(x)$ 使得对于给定的函数 $F(x)$ 有 $F(G(x))\equiv 0\pmod {x^n}$。我们假设已经知道了 $G_0(x)$ 使得 $F(G_0(x))\equiv 0\pmod{x^{\frac{n}{2}}}$，此时我们在 $G_0(x)$ 处将 $F(G(x))$ 泰勒展开：

由于 $G_0(x)\equiv G(x)\pmod{x^{\frac{n}{2}}}$，所以 $G_0^2(x)\equiv G^2(x)\pmod {x^n}$。

而根据二项式展开定理，

于是得到泰勒展开式第三项往后都为 $0$，因此得到

最后得到了牛顿迭代公式

有了这个公式之后我们就可以解决我们的问题啦！

本题中 $g(x)\equiv \mathrm e^{f(x)}$ 等价于 $f(x)\equiv\ln g(x)$，即 $\ln g(x)-f(x)\equiv 0$。

因此构造函数 $A(g(x))=\ln g(x)-f(x)$，套用牛顿迭代公式可知

然后就会了，一次迭代时间复杂度 $O(n\log n)$，总时间复杂度还是 $O(n\log n)$。

边界条件是 $n=1$，此时由于 $f_0=0$，所以 $\ln g(x)=0$，$g(0)=1$，直接返回即可。

核心代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
il void getExp(ll len,ll *a,ll *expa,ll *lna,ll *inva) {
	if(len==1) {
		expa[0]=1;
		return;
	}
	getExp((len+1)>>1,a,expa,lna,inva);
	getLn(len,expa,lna,inva);
	ll lim=1,cnt=0;
	while(lim<(len<<1)) {
		lim<<=1;
		cnt++;
	}
	for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(cnt-1));
	for(int i=0;i<len;i++) a1[i]=a[i];
	for(int i=len;i<lim;i++) lna[i]=a1[i]=0;
	for(int i=0;i<lim;i++) a1[i]=((a1[i]-lna[i])%mod+mod)%mod;
	a1[0]=(a1[0]+1)%mod;
	NTT(a1,lim,1);NTT(expa,lim,1);
	for(int i=0;i<lim;i++) expa[i]=expa[i]*a1[i]%mod;
	NTT(expa,lim,-1);
	for(int i=len;i<lim;i++) expa[i]=0;
}

应用：多项式快速幂

给定一个多项式 $f(x)$，求模 $x^n$ 意义下 $f^k(x)$ 的值。

暴力做法是做 $k$ 次 NTT 乘法，时间复杂度是 $O(nk\log n)$。

次暴力做法是快速幂上套 NTT，时间复杂度是 $O(n\log n\log k)$，但是 $k\le10^{10^5}$ 的量级让这玩意也过不去。

我们想一下能不能加速一下这个过程。

有一个公式：

带入到多项式中也是一样的，即

注意到这题 $f_0=1$，因此直接 $\ln,\exp$ 即可。

对于 $k$ 要在读入的时候取模的事情，因为

直接把 $k\ln f(x)$ 带入得到

这东西显然可以取模。

然而加强版把这个限制去掉了，怎么办呢？

很简单，让整个多项式的系数除以 $f_0$，这样 $f_0$ 就是 $1$ 了。

你说得对，但是数据里有 $f_0=0$ 的点。那就整体除以 $f_1x$。你说得对，但是数据里有 $f_0=f_1=0$ 的点。那就整体除以 $f_2x^2$……

直到得到一个常数项为 $1$ 的多项式之后直接做就可以了。时间复杂度是 $O(n\log n)$。

生成函数

前置知识：

  • 多项式。
  • 组合数学。
  • 推式子。

什么是生成函数

考虑如下的一个数列

当然可以写成这样的形式，但是另一种方式是将其写成生成函数的形式（普通生成函数），即

这里 $z$ 起到了一个占位符的作用，与一般意义上函数自变量不同。当然也可以当函数自变量理解啦

对于一个生成函数 $A(z)$，取其 $z^n$ 项上的系数 $a_n$，记作 $[z^n]A(z)$。

生成函数相加是可行的，即

但是这种形式幂级数不方便计算，一般来说我们要将其写成封闭形式。

例如，对于序列 $a=\left<1,1,1,\dots\right>$ 的生成函数 $A(z)=\sum\limits_{k=0}^\infty z^i$，我们可以发现

因此 $(1-z)A(z)=1$，得到 $A(z)$ 的封闭形式

类似的，我们可以得到如下公式：

事实上 $A(z)$ 的意义不止于此。考虑另一个生成函数 $B(z)=\sum\limits_{k=0}^\infty b_kz^k$，如果将这两个函数乘起来，可以得到一个新的函数 $C(z)$。看看这个 $C(z)$ 是什么：

所以说 $C(z)$ 的系数是 $B(z)$ 的系数前缀和。这是相当有用的。

另外我们也可以发现 $A(z)$ 的乘法逆 $A^{-1}(z)=1-z$ 就是差分运算。

根据这个例子，我们发现生成函数的卷积就是将 $A$ 中的 $z^0,z^1,\dots,z^n$ 项与 $B$ 中的 $z^n,z^{n-1},\dots,z^0$ 项相搭配得到新的系数。

生成函数能干什么

推数列通项

不妨拿斐波那契数列举例子。

众所周知，$f_1=f_2=1,f_n=f_{n-1}+f_{n-2}\ (n\ge3)$ 这个数列就是斐波那契数列。现在我们想要它的通项公式。

如果你背过了，就知道

但这是怎么推出来的呢？

不妨把斐波那契数列的生成函数写出来：

我们希望先找到它的封闭形式。

然后考虑取 $[z^n]F(z)$。我们考虑使用以前知道的公式进行求解。

我们不妨设 $F(z)=\dfrac{c}{1-az}+\dfrac{d}{1-bz}$。

开始推式子：

对照系数可知：

解得

这里不把 $c,d$ 写出来，一部分是因为难算，还有一部分原因是后面会消去。

将结果代入得到：

因此

最后把 $a,b$ 代入就得到了大名鼎鼎的公式：

有意思的是，一道计算题也可以用 $F(z)$ 做。

求证：

的结果是有理数，并求出这个结果。

其实题目就是让求

代入公式得到 $F\left(\dfrac{1}{10}\right)=\dfrac{10}{89}$。答案确实是这个：

$\dfrac{10}{89}\approx0.11235955056179775280898876404494\dots$

数数

比如说你有 $3$ 种水果。苹果最多吃 $3$ 个，香蕉最多吃 $1$ 个，橘子必须吃 $4$ 的倍数个，问吃 $8943477323$ 个水果有多少种方法。怎么做？

不是老弟，吃 $n$ 个水果的话 $n$ 要和苹果香蕉总个数模 $4$ 同余，除了吃 $0$ 个以外都只有 $2$ 种情况啊？

额。

我们写出三种水果的数量生成函数：

直接卷起来，得到答案生成函数：

展开得到：

bruh.

什么考场上生成函数太大了没法得到封闭形式？

这时候多项式算法就派上用场了。你直接写出形式幂级数，对 $z^n$ 取模就可以得到 $z^0,z^1,\dots,z^{n-1}$ 的系数。而这东西拿多项式写是 $O(n\log n)$ 的。也就是说你可以直接操纵生成函数了！

例题：P4389

完全背包，但是 $n,m,v_i\le10^5$，对每个 $1\le s\le m$ 求装到 $s$ 体积的方案数。模 $998244353$。

首先跑个 dp $O(nm)$。

注意到对于一个体积为 $v$ 的物品它的生成函数是

因此我们可以简单地将 $1-z^{v_i}$ 卷积起来再求个逆，时间复杂度是……$O(nm\log m)$。: (

直接卷积不可行，考虑先求 $\ln$，系数加起来，再 $\exp$ 回去。这样的时间复杂度是……$O(nm\log m)$。:<

我们看看能不能直接求出 $\ln(1-z^c)$。

发现这个东西模 $x^{m+1}$ 意义下只有 $O\left(\dfrac{m}{c}\right)$ 个点有值，就做完了，直接预处理出 $\ln ans$，时间复杂度 $O(m\ln m)$。最后 $\exp$ 回去即可得到答案。