多项式可以加减乘,这是基本常识。然而,多项式还支持除法、乘方、开方、求对数等等运算,这可实在是有些神奇。于是在此小结一下多项式的诸多运算。在接下来的叙述中,定义多项式。
为方便叙述,在下列所有运算中,多项式的运算均定义在意义下,整数的运算均定义在
意义下,其中
是一个大质数。在程序实现中,取
。
加法/减法
给定多项式
,求
或
。
将多项式的对应项直接相加减即可。时间复杂度。
1 2 3 4 5 6 7 8 9 |
//输入:A,B,n,m 输出:C void add(LL* A,LL* B,int n,int m,LL* C) { for(int i=0;i<max(n,m);i++) C[i]=(A[i]+B[i])%mo; } void sub(LL* A,LL* B,int n,int m,LL* C) { for(int i=0;i<max(n,m);i++) C[i]=(A[i]-B[i]+mo)%mo; } |
乘法
给定多项式
,求
。
对于普通的多项式乘法,可以用经典的FFT或分治乘法实现。在模特殊的大质数时可以使用NTT,一般要求
,如UOJ模数
。时间复杂度
。
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 |
const int maxn=(int)(1e5)+5; const double pi=acos(-1.0); struct cpx { double x,y; cpx(double a=0,double b=0):x(a),y(b) {} }; cpx operator + (cpx a,cpx b) {return cpx(a.x+b.x,a.y+b.y);} cpx operator - (cpx a,cpx b) {return cpx(a.x-b.x,a.y-b.y);} cpx operator * (cpx a,cpx b) {return cpx(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);} void DFT(cpx* a,int n,int d=1) { for(int i=(n>>1),j=1;j<n;j++) { if(i<j) swap(a[i],a[j]); int k;for(k=(n>>1);i&k;i^=k,k>>=1);i^=k; } for(int m=2;m<=n;m<<=1) { cpx w=cpx(cos(2*pi/m*d),sin(2*pi/m*d)); for(int i=0;i<n;i+=m) { cpx s=cpx(1,0); for(int j=i;j<(i+(m>>1));j++) { cpx u=a[j],v=a[j+(m>>1)]*s; a[j]=u+v;a[j+(m>>1)]=u-v; s=s*w; } } } if(d==-1) for(int i=0;i<n;i++) a[i].x/=n; } //输入:A,B,n,m 输出:C void mul(cpx* A,cpx* B,int n,int m,cpx* C) { int len=0; while((1<<len) < max(n,m)<<1) len++; DFT(A,1<<len);DFT(B,1<<len); for(int i=0;i<(1<<len);i++) C[i]=A[i]*B[i]; DFT(A,1<<len,-1);DFT(B,1<<len,-1);DFT(C,1<<len,-1); } |
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 |
const int maxn=(int)(1e5)+5; const int mo=998244353; const int ge=3; LL pow_mod(LL b,LL p,LL k) { LL res=1; for(;p;p>>=1) { if(p&1) res=res*b%k; b=b*b%k; } return res; } void NTT(LL* a,int n,int d=1) { for(int i=(n>>1),j=1;j<n;j++) { if(i<j) swap(a[i],a[j]); int k;for(k=(n>>1);i&k;i^=k,k>>=1);i^=k; } for(int m=2;m<=n;m<<=1) { LL w=pow_mod(ge,(mo-1+d*mo/m)%(mo-1),mo); for(int i=0;i<n;i+=m) { LL s=1; for(int j=i;j<(i+(m>>1));j++) { LL u=a[j],v=a[j+(m>>1)]*s%mo; a[j]=(u+v)%mo;a[j+(m>>1)]=(u+mo-v)%mo; s=s*w%mo; } } } LL xinv=pow_mod(n,mo-2,mo); if(d==-1) for(int i=0;i<n;i++) a[i]=a[i]*xinv%mo; } //输入:A,B,n,m 输出:C void mul(LL* A,LL* B,int n,int m,LL* C) { int len=0; while((1<<len) < max(n,m)<<1) len++; NTT(A,1<<len);NTT(B,1<<len); for(int i=0;i<(1<<len);i++) C[i]=A[i]*B[i]; NTT(A,1<<len,-1);NTT(B,1<<len,-1);NTT(C,1<<len,-1); for(int i=n+m;i<(1<<len);i++) C[i]=0; } |
求导/积分
给定多项式
,求
或
。
根据定义逐项求导/积分即可,具体过程见微积分基本定理。时间复杂度。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
//输入:A,n 输出:D/I LL pow_mod(LL b,LL p,LL k) { LL res=1; for(;p;p>>=1) { if(p&1) res=res*b%k; b=b*b%k; } return res; } void derivation(LL* A,int n,LL* D) { for(int i=0;i<n-1;i++) D[i]=A[i+1]*(i+1)%mo; } void integration(LL* A,int n,LL* I) { I[0]=1; for(int i=1;i<n+1;i++) I[i]=A[i-1]*pow_mod(i,mo-2,mo)%mo; } |
多项式的一些比较正常的运算就到此结束了。接下来是一些有趣的姿势。
逆元
给定多项式
,求
。
称为
的逆元,满足
。
首先考虑,显然有
。
考虑倍增。如果我们知道了,怎样得到
呢?
设,显然有
则有
可以推得
上述过程可以用FFT或NTT实现。如此,就可以在次迭代后求出多项式的逆元。根据主定理有
于是时间复杂度是。在此膜拜一下Miskcoo神犇的Blog。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
const int maxn=(int)(1e5)+5; const int mo=998244353; int stk[maxn],top; //输入:A,n 输出:B 辅助空间:T void inv(LL* A,int n,LL* B,LL* T) { int tot=n; for(top=0;tot>1;stk[++top]=tot,tot=(tot+1)/2); B[0]=pow_mod(A[0],mo-2,mo); for(int i=top;i>=1;i--) { int k=stk[i]; int len=0; while((1<<len) < (k<<1)) len++; for(int i=0;i<k;i++) T[i]=A[i]; for(int i=k;i<(1<<len);i++) T[i]=0; NTT(T,1<<len);NTT(B,1<<len); for(int i=0;i<(1<<len);i++) B[i]=(2+mo-(T[i]*B[i]%mo))*B[i]%mo; NTT(B,1<<len,-1); for(int i=k;i<(1<<len);i++) B[i]=0; } } |
除法/模
给定
次多项式
,
次多项式
,求
。
//似乎没有这么简单....我重更一下。
题目所求即为,我们只需要求
的逆元即可。时间复杂度
。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
const int maxn=(int)(1e5)+5; LL rA[maxn],rB[maxn],X[maxn],T[maxn]; //输入:A,B,n,m 输出:D,R void div_mod(LL* A,LL* B,int n,int m,LL* D,LL* R) { for(int i=0;i<n;i++) rA[i]=A[n-i-1]; for(int i=0;i<m;i++) rB[i]=B[m-i-1]; inv(rB,n-m+1,X,T);mul(rA,X,n,m,D); for(int i=n-m+1;i<n+m;i++) D[i]=0; for(int i=0;i<n-m-i;i++) swap(D[i],D[n-m-i]); mul(B,D,n,n-m+1,X);sub(A,X,m-1,m-1,R); for(int i=n-m+1;i<n+m-1;i++) D[i]=0; } |
对数
给定多项式
,求
。
设,根据链式法则有
于是只需要求一次的逆元,与
相乘得到
,再积分还原
即可。时间复杂度
。
注意,在求对数时,必须保证。
指数
给定多项式
,求
,即
。
设。根据定义,可以列出方程
。
考虑牛顿迭代法。首先,易知。
接下来我们考虑怎样由推出
。
显然有。在
处泰勒展开,可得
则有
注意到泰勒展开时,由于,大于等于二次的项均可以被消去,于是剩下的式子就十分简洁明了。将方程
代入上式,得
根据主定理有,于是时间复杂度
。
乘方
给定多项式
,求
的前
次项。
方法1:快速幂
考虑像一般整数的快速幂一样,每一轮迭代FFT一次并舍去多余项。时间复杂度。
方法2:指对变换
由指数函数和对数函数的性质,有
于是时间复杂度就是。注意,在这种方法中,由于需要对多项式求对数,则须
。若
怎么办呢?我们可以稍微转化一下。设
的最低次项为
,则有
此时括号内的式子一定满足,满足对数运算条件。
开方
给定多项式
,求
。
类似乘方运算的方法2,我们有
时间复杂度同样是。对于常数项不为1的情况,同样有
注意到必须满足才可以。
文中许多思路和推导过程摘自JCVB的论文《生成函数的运算和组合计数问题》。