【算法专题】多项式运算与生成函数

【快速傅里叶变换】FFT

参考：从多项式乘法到快速傅里叶变换 by miskcoo

FFT 学习笔记 by Menci

（一）多项式的表示法

系数表示法：f(x)=a[n-1]*x^(n-1)+...+a[0]，称为n-1次多项式。

点值表示法：一个n-1次多项式在复数域中有n个根，即n个(x,y)可以唯一确定一个n-1次多项式。

对于一个多项式，从其系数表示法到其点值表示法的变换称为离散傅里叶变换（DFT），反之称为傅里叶逆变换（IDFT）。

朴素的离散傅里叶变换，枚举实现的复杂度为O(n^2)。

快速傅里叶变换是指以O(n log n)的复杂度实现IDF和IDFT的算法，常用Cooley-Tukey算法。

（二）复数

复数是形如a+bi的数，当b=0时为实数。

定义一个平面为复平面，那么平面内的每个点(a,b)唯一对应一个复数a+bi，i可以理解为y轴上的单位长度，正如1是x轴上的单位长度。

i的本质是在数轴上定义旋转变换，i是1逆时针旋转90°，那么i^2=-1。

复数相加，遵循平行四边形定则。

复数相乘，模长相乘，幅角相加。

（三）单位根

以圆点为起点，以复平面单位圆的n等分点为终点，作n个向量，设所得幅角为正且最小的向量对应的复数为ω(1,n)，即n次单位根。（括号左为上标，右为下标）。

#include<cstdio>
#include<cstring>
#include<cctype>
#include<cmath>
#include<queue>
#include<stack>
#include<set>
#include<vector>
#include<complex>
#include<algorithm>
#define ll long long
#define lowbit(x) x&-x
using namespace std;
int read(){
char c;int s=0,t=1;
while(!isdigit(c=getchar()))if(c=='-')t=-1;
do{s=s*10+c-'0';}while(isdigit(c=getchar()));
return s*t;
}
int min(int a,int b){return a<b?a:b;}
int max(int a,int b){return a<b?b:a;}
int ab(int x){return x>0?x:-x;}
//int MO(int x){return x>=MOD?x-MOD:x;}
//void insert(int u,int v){tot++;e[tot].v=v;e[tot].from=first[u];first[u]=tot;}
/*------------------------------------------------------------*/
const int inf=0x3f3f3f3f,MOD=1e9+7,maxn=300010;//
const double PI=acos(-1);
namespace fft{
complex<double>o[maxn],oi[maxn];
void init(int n){
for(int k=0;k<n;k++){o[k]=complex<double>(cos(2*PI*k/n),sin(2*PI*k/n));oi[k]=conj(o[k]);}
}
void transform(complex<double>*a,int n,complex<double>*o){
int k=0;
while((1<<k)<n)k++;
for(int i=0;i<n;i++){
int t=0;
for(int j=0;j<k;j++)if(i&(1<<j))t|=(1<<(k-j-1));
if(i<t)swap(a[i],a[t]);
}
for(int l=2;l<=n;l*=2){
int m=l/2;
for(complex<double>*p=a;p!=a+n;p+=l){
for(int i=0;i<m;i++){
complex<double>t=o[n/l*i]*p[i+m];
p[i+m]=p[i]-t;
p[i]+=t;
}
}
}
}
void dft(complex<double>*a,int n){transform(a,n,o);}
void idft(complex<double>*a,int n){
transform(a,n,oi);
for(int i=0;i<n;i++)a[i]/=n;
}
}
int n,N,m,kind,f[maxn][30],h[30],F[maxn],g[2][maxn];
complex<double>a[maxn],b[maxn],c[maxn],d[maxn],v[maxn];
ll ans1[maxn],ans2[maxn],ans3[maxn],ans4[maxn];//
void multply(complex<double>*x,complex<double>*y,ll *z){
for(int i=0;i<N;i++)v[i]=x[i]*y[i];
fft::idft(v,N);
for(int i=0;i<n;i++)z[i]=(ll)(v[i].real()+0.5);
}
void MTT(int *x,int *y,int *z){
for(int i=0;i<N;i++)a[i]=b[i]=c[i]=d[i]=complex<double>(0,0);//
for(int i=0;i<n;i++)a[i].real(x[i]>>15);
for(int i=0;i<n;i++)b[i].real(x[i]&32767);
for(int i=0;i<n;i++)c[i].real(y[i]>>15);
for(int i=0;i<n;i++)d[i].real(y[i]&32767);
fft::dft(a,N);fft::dft(b,N);fft::dft(c,N);fft::dft(d,N);
multply(a,c,ans1);multply(a,d,ans2);
multply(b,c,ans3);
multply(b,d,ans4);
for(int i=0;i<n;i++)z[i]=(ans1[i]*32768%MOD*32768%MOD+(ans2[i]+ans3[i])*32768%MOD+ans4[i])%MOD;
}
int main(){
n=read();m=read();kind=read();
f[0][0]=1;h[0]=1;
int mx=0;
for(mx=1;h[mx-1]<n;mx++)h[mx]=h[mx-1]*m;mx--;
for(int i=1;i<=n;i++){
for(int j=0;j<=mx;j++)if(i-h[j]>=0){
for(int k=0;k<=j;k++)f[i][j]=(f[i][j]+f[i-h[j]][k])%MOD;
}
}
for(int i=0;i<=n;i++)for(int j=0;j<=mx;j++)F[i]=(F[i]+f[i][j])%MOD;
int x=0;
for(int i=0;i<=n;i++)g[x][i]=F[i];
n++;
N=1;//
while(N<n+n)N*=2;fft::init(N);
for(int k=1;k<kind;k++)MTT(g[x],F,g[1-x]),x=1-x;
int sum=0;
for(int i=0;i<=n;i++)sum=(sum+g[x][i])%MOD;
printf("%d",sum);
return 0;
}

View Code 注意清空complex的时候实部和虚部(imag)要一起清空。

好看的模板：L_0_Forever_LF

 

【快速数论变换】NTT

只记录重要概念，证明略去。

（一）原根

当(a,m)=1时，对于满足a^x=1(%m)的最小正整数x，称x为a模m的阶。

根据欧拉定理a^φ(m)=1(%m)，当x=φ(m)时，称a为m的原根。

以下只讨论m为素数的情况，则当a为m的原根时，a^0~a^(p-2)取遍1~p-1所有值。

模m有原根的充要条件：m=2,4,p^e,2*p^3，p是奇素数。（也就是说，m为素数时一定有原根）

求m的原根：p-1= p1^a1 * p2^a2 * pk^ak，g是p的原根当且仅当对于所有的pi满足g^[ (p-1)/pi ] ≠ 1 (%p)

例题：【51NOD】1135 原根

（二）快速数论变换

当模数为形如p=r*2^k+1的素数（费马素数）时，则有n|p-1，可以进行NTT。

先找到p的原根g（p=998244353 || 1004535809，g=3）

在原来FFT的基础上，omega[i]=g^[ (p-1)/n*i ] % p，倒数为逆元。

IDFT时，除以n改为乘n的逆元。

（三）模数任意的NTT

找到三个费马素数满足相乘结果>n*(m-1)^2，分别进行NTT后用CRT合并。

p=998244353,1004535809,469762049，g=3。

我写的是FFT合并。

（四）离散对数

当(a,p)=1时，若满足a^x=b (%p)，则称在模p意义下，x是b以a为底的离散对数，即logab=x（单个快速求解可用BSGS算法）。

1.对于x*y=z(%p)，有log x+log y=log z(%p-1)，因此离散对数常用于乘法转加法（生成函数）。

2.对于x^y=z(%p)，有y*log x=log z(%p-1)。

其中log以p的原根g为底。

例题：【BZOJ】3992: [SDOI2015]序列统计 NTT+生成函数

【生成函数】母函数

生成函数的三大要素：①选择项，②大小，③元素个数。一般最终要求某个“大小”的元素个数。

对于一类组合对象构成的集合A：

1.每个元素a∈A都定义了一个非负整数的”大小“，记为|a|。

2.大小为n的元素个数记为$A_n$。

那么A的一般生成函数是

$$A(x)=\sum_{i=0}^{n}A_ix^i$$

在这里每个元素都抽象成”大小“，元素a可以理解为有|a|的单位元素的元素。

组合对象集合D为A和B的笛卡尔积，即D中的每个元素都是A中某个元素a和B中某个元素b组成的有序二元组(a,b)，那么显然有D(x)=A(x)B(x)。

若干一般生成函数的乘积中，第n项的含义是：每个选择项取一个元素，大小相加为n的元素个数。

每个生成函数本质上是一个集合，那么若干生成函数的乘积就是★每个集合取一个元素的组合。例如生成函数A,B,C，A*B*C的每个元素就是有序三元组(a,b,c)。

指数型生成函数是

$$A(x)=\sum_{i=0}^{n}A_i\frac{x^i}{i!}$$

这样就会有：

$$D_n=\sum_{i+j=n}A_iB_j\frac{(i+j)!}{i!j!}=\sum_{i+j=n}A_iB_j\binom{i+j}{i}$$

这里乘(i+j)!是因为这只是系数，后面要除以(i+j)!。

理解为每个元素内部有序就可以了，这样元素内部是排列。

生成函数都是处理对于n个选择项各选一个组成对应”大小“的元素个数，而一般生成函数元素内部是组合，指数型生成函数元素内部是排列。

一般生成函数还有个化简公式，令x∈[-1,1]时套等比数列公式即可收敛：

$$\sum_{i=0}^{n}x^i=\frac{1}{1-x}$$

指数型生成函数也有个化简公式——泰勒展开：

$$\sum_{i=0}^{n}\frac{x^i}{i!}=e^x$$

例题：

1.热身：苹果只能取偶数个，橘子只能取1~4个，求拿n个水果的方案数

题解：定义每种水果为一个集合（每个集合选一个），“大小”为水果个数，最后求总集合”大小“为n的元素个数。

f(x)=(1+x^2+x^4+...)*(1+x+x^2+x^3+x^4)。

如果觉得集合很难理解，不妨用”选择项“这个词。

2.【BZOJ】3771: Triple FTT+生成函数

题意：给定n个物品，价值为ai，物品价格互不相同，求选一个或两个或三个的价值为x的方案数，输出所有存在的x和对应方案数。ai<=40000。

题解：要求什么就定义什么为”大小“，所以定义”大小“为价值，[第一个物品][第二个物品][第三个物品]为三个选择项。

那么每个选择项的每个系数记录对应价值的物品数量（1个）。

这样拼起来就好了吗？不是，物品不能重复取，所以拼起来之后再容斥掉选相同的。

我们可以直接写出选一个物品的集合的生成函数f，两个相同物品的g和三个相同物品的h。

考虑有AAB,ABA,BAA,AAA四种不合法情况，答案就是f^3-3f*g+2h。最后这个求得排列数，需要/3!。选1个或2个的随便推推也一样。

3.【BZOJ】3992: [SDOI2015]序列统计 NTT+生成函数

题意：给定一个[0,m-1]范围内的数字集合S，从中选择n个数字（可重复）构成序列。给定x，求序列所有数字乘积%m后为x的序列方案数%1004535809。1<=n<=10^9，3<=m<=8000，m为素数，1<=x<=m-1。

题解：要求乘积，定义”大小“为数字的乘积。但是我们不能加减”大小“啊？

换成离散对数就可以了，然后定义每个数字为选择项，答案就是f^n。

4. [母函数]HDU 1521——排列组合 

题意：有n种物品，并且知道每种物品的数量。要求从中选出m件物品的排列数。例如有两种物品A,B，并且数量都是1，从中选2件物品，则排列有”AB”,”BA”两种。

题解：定义”大小“为物品数量，选择项为每种物品，那么组合数就是一般生成函数（元素内部有序，严格按物品编号排），排列数就是指数型生成函数（元素内部带标号，可以打乱）。

这里指数型生成函数，指的是方案就是元素内部带标号的方案数，这样再计算过程中那个公式自动/i!后进行计算再乘回(i+j)!的。

暴力枚举求解这个生成函数。

 

另一部分知识：生成函数公式的化简

参考：什么是生成函数？ by M67

首先根据二项式定理：

$$(1+x)^n=\sum_{i=0}^{n}\binom{n}{i}x^i$$

扩展到负数和实数，即广义二项式定理，只要将组合数表示成下降幂即可：

$$(1+x)^{n}=\sum_{i=0}^{\infty}\frac{n^{\underline{i}}}{i!}x^i$$

这里有一个特殊的变换，当n>0时：

$$(1-x)^{-n}=\sum_{i=0}^{\infty}\frac{(-n)^{\underline{i}}}{i!}(-x)^i$$

$$(1-x)^{-n}=\sum_{i=0}^{\infty}\frac{(-1)^i*n^{\overline{i}}}{i!}*(-1)^i*x^i$$

$$(1-x)^{-n}=\sum_{i=0}^{\infty}\binom{n+i-1}{n-1}*x^i$$

所以这是将i个相同的数分割n个非空部分的方案数的生成函数。

另外还常用等比数列递推公式来收敛等比数列：

$$\sum_{i=0}^{\infty}(x^q)^i=\frac{1}{1-x^q}$$

另外对于有限的等比数列还常用类似错位相减的方法，左边乘上(1-x)就会变成右边分子，即：

$$\sum_{i=0}^{n}x^i=\frac{1-x^{n+1}}{1-x}$$

例题：

1.求：

$$g(x)=(1+x^2+x^4+...)(1+x^5+x^{10}+...)(1+x+x^2+x^3+x^4)(1+x)$$

用上面提到的技巧即可得到：

$$g(x)=\frac{1}{1-x^2}*\frac{1}{1-x^5}*\frac{1-x^5}{1-x}*(1+x)$$

不断地约分，最后用平方差公式化简可得：

$$g(x)=\frac{1}{(1-x)^2}=(1-x)^{-2}=\sum_{i=0}^{\infty}(i+1)x^i$$ 

再代入上面那个结论就可以得到多项式了。

2.食物：用上面的技巧化简约分最后剩个小组合数，n=10^500用读入取模。

 

【多项式求逆】

核心原理是$\%x^{\frac{n}{2}}$的多项式平方后可以转化为$\%x^n$。

已知多项式f(x)，求多项式g(x)，满足f(x)g(x)=1(%x^n)。

$$f(x)g(x)=1(\%x^{\frac{n}{2}})$$

$$f(x)h(x)=1(\%x^n)$$

现在已知g(x)，求h(x)。

$$f(x)g(x)-1=0(\%x^{\frac{n}{2}})$$

将1移到左边后就可以平方了。

$$f(x)^2g(x)^2-2f(x)g(x)+1=0(\%x^n)$$

将1换成f(x)*g(x)，从而将h(x)代入，提出f(x)消去。

$$f(x)g(x)^2-2g(x)+h(x)=0(\%x^n)$$

最终得到：

$$h(x)=g(x)(2-f(x)g(x))(\%x^n)$$

然后就可以递归求解，边界条件:当n=1时，f(x)g(x)=1(%x)，g(0)是f(0)在%MOD意义下的数论逆元。

注意：每次n不同都要重新预处理Omega[]。                                                                                                                              

例题：【BZOJ】4555: [Tjoi2016&Heoi2016]求和 排列组合+多项式求逆 或 斯特林数+NTT

 

【拉格朗日插值法】

参考：拉格朗日插值法(图文详解) by Angel_Kitty

对于一个n次多项式，如果已知n+1个点，可以构造拉格朗日多项式L(x)：

$$L(x)=\sum_{i=0}^{n}y_il_i(x)$$

其中$l_i(x)$为插值基函数：

$$l_i(x)=\prod_{j=0,j\neq i}^{n}\frac{x-x_j}{x_i-x_j}$$

通过代入需要的x即可得到答案。

每次插值的复杂度为O(n^2)。

例题：【BZOJ】4559: [JLoi2016]成绩比较 计数DP+排列组合+拉格朗日插值

当横坐标连续时，上式可以表示为：

$$f(x)=\sum_{i=0}^{n}f(i)*\prod_{j=0,j\neq i}^{n}\frac{x-j}{i-j}$$

预处理阶乘的逆元$\frac{1}{i!}$，预处理$v=\prod x-j$，每次乘上v/(i-j)，分母是两段阶乘，再根据负数的个数判断正负性。

复杂度O(n)。

例题：【BZOJ】3453: tyvj 1858 XLkxc 拉格朗日插值（自然数幂和）