多项式的操作大赏

该文章已整合到多项式全家桶，模板封装，并且修改了很多错误，笔芯~

---

写在前面

         \ \ \ \ \ \ \,       今年WC，t1就考了多项式呢，是laofu出的，当然了我这么才并没有做出来，骗了36分就跑路了。

         \ \ \ \ \ \ \,       不过还是看出来是考的多项式了，记得laofu在讲课的时候，前半场还知道怎么操作，后面就完全不知道了。

         \ \ \ \ \ \ \,       比如说exp啊，不仅不会还不知道是什么东西呢。

         \ \ \ \ \ \ \,       因为以前考的多项式比较少，就没怎么注意去学，就学了求逆和卷积这种简单的。于是后来就找时间学了一下感觉比较重要的操作。就以洛咕板子为目录，在此记录一下。感觉多项式还是挺好理解的，就是要用很多所学知识，还有临时变量名没有清干净的话会特别麻烦。

         \ \ \ \ \ \ \,       那就这样开始吧，我会贴上一部分我丑陋的代码。

---

• 乘法(卷积)

         \ \ \ \ \ \ \,       这个东西嘛，说实话吧，没有那么好理解，但是却很基础很重要，算是一个大前提吧，具体一点的我写在了这里【求多项式卷积的变换】。

         \ \ \ \ \ \ \,       那么这里就贴一下板子，下面的模板会用到的这个样子，因为很多时候都会用到逆元，所以说多项式的题，NTT见得多：

const int mod=998244353,mod_g=3;
int R[N];
int power(int a,int b){
  	int ans=1;
  	for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;
  	return ans;
}
#define Inv(x) power(x,mod-2)
int Polynomial_init(int n){
	int len;for(len=1;len<=n;len<<=1);
	return len;
}
void NTT(int *a,int f,int la){
	int n=la;
   	for(int i=0;i<n;i++)if(i<R[i])swap(a[i],a[R[i]]);
   	for(int i=1;i<n;i<<=1){
    	int gn=power(mod_g,(mod-1)/(i<<1));
    	for(int j=0;j<n;j+=(i<<1)){
      		int g=1;
      		for(int k=0;k<i;k++,g=1ll*g*gn%mod){
        		int x=a[j+k],y=1ll*g*a[i+j+k]%mod;
        		a[j+k]=(x+y)%mod;a[i+j+k]=(x-y+mod)%mod;
      		}
    	} 
  	}
  	if(f==-1){
   		reverse(a+1,a+n);
    	int inv=Inv(n);
    	for(int i=0;i<n;i++)a[i]=1ll*a[i]*inv%mod;
  	}
}

         \ \ \ \ \ \ \,       然后会用下面两种方式用到卷积：

         \ \ \ \ \ \ \,       $Convolution$函数是把 $a$，$b$两个多项式一起卷到$a$里，会破坏 $a$，$b$原来的形态。

int Convolution(int *a,int *b,int la,int lb){
	int n=la,m=lb;
	int L=0;for(m+=n,n=1;n<=m;n<<=1)L++;
	for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	NTT(a,1,n);NTT(b,1,n);
	for(int i=0;i<=n;i++)a[i]=1ll*a[i]*b[i]%mod;
	NTT(a,-1,n);
	return m;
}

         \ \ \ \ \ \ \,       $Multiplication$函数是把 $a$，$b$两个多项式一起另外卷到$c$里，不会破坏 $a$，$b$原来的形态。

int X[N],Y[N];
void Multiplication(int *a,int *b,int *c,int la,int lb){
  	int n=la,m=lb;
  	int L=0;for(m+=n,n=1;n<=m;n<<=1)L++;
  	for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	for(int i=0;i<=la;++i) X[i]=a[i];
	for(int i=0;i<=lb;++i) Y[i]=b[i];
  	NTT(X,1,n);NTT(Y,1,n);
  	for(int i=0;i<=n;i++)c[i]=1ll*X[i]*Y[i]%mod,X[i]=Y[i]=0;
  	NTT(c,-1,n);
}

---

• 求逆

         \ \ \ \ \ \ \,       求一个多项式 $A$ 的模 $x^n$ 的逆元$B$ 时，假设先求出了模$x^{\frac{n}{2}}$ 的逆元 $B_2$，既：

$$A\ast B_2\equiv 1(mod{x}^{\frac{n}{2}})$$

$$A\ast B\equiv 1(mod{x}^n)$$

         \ \ \ \ \ \ \,       那么显然存在：

$$A\ast B\equiv 1(mod{x}^{\frac{n}{2}})=A\ast B_2$$

$$B-B_2\equiv 0(mod{x}^{\frac{n}{2}})$$

         \ \ \ \ \ \ \,       两边同时平方：

$$B^2-2B{B}_2+B_2^2\equiv 0(mod{x}^n)$$

         \ \ \ \ \ \ \,       再把 $A$ 乘回去：

$$(A\ast B)\ast B-(A\ast B)\ast 2B_2+A\ast B_2^2\equiv A\ast 0(mod{x}^n)$$

$$B-2B_2+A\ast B^2\equiv 0(mod{x}^n)$$

$$B\equiv 2B_2-A\ast B_2^2(mod{x}^n)$$

         \ \ \ \ \ \ \,       我们就可以递归来处理他了，边界是$B_0\equiv A_0^{-1}(mod{x}^1)$

         \ \ \ \ \ \ \,       这里代码是把 $a$ 求逆到 $b$,不会破坏 $a$ 原来的形态，其中的 $len$ 是 $a$ 的长度在二进制下向上取的整数：

int C[N];
void Inverse(int *a,int *b,int len){
  	if(len==1){b[0]=Inv(a[0]);return;}
  	Inverse(a,b,(len+1)>>1);
  	int L=0,n=1;
  	for(;n<(len<<1);n<<=1)L++;
  	for(int i=1;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
  	for(int i=0;i<len;i++)C[i]=a[i];
  	for(int i=len;i<n;i++)C[i]=0;
  	NTT(C,1,n);NTT(b,1,n);
  	for(int i=0;i<=n;i++)b[i]=1ll*(2ll-1ll*C[i]*b[i]%mod+mod)%mod*b[i]%mod;
  	NTT(b,-1,n);
  	for(int i=len;i<n;i++)b[i]=0;
}

---

• 除法

         \ \ \ \ \ \ \,       会了乘法和逆元当然就会除法啦，就是套板子的事情了呢。

         \ \ \ \ \ \ \,       我们求多项式 $A$ 除以 $B$ ，那么商 $Q$就等于：

$$A\ast Inverse(B)$$

         \ \ \ \ \ \ \,       余 $R$ 呢就等于：

$$A-(G\ast B)$$

         \ \ \ \ \ \ \,       很简单吧，就不解释了，这里模板指的是$a÷b$，$la$是$a$的长度，$lb$是$b$的长度：

int H[N];
void Division(int *a,int *b,int la,int lb,int *Quotient,int *Remainder){
	int n=la,m=lb;
	reverse(a,a+1+n),reverse(b,b+1+m);
  	int len=1;while(len<=n-m)len<<=1;
  	memset(C,0,sizeof(C));
  	Inverse(b,H,len);
  	Multiplication(a,H,Quotient,n-m,n-m);
  	reverse(Quotient,Quotient+n-m+1);
  	reverse(a,a+n+1),reverse(b,b+m+1);
  	Multiplication(Quotient,b,Remainder,n-m,m);
  	for(int i=0;i<m;++i)Remainder[i]=(a[i]-Remainder[i]+mod)%mod;
}

         \ \ \ \ \ \ \,       顺带一提，洛咕这道模板题还是蛮恶心的，主要是变量名太多了，调起来巨麻烦。

下面来两个简单的，求导和求积分

---

• 求导

         \ \ \ \ \ \ \,       按照公式来，设多项式 $A$ 的导数为 $A^{\prime }$。

         \ \ \ \ \ \ \,       那么有：

$$x^{A^{\prime }}=A{x}^{A-1}$$

         \ \ \ \ \ \ \,       既：

$$A_i^{\prime }=i\times A_{i+1}$$

         \ \ \ \ \ \ \,       模板是把 $a$ 求导到 $b$，其中的 $len$ 是 $a$ 的长度在二进制下向上取的整数。

void Derivation(int *a,int *b,int n){
	for(int i=1;i<n;i++)
  	b[i-1]=1ll*i*a[i]%mod;
	b[n-1]=0;
}

---

• 求积分

         \ \ \ \ \ \ \,       还是按照公式来，设多项式 $A$ 的积分为 $A^{\prime }$。

         \ \ \ \ \ \ \,       那么有：

$$\int x^{A^{\prime }}dx=\frac{1}{A+1}{x}^{A-1}$$

         \ \ \ \ \ \ \,       既：

$$A_i^{\prime }=\frac{A_{i-1}}{i}$$

         \ \ \ \ \ \ \,       模板还是把 $a$ 求积分到 $b$，其中的 $len$ 还是 $a$ 的长度在二进制下向上取的整数。

void Integral(int *a,int *b,int n){
	for(int i=1;i<n;i++)
  	b[i]=1ll*Inv(i)*a[i-1]%mod;
	b[0]=0;
}

然后看看这两个在求对数上的实际应用吧

---

• 对数

         \ \ \ \ \ \ \,       嘛，微积分这种东西要我讲我肯定也不懂的，但是我知道公式呀：

$$\ln (A)=\int \frac{A^{\prime }}{A}dx$$

         \ \ \ \ \ \ \,       就是它的逆卷上它的导，然后求个积分，基本上前面说的全都用上了呢，模板还是把 $a$ 求对数到 $b$，$len$ 还是 $a$ 的长度在二进制下向上取的整数：

int A[N],B[N];
void Logarithmic(int *a,int *b,int len){
	memset(A,0,sizeof(A));
	memset(B,0,sizeof(B));
  	Derivation(a,A,len);
  	memset(C,0,sizeof(C));
  	Inverse(a,B,len);
  	Convolution(A,B,len,len);
  	Integral(A,b,len);
}

---

• 指数

         \ \ \ \ \ \ \,       这玩意说简单也简单，说复杂也挺复杂的，我们得引入一个新玩意，【牛顿迭代】，我在这里有简略的讲过：【导数和牛顿迭代】，我也在洛谷出过一个牛顿迭代的裸题，感兴趣可以去体验一下牛顿迭代的神奇：【P4986 逃离】。

         \ \ \ \ \ \ \,       说远了，现在我们来康康怎么求指数函数。

         \ \ \ \ \ \ \,       令我们要求的是 $A$ 的指数函数 $B$，既是：

$$B=e^A$$

         \ \ \ \ \ \ \,       变形得：

$$\ln (B)-A=0$$

         \ \ \ \ \ \ \,       哇，函数零点，哭了啊，马上牛顿迭代呀！

         \ \ \ \ \ \ \,       不过这种套多项式的函数怎么牛顿迭代啊？我们不妨先写出来化简康康，（其实多项式也可以看做数值呀）。

         \ \ \ \ \ \ \,       设$F(B)=\ln (B)-A$，我们要求 $F$的零点，根据牛顿迭代的公式可得（注意这里B后面的括号的迭代版本的意思，不是多项式的项）：

$$B(x)=B(x-1)-\frac{F\left(B(x-1)\right)}{F^{\prime }\left(B(x-1)\right)}$$

         \ \ \ \ \ \ \,       而根据导数的定义，$F^{\prime }(B)=\frac{1}{B}$，（【导数和牛顿迭代】里面有提到过一点，这里是把$A$当做常数了）

         \ \ \ \ \ \ \,       那我们现在把牛顿迭代的公式化简：

$$\begin{array}{rl}B(x)& =B(x-1)-B(x-1)\ast F\left(B(x-1)\right)\\ & =B(x-1)-B(x-1)\ast \ln \left(B(x-1)\right)-A\\ & =B(x-1)\ast \left(1-\ln \left(B(x-1)\right)-A\right)\end{array}$$

         \ \ \ \ \ \ \,       再次强调B后面的括号的迭代版本的意思，不是多项式的项。

         \ \ \ \ \ \ \,       现在看似分治可做，我们用两个容器相互装版本。每次老版本一卷，新版本长度就会倍增，所以说我们做 $\log n$ 次迭代就好。我们的操作相当于把式子拆了求收敛值，所以不会有精度的问题，求出来就好了。

         \ \ \ \ \ \ \,       那么现在的问题是，第一个版本怎么样，洛咕模板给的是保证 $A_0=0$，所以说 $B_0=1$，也就是常数项为 $1$。

         \ \ \ \ \ \ \,       看到这里感觉很难很绕，不过没关系，模板挺短的，光记模板也可以吧。

         \ \ \ \ \ \ \,       模板还是把 $a$ 求指数到 $b$，$len$ 还是 $a$ 的长度在二进制下向上取的整数：

int D[N];
void Exponential(int *a,int *b,int len){
  	if(len==1){b[0]=1;return;}
  	Exponential(a,b,len>>1),Logarithmic(b,D,len);
  	D[0]=(1ll*a[0]+1ll-D[0]+mod)%mod;
  	for(int i=1;i<len;++i) D[i]=(1ll*a[i]-D[i]+mod)%mod;
  	Convolution(b,D,len<<1,len<<1);
  	for(int i=len;i<(len<<1);++i) b[i]=D[i]=0;
}

有了指数和对数这两个工具，多项式可以干的事情就太多了，比如说开$k$次方根，快速幂什么的。

---

• 开根

         \ \ \ \ \ \ \,       我们知道对于一个数 $a^b$ 开$k$次方
$$\sqrt[k]{a^b}=a^{\frac{b}{k}}$$

         \ \ \ \ \ \ \,       那么对于一个多项式 $A$ 开 $k$ 次方也同理：

$$\sqrt[k]{A}=e^{\frac{\ln (A)}{k}}$$

         \ \ \ \ \ \ \,       好像有更快的方法，不过这个方法挺显然简单的，模板也很简单了，模板还是把 $a$ 求指数到 $b$，$len$ 还是 $a$ 的长度在二进制下向上取的整数 （说吐了）：

         \ \ \ \ \ \ \,       对了，求 $EXP$ 时，常数项还需注意单独求，复杂的话会涉及高次剩余，就不多赘述了，$solve(a,K)$ 函数就是求高次剩余的，表示$a$在模$mod$的情况下的$k$次剩余，也可以BSGS水，感兴趣可以看看这里：题解 P5277 【【模板】多项式开根（加强版）】

int E[N];
void Kth_root(int *a,int *b,int len,int k){
	Logarithmic(a,E,len);
	int bot=Inv(k);
	for(int i=1;i<=len;i++)E[i]=1ll*E[i]*bot%mod;
	b[0]=solve(a[0],K);
	Exponential(E,b,len);
}

---

• 快速幂

         \ \ \ \ \ \ \,       快速幂和开根同理啦，直接给代码咯，就是注意求 $EXP$ 时，常数项不是 $1$ 了，注意改一改，为什么可以这样搞我放在这里了，有兴趣可以看一看： 题解 P5273 【【模板】多项式幂函数 (加强版)】

int P[N];
void Fast_Power(int *a,int *b,int len,int k){
	Logarithmic(a,P,len);
	for(int i=1;i<=len;i++)P[i]=1ll*E[i]*k%mod;
  	b[0]=power(a[0],k);
	Exponential(P,b,len);
}

---

• 多项式三角函数

         \ \ \ \ \ \ \,       先有一个神奇的欧拉公式：

$$e^{ix}=\cos x+i\cdot \sin x$$

         \ \ \ \ \ \ \,       容易变换得到：

$$\sin x=\frac{e^{ix}-e^{-ix}}{2i}$$

$$\cos x=\frac{e^{ix}+e^{-ix}}{2}$$

         \ \ \ \ \ \ \,       那么对于一个多项式 $A$：

$$\sin x=\frac{e^{i\cdot A}-e^{-i\cdot A}}{2i}$$

$$\cos x=\frac{e^{i\cdot A}+e^{-i\cdot A}}{2}$$

         \ \ \ \ \ \ \,       然后就是多项式求指数函数了，贴板，问题不大……

         \ \ \ \ \ \ \,       诶等等，这个 $i$ 怎么办啊？怎么取模？

         \ \ \ \ \ \ \,       不慌，我们来回忆一下 $i$ 是个什么东西：

$$i=\sqrt{-1}$$

$$i^2=-1$$

         \ \ \ \ \ \ \,       所以说，在模的意义下：

$$i^2\equiv mod-1(\%mod)$$

         \ \ \ \ \ \ \,       哇！二次剩余诶！（可是我不会，死了）

         \ \ \ \ \ \ \,       有原根的模数，还需要求二次剩余吗？

         \ \ \ \ \ \ \,       当然是$BSGS$乱搞啊：

         \ \ \ \ \ \ \,       令已知原根是$g$，设 $i=g^x$ 所以原式化为:

$$g^{2x}\%mod=mod-1$$

         \ \ \ \ \ \ \,       $BSGS$即可求出 $2x$ 的取值，便可以求出 $i$ 。

         \ \ \ \ \ \ \,       不过在模数是$998244353$的情况下，我们已经算出来是下面这样的，可以直接用了：

const int img=86583718;

         \ \ \ \ \ \ \,       然后就是贴板子了，$a$,$b$,$len$的意义同上，不重复说了,$tan$的话，再加个求逆，不说了：

int Sin_x[N],Sin_y[N],Sin_z[N];
void Sin(int *a,int *b,int len){
	for(int i=0;i<=len;i++)Sin_x[i]=1ll*img*a[i]%mod;
  	memset(Exp_x,0,sizeof(Exp_x));
	Exponential(Sin_x,Sin_y,len);
	for(int i=0;i<=len;i++)Sin_x[i]=(1ll*mod-(1ll*img*a[i]%mod))%mod;
  	memset(Exp_x,0,sizeof(Exp_x));
	Exponential(Sin_x,Sin_z,len);
	int bot=Inv(2ll*img%mod);
	for(int i=0;i<=len;i++)b[i]=1ll*(Sin_y[i]-Sin_z[i]+mod)%mod*bot%mod;
}
int Cos_x[N],Cos_y[N],Cos_z[N];
void Cos(int *a,int *b,int len){
	for(int i=0;i<=len;i++)Cos_x[i]=1ll*img*a[i]%mod;
  	memset(Exp_x,0,sizeof(Exp_x));
	Exponential(Cos_x,Cos_y,len);
	for(int i=0;i<=len;i++)Cos_x[i]=(1ll*mod-(1ll*img*a[i]%mod))%mod;
  	memset(Exp_x,0,sizeof(Exp_x));
	Exponential(Cos_x,Cos_z,len);
	int bot=Inv(2);
	for(int i=0;i<=len;i++)b[i]=1ll*(Cos_y[i]+Cos_z[i])%mod*bot%mod;
}

---

• 多项式反三角函数

         \ \ \ \ \ \ \,       嘛，这个还是套模板呢，先给公式吧，一般也搜得到下面公式的推导，是关于多项式$A$的三角函数：

         \ \ \ \ \ \ \,       根据三角函数求导公式：

$$\frac{d}{dx}asin(x)=\frac{1}{\sqrt{-x^2+1}}$$

$$\frac{d}{dx}atan(x)=\frac{1}{x^2+1}$$

         \ \ \ \ \ \ \,       所以显然有

$$asin(A)=\int \frac{A^{\prime }}{\sqrt{-A^2+1}}dx$$

$$atan(A)=\int \frac{A^{\prime }}{A^2+1}dx$$

         \ \ \ \ \ \ \,       套模板：$a$,$b$,$len$的意义同上：

int Asin_x[N],Asin_y[N],Asin_z[N];
void Asin(int *a,int *b,int len){
  	Derivation(a,Asin_x,len);
	int n=len,m=len;
  	int L=0;for(m+=n,n=1;n<=m;n<<=1)L++;
  	for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	for(int i=0;i<=len;++i) Asin_y[i]=a[i];
  	NTT(Asin_y,1,n);
  	for(int i=0;i<=n;i++)Asin_y[i]=1ll*Asin_y[i]*Asin_y[i]%mod;
  	NTT(Asin_y,-1,n);
  	for(int i=0;i<=n;i++)Asin_y[i]=(mod-Asin_y[i]);
	Asin_y[0]=(1+Asin_y[0])%mod;
  	Kth_root(Asin_y,Asin_z,len,2);
  	memset(Asin_y,0,sizeof(Asin_y));
  	Inverse(Asin_z,Asin_y,len);
  	Convolution(Asin_x,Asin_y,len,len);
  	Integral(Asin_x,b,len);
}
int Atan_x[N],Atan_y[N],Atan_z[N];
void Atan(int *a,int *b,int len){
	Derivation(a,Atan_x,len);
	int n=len,m=len;
  	int L=0;for(m+=n,n=1;n<=m;n<<=1)L++;
  	for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	for(int i=0;i<=len;++i) Atan_y[i]=a[i];
  	NTT(Atan_y,1,n);
  	for(int i=0;i<=n;i++)Atan_y[i]=1ll*Atan_y[i]*Atan_y[i]%mod;
  	NTT(Atan_y,-1,n);
	Atan_y[0]=(1+Atan_y[0])%mod;
  	Inverse(Atan_y,Atan_z,len);
  	Convolution(Atan_x,Atan_z,len,len);
  	Integral(Atan_x,b,len);
}

         \ \ \ \ \ \ \,        你问我 $acos$ ? 大同小异就不放了 （其实是洛咕没有这板子）

---

         \ \ \ \ \ \ \,       差不多基础操作就这么多了吧。多点求值的话，正解两个$\log$还带大常数，周道觉得不如写暴力卡卡常，快速插值呢，没有计划去学，拉格朗日差值可以做到$O(n^2)$的复杂度，有兴趣可以去看看【拉格朗日插值法】。虽然感觉很多时候还不够用，不过目前应该不会怎么遇到差值的题吧（奶。

         \ \ \ \ \ \ \,       明天找时间开这个东西：#150. 挑战多项式

         \ \ \ \ \ \ \,        2.22 Update：太恶心了，放弃了，有机会找时间慢慢调。

         \ \ \ \ \ \ \,        4.1 Update：周道终于写出来了！！！！

【没有结构体，没有封装，大常数预警（代码链接）】

         \ \ \ \ \ \ \,       截止周道提交时间，周道喜提全网最慢代码，全网第三大空间代码！