【算法专题】多项式运算与生成函数
【快速傅里叶变换】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次单位根。（括号左为上标，右为下标）。
图片来源：其中B点是单位根ω(1,n)，逆时针依次为ω(2,n)，ω(3,n)...，ω(n,n)=ω(n,0)=1。
计算公式：ω(k,n)=cos ( 2kπ/n ) + i*sin ( 2kπ/n )
单位根的性质：
(1)消去：ω(2n,2k)=ω(n,k)
(2)折半：ω(n,k+n/2)=-ω(n,k)
将ω(n,0)~ω(n,n-1)这n个单位根作为代表n-1次多项式的n个点的横坐标，可以得到很好的性质。
（四）快速傅里叶变换（FFT解决DFA）
这部分因为不会操作数学公式，直接粘贴Menci博客QAQ。
将n-1次多项式A(x)的系数奇偶分成两个多项式A1(x)和A2(x)，则A(x)=A1(x^2)+x*A2(x^2)。
对于k<n/2，有A(ω(n,k))=A1(ω(n/2,k)) + ω(n,k)*A2(ω(n/2,k))
同时，有A(ω(n,k+n/2))=A1(ω(n/2,k)) - ω(n,k)*A2(ω(n/2,k))
对于一个k次多项式，通过奇偶分项得到两个k/2次多项式，分别计算后再调用其值解决k次多项式，即分治解决。
（五）傅里叶逆变换（IDFA）
对于n-1次多项式，其n-1维系数向量{a0,a1...an-1}通过DFA得到点值向量{b0,b1...bn-1}，反之操作称为IDFA。
将点值向量作为系数，以单位根的倒数进行FFT，得到的每个数除以n，就是IDFA的结果。
（六）迭代实现FFT
对于多项式A(x)，已知系数向量a[]，求横坐标为ω(n,0)~ω(n,n-1)的点值向量b[]。
将多项式奇偶分项后，对于k<n/2，有A(ω(n,k))=A1(ω(n/2,k)) + ω(n,k)*A2(ω(n/2,k))，同时有A(ω(n,k+n/2))=A1(ω(n/2,k)) - ω(n,k)*A2(ω(n/2,k))，分治边界是a[i]*ω(0,0)即a[i]。
边界元素：FFT递归边界的数组排布恰好是原数组每个位置二进制反转后的数字，例如：
原：00 01 10 11
终：00 10 01 11
蝴蝶操作：为了在合并时不引入新数组，进行一下操作。ω(l,k)=ω(n,n/l*k)，预处理以n为底的ω[]，IDFT时预处理倒数。
t=ω(n/l*k)*a[l/2+k]
a[l/2+k]=a[k]-t
a[k]+=t
（七）多项式乘法
多项式的点值表示法易于进行乘法，因为对于fc(x)=fa(x)*fb(x)，每个点x在多项式A,B中的点值相乘即可得到在多项式C中的点值。
将n-1次多项式A和m-1次多项式B均视为n+m-2次多项式（高位补0），进行DFT后相乘再通过IDFT即可得到多项式C。
#include<cstdio> #include<cstring> #include<cctype> #include<cmath> #include<queue> #include<stack> #include<set> #include<vector> #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; const int maxn=300010;//2^18!!! const double PI=acos(-1); int n,m; struct cp{ double x,y; cp(double a,double b){x=a;y=b;} cp(){x=y=0;}; cp operator + (cp a){return cp(x+a.x,y+a.y);} cp operator - (cp a){return cp(x-a.x,y-a.y);} cp operator * (cp a){return cp(x*a.x-y*a.y,x*a.y+y*a.x);} }a[maxn],b[maxn]; void fft(cp *a,int n,int f){ int k=0; for(int i=0;i<n;i++){ if(i>k)swap(a[i],a[k]); for(int j=n>>1;(k^=j)<j;j>>=1); } for(int l=2;l<=n;l<<=1){ int m=l/2; cp wn(cos(2*PI*f/l),sin(2*PI*f/l)); for(cp *p=a;p!=a+n;p+=l){ cp w(1,0); for(int i=0;i<l/2;i++){ cp t=w*p[i+m]; p[i+m]=p[i]-t; p[i]=p[i]+t; w=w*wn; } } } if(f==-1){for(int i=0;i<n;i++)a[i].x/=n;}// } int main(){ scanf("%d%d",&n,&m);n++;m++; for(int i=0;i<n;i++){int u;scanf("%d",&u);a[i]=cp(u,0);} for(int i=0;i<m;i++){int u;scanf("%d",&u);b[i]=cp(u,0);} int N=1;while(N<n+m)N*=2; fft(a,N,1);fft(b,N,1); for(int i=0;i<N;i++)a[i]=a[i]*b[i]; fft(a,N,-1); for(int i=0;i<n+m-1;i++)printf("%d ",(int)(a[i].x+0.1)); return 0; }
注意：
1.数组空间必须是大于两个卷积数组长度和的2的幂（记为N），当n为1e5时数组空间为300000。
2.DFT后的操作一定要以N为单位，例如点值相乘。这个问题在多重fft的题目中很容易写错，多重fft必须把上一个的结果后半部分清零再继续。
3.代码需要手写复数模板来减小常数，NTT还要预处理omega。枚举反二进制的方法是从高位开始模拟进位。
（八）卷积
对于函数f(n)和g(n)，定义其卷积为函数(f⊗g)
(fg)(n)=Σf(i)g(n-i)，i=0~n。——形式幂级数
卷积的形式和多项式乘法类似，(fg)的生成函数是f(n)和g(n)的生成函数的乘积。
卷积是和为定值的形式，若差为定值，将其中一个数组反转后即可卷积。卷积中多余的部分数组为0，不影响答案。
高精度乘法：将数字从低位到高位编号0~len-1，每一位代入多项式系数，那么数字乘法就是多项式乘法，最后从低到高处理进位。
（九）模数任意——fft合并
将每个数字拆成a*32768+b，然后四次dft后4次idft合并。
#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; }
注意清空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)
（二）快速数论变换
当模数为形如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)。
如果觉得集合很难理解，不妨用”选择项“这个词。
题意：给定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。
题意：有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 拉格朗日插值（自然数幂和）
