多项式全家桶——Part.3 多项式求逆、除法、开根号

多项式全家桶正式进入正片。

有一说一，感觉之前写的其实都是写不重要或特别基础的东东。
现在才是一些比较有难度的东东。

多项式求逆

介绍

多项式求逆即为求一个多项式的乘法逆元。
它长这样：
$$A(x)\ast B(x)\equiv 1(mod{x}^n)$$
然后多项式$B(x)$就叫做多项式$A(x)$的逆元。
这玩意和普通的乘法逆元有一点不同的就是，它的模是一个$x^n$，其意义为把结果中次数大于等于n的项都给忽略掉。

过程

引理一： 考虑一个乘法逆元满足：$A(x)\ast B(x)\equiv 1(mod{x}^n)$，那么我们可以得到另一个柿子：$A(x)\ast B(x)\equiv 1(mod{x}^m)$其中$1<=m<=n$
证明？把卷积形式写出来就知道了。

考虑设$B(x)$表示$A(x)$在模$x^n$意义下的逆元。
设$G(x)$表示$A(x)$在模$x^{\lceil \frac{n}{2}\rceil }$意义下的逆元。
由引理得：
$$\begin{gathered} B(x)\equiv G(x)(mod{x}^{\lceil \frac{n}{2}\rceil }) \\ B(x)-G(x)\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil }) \end{gathered}$$
开方（这一步我也不能严谨证明，但是打表发现这是对的）
$$B(x{)}^2+G(x{)}^2-2\ast B(x)\ast G(x)\equiv 0(mod{x}^n)$$
乘个$A(x)$进去，然后由于是模$x^n$意义下，所以可以消去一些$B(x)$:
$$\begin{gathered} B(x)+A(x)\ast G(x{)}^2-2\ast G(x)\equiv 0(mod{x}^n) \\ B(x)\equiv G(x)(2-A(x)\ast G(x))(mod{x}^n) \end{gathered}$$
这时，就可以用$G(x)$来求出$B(x)$。具体过程就递归求解即可（当然不觉得麻烦的话打分治FFT也是可以），中间需要用到NTT。

学习资料：
https://blog.csdn.net/a_forever_dream/article/details/102483602
https://www.cnblogs.com/LanrTabe/p/11359952.html
https://www.cnblogs.com/guoshaoyang/p/11296025.html

代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
const long long mo=998244353;
const int maxn=300010;

long long a[maxn],b[maxn],rev,A[maxn],G[maxn];
int jl,n,rec[maxn],w[maxn],up,dep;

long long qsm(long long a,long long b)
{
	long long t=1;
	long long y=a;
	while (b>0)
	{
		if ((b&1)==1) t=t*y%mo;
		y=y*y%mo;
		b/=2;
	}
	return t;
}

void FFT(long long a[],int n,int inv)
{
	if (n==1) return;
	int mid=n/2;
	for (int i=0;i<n;i++)
	{
		rec[i]=(rec[i>>1]>>1)|((i&1)<<(dep-1));
		if (i<rec[i]) swap(a[i],a[rec[i]]);
	}
	for (int len=2;len<=n;len<<=1)
	{
		int mid=len/2;
		for (int j=0;j<mid;j++)
		{
			for (int k=j;k<n;k+=len)
			{
				int p,q;
				q=w[j*(n/len)]*a[k+mid]%mo,p=a[k];
				a[k+mid]=(p-q+mo)%mo;
				a[k]=(p+q+mo)%mo;
			}
		}
	}
}

void solve(int n)
{
	if (n==1)
	{
		b[0]=qsm(a[0],mo-2);
	}
	else
	{
		double op=n;
		int oq=ceil(op/2);
		solve(oq);
		
		up=1;dep=0;
		while (up<=n+jl) up=up*2,dep++;
		w[0]=1;
		rev=qsm(3,(mo-1)/up);
		for (int i=1;i<=up;i++) w[i]=w[i-1]*rev%mo;
		for (int i=0;i<up;i++)
		{
			G[i]=0;
			A[i]=a[i];
			if (i<oq) G[i]=b[i];
		}
		FFT(A,up,1);
		FFT(G,up,1);
		for (int i=0;i<up;i++)
		{
			A[i]=G[i]*(2-A[i]*G[i]%mo+mo)%mo;
		}
		for (int i=0;i<=up/2;i++)
		{
			swap(w[i],w[up-i]);
		}
		FFT(A,up,-1);
		for (int i=0;i<n;i++) b[i]=A[i]*qsm(up,mo-2)%mo;
		int kk=0;
	}
}

int main()
{
	scanf("%d",&n);
	jl=n;
	for (int i=0;i<n;i++)
	{
		scanf("%lld",&a[i]);
	}
	solve(n);
	for (int i=0;i<n;i++)
	{
		printf("%lld ",b[i]);
	}
	printf("\n");
}

多项式除法

介绍

这玩意是求一个长成这样的柿子的：
$$A(x)=B(x)\ast D(x)+R(x)$$
其中$A(x)$和$B(x)$已经给出，且$A(x)$次数为n，$B(x)$次数为m，求$D(x)$和$R(x)$
也就是除法运算中求商和余数，话说如果上过初中就应该明白长除法。

过程

引理二：
考虑一个反转操作为将一个多项式的系数调换，用数学形式表达即为：
$$Pr(x)=x^n\ast P(\frac{1}{x})=\sum _{i=0}^n{a}_{n-i}\ast x^i$$
证明显然。

那么现在就推柿子：
$$\begin{gathered} A(x)=B(x)D(x)+R(x) \\ A(\frac{1}{x})=B(\frac{1}{x})D(\frac{1}{x})+R(\frac{1}{x}) \\ {x}^{n}A(\frac{1}{x})=x^{n}B(\frac{1}{x})D(\frac{1}{x})+x^{n}R(x) \\ {x}^{n}A(\frac{1}{x})=x^{m}B(\frac{1}{x})x^{n-m}D(\frac{1}{x})+x^{n-m+1}{x}^{m-1}R(x) \end{gathered}$$
然后再引入引理：
$$Ar(x)=Br(x)Dr(x)+x^{n-m+1}Rr(x)$$
由于我们的$Dr(x)$的次数是小于等于$n-m$的，所以如果对原式取模只要大于之即可，而$A(x)$、$B(x)$都已知，所以不必理会。
$$Ar(x)\equiv Br(x)Dr(x)(mod{x}^{n-m+1})$$
然后我们就可以直接求$Br(x)$的逆元，再乘$Ar(x)$即可得到$Dr(x)$，再把系数反过来即可求出$D(x)$。

当然，求出$D(x)$之后，再求$R(x)$就很简单了。

学习资料
https://www.cnblogs.com/knife-rose/p/12056617.html
https://www.cnblogs.com/owenyu/p/6724611.html

代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
const long long mo=998244353;
const int maxn=300010;

long long a[maxn],b[maxn],ar[maxn],br[maxn],d[maxn],D[maxn],r[maxn],rev,A[maxn],G[maxn];
int jl,n,m,rec[maxn],w[maxn],up,dep;

long long qsm(long long a,long long b)
{
	long long t=1;
	long long y=a;
	while (b>0)
	{
		if ((b&1)==1) t=t*y%mo;
		y=y*y%mo;
		b/=2;
	}
	return t;
}

void FFT(long long a[],int n,int inv)
{
	if (n==1) return;
	int mid=n/2;
	for (int i=0;i<n;i++)
	{
		rec[i]=(rec[i>>1]>>1)|((i&1)<<(dep-1));
		if (i<rec[i]) swap(a[i],a[rec[i]]);
	}
	for (int len=2;len<=n;len<<=1)
	{
		int mid=len/2;
		for (int j=0;j<mid;j++)
		{
			for (int k=j;k<n;k+=len)
			{
				int p,q;
				q=w[j*(n/len)]*a[k+mid]%mo,p=a[k];
				a[k+mid]=(p-q+mo)%mo;
				a[k]=(p+q+mo)%mo;
			}
		}
	}
}

void solve(int n)
{
	if (n==1)
	{
		d[0]=qsm(br[0],mo-2);
	}
	else
	{
		double op=n;
		int oq=ceil(op/2);
		solve(oq);
		
		up=1;dep=0;
		while (up<=n+jl) up=up*2,dep++;
		w[0]=1;
		rev=qsm(3,(mo-1)/up);
		for (int i=1;i<=up;i++) w[i]=w[i-1]*rev%mo;
		for (int i=0;i<up;i++)
		{
			G[i]=0;
			A[i]=br[i];
			if (i<oq) G[i]=d[i];
		}
		FFT(A,up,1);
		FFT(G,up,1);
		for (int i=0;i<up;i++)
		{
			A[i]=G[i]*(2-A[i]*G[i]%mo+mo)%mo;
		}
		for (int i=0;i<=up/2;i++)
		{
			swap(w[i],w[up-i]);
		}
		FFT(A,up,-1);
		for (int i=0;i<n;i++) d[i]=A[i]*qsm(up,mo-2)%mo;
		int kk=0;
	}
}

int main()
{
	scanf("%d%d",&n,&m);
	for (int i=0;i<=n;i++)
	{
		scanf("%d",&a[i]);
		ar[n-i]=a[i]; 
	}
	for (int i=0;i<=m;i++)
	{
		scanf("%d",&b[i]); 
		br[m-i]=b[i];
	} 
	jl=m;
	solve(n-m+1);
	
	up=1;dep=0;
	while (up<=2*n-m) up=up*2,dep++;
	w[0]=1;
	rev=qsm(3,(mo-1)/up);
	for (int i=1;i<=up;i++) w[i]=w[i-1]*rev%mo;
	FFT(ar,up,1);
	FFT(d,up,1);
	for (int i=0;i<up;i++)
	{
		d[i]=d[i]*ar[i]%mo;
	}
	for (int i=0;i<=up/2;i++)
	{
		swap(w[i],w[up-i]);
	}
	FFT(d,up,-1);
	for (int i=0;i<=n-m;i++) d[i]=d[i]*qsm(up,mo-2)%mo;
	
	int j=0;
	for (int i=n-m;i>=0;i--)
	{
		printf("%d ",d[i]);
		D[j++]=d[i];
	}
	printf("\n");
	
	up=1;dep=0;
	while (up<=n) up=up*2,dep++;
	w[0]=1;
	rev=qsm(3,(mo-1)/up);
	for (int i=1;i<=up;i++) w[i]=w[i-1]*rev%mo;
	FFT(b,up,1);
	FFT(D,up,1);
	for (int i=0;i<up;i++)
	{
		b[i]=b[i]*D[i]%mo;
	}
	for (int i=0;i<=up/2;i++)
	{
		swap(w[i],w[up-i]);
	}
	FFT(b,up,-1);
	for (int i=0;i<n;i++) r[i]=(a[i]-b[i]*qsm(up,mo-2)%mo+mo)%mo;
	for (int i=0;i<m;i++)
	{
		printf("%d ",r[i]);
	}
	printf("\n");
} 

多项式开根

介绍

这玩意是求一个这样的东东：
$$B(x{)}^2\equiv A(x)(mod{x}^n)$$
其中$A(x)$次数为$n-1$，要求一个$B(x)$。注意！由于是模意义下的多项式，所以我们的$B(x)$的次数不一定为$\frac{n-1}{2}$。（废话）
同样考虑用什么递归分治的思想来搞。

前置芝士

这里就需要一点额外的芝士了。
就这题而言，有极大的概率会用到二次剩余。

那么来康康怎么搞。
首先，二次剩余的定义为：
$$\exists 整数x，使得x^2\equiv n(modp)则称n为p的二次剩余$$
再定义一个东东：
$$(\frac{n}{p})=\{\begin{array}{ccc}1(n为p二次剩余)& & \\ -1(n不为p的二次剩余)& & \\ 0(p|n)& & \end{array}$$
这个东东叫做勒让德符号（听说还有个叫做高斯记号的东东）。当p为奇质数的时候（$p=2$的时候每个整数就都是其二次剩余了），还有个性质：
$$(\frac{n}{p})=n^{\frac{p-1}{2}}$$

那么我们再引入两个东东：（现在所有的p都为奇质数）

定理一： 对于$x^2\equiv n(modp)$中有$\frac{p-1}{2}$个n是p的二次剩余。

证明？我不费！

定理二： 首先我们根据定理一，可以考虑random一个数$a\in [0,p-1]$。设$w=a^2-n$。若找到一个$w$满足：$(\frac{w}{p})=-1$或者写成这样：$w^{\frac{p-1}{2}}\equiv -1$，则$(a+\sqrt{w}{)}^{\frac{p+1}{2}}$为一组x值。

证明？这个我费！
$$(a+\sqrt{w}{)}^p=\sum _{i=0}^p{a}^i{\sqrt{w}}^{p-i}\ast C_p^i$$
在模意义下这玩意只有当$i=0$或$i=p$时有值。
那么得到：$$(a+\sqrt{w}{)}^p\equiv a^p+{\sqrt{w}}^p$$
根据费马小定理：$a^{p-1}\equiv 1$
根据条件：$w^{\frac{p+1}{2}}\equiv 1$
所以可以得到：$$a^p+{\sqrt{w}}^p\equiv a-\sqrt{w}$$
那么把原式乘上这玩意即可得到：
$$(a+\sqrt{w}{)}^{p+1}\equiv (a+\sqrt{w})(a-\sqrt{w})\equiv a^2-w\equiv n$$
证毕。

那么得到定理后，我们可以去计算x值了。
一点注意的地方就是在跑快速幂时，我们需要分实数部分和整数部分来计算，因为$\sqrt{w}$我们并不能很好表示，所以要带进去一起算。

代码：待填。QWQ

过程

n=1

既然我们知道了二次剩余，那么我们可以快速得到$n=1$时的情况。
$$B{0}^2\equiv A0(modx)$$
其实有时候并不需要什么二次剩余，因为可能良心出题人给的$A0$是个比较好看的数，比如说1。

n>1

现在才是正轨。
考虑原式：
$$B(x{)}^2\equiv A(x)(mod{x}^n)$$
用同样的套路设$G(x)$表示$A(x)$在模$x^{\lceil \frac{n}{2}\rceil }$意义下的开方。
因此得到：
$$\begin{gathered} G(x{)}^2\equiv A(x)(mod{x}^{\lceil \frac{n}{2}\rceil }) \\ G(x{)}^2-A(x)\equiv 0(mod{x}^{\lceil \frac{n}{2}\rceil }) \\ (G(x{)}^2-A(x){)}^2\equiv 0(mod{x}^n) \\ G(x{)}^4-2G(x{)}^{2}A(x)+A(x{)}^2\equiv 0(mod{x}^n) \\ G(x{)}^4+2G(x{)}^{2}A(x)+A(x{)}^2\equiv 4G(x{)}^{2}A(x)(mod{x}^n) \\ (G(x{)}^2+A(x){)}^2\equiv 4G(x{)}^{2}A(x)(mod{x}^n) \\ \frac{(G(x{)}^2+A(x){)}^2}{4G(x{)}^2}\equiv A(x)(mod{x}^n) \end{gathered}$$

再看到原式，可以得到：
$$B(x)=\frac{G(x{)}^2+A(x)}{2G(x)}$$
愉快递归求解吧~

学习资料：
百度百科
https://blog.csdn.net/stevensonson/article/details/85845334
https://www.cnblogs.com/zhangleo/p/11010302.html
https://www.cnblogs.com/xiefengze1/p/9373272.html

代码

话说递归版的常数是真滴大，洛谷反正我是TLE了QWQ。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
const long long mo=998244353;
const long long inv2=499122177;
const int maxn=300010;

int a[maxn],b[maxn],d[maxn],rev,A[maxn],G[maxn],D[maxn];
int jl,n,rec[maxn],w[maxn],up,dep;

int jia(int a,int b)
{
	return (a+b)%mo;
}

int cheng(int a,int b)
{
	return 1ll*a*b%mo;
}

inline long long qsm(long long a,long long b)
{
	long long t=1;
	long long y=a;
	while (b>0)
	{
		if ((b&1)==1) t=t*y%mo;
		y=y*y%mo;
		b/=2;
	}
	return t;
}

inline void FFT(int a[],int n,int inv)
{
	if (n==1) return;
	int mid=n/2;
	for (int i=0;i<n;i++)
	{
		rec[i]=(rec[i>>1]>>1)|((i&1)<<(dep-1));
		if (i<rec[i]) swap(a[i],a[rec[i]]);
	}
	for (int len=2;len<=n;len<<=1)
	{
		int mid=len/2;
		for (int j=0;j<mid;j++)
		{
			for (int k=j;k<n;k+=len)
			{
				int p,q;
				q=cheng(w[j*(n/len)],a[k+mid]),p=a[k];
				a[k+mid]=jia(p+mo,-q)%mo;
				a[k]=jia(p,q)%mo;
			}
		}
	}
}

inline void solve(int n)
{
	if (n==1)
	{
		d[0]=qsm(b[0],mo-2);
	}
	else
	{
		double op=n;
		int oq=ceil(op/2);
		solve(oq);
		
		up=1;dep=0;
		while (up<=n+jl) up=up*2,dep++;
		w[0]=1;
		rev=qsm(3,(mo-1)/up);
		for (register int i=1;i<=up;i++) w[i]=cheng(w[i-1],rev)%mo;
		for (register int i=0;i<up;i++)
		{
			G[i]=0;A[i]=0;
			if (i<n) A[i]=b[i];
			if (i<oq) G[i]=d[i];
		}
		FFT(A,up,1);
		FFT(G,up,1);
		for (register int i=0;i<up;i++)
		{
			A[i]=cheng(G[i],(2-cheng(A[i],G[i])+mo)%mo);
		}
		for (register int i=0;i<=up/2;i++)
		{
			swap(w[i],w[up-i]);
		}
		FFT(A,up,-1);
		for (register int i=0;i<n;i++) d[i]=cheng(A[i],qsm(up,mo-2));
		for (register int i=n;i<=up;i++) d[i]=0;
		int kk=0;
	}
}

inline void solve1(int n)
{
	if (n==1)
	{
		b[0]=1;
	}
	else
	{
		double op=n;
		int oq=ceil(op/2);
		solve1(oq);
		solve(n);
		
		up=1;dep=0;
		while (up<=n+jl) up=up*2,dep++;
		w[0]=1;
		rev=qsm(3,(mo-1)/up);
		for (register int i=1;i<=up;i++) w[i]=cheng(w[i-1],rev);
		for (register int i=0;i<up;i++)
		{
			A[i]=0;
			if (i<n) A[i]=a[i];
		}
		FFT(A,up,1);FFT(b,up,1);FFT(d,up,1);
		for (register int i=0;i<up;i++) b[i]=cheng(jia(b[i],cheng(A[i],d[i])),inv2);
		for (register int i=0;i<=up/2;i++) swap(w[i],w[up-i]);
		FFT(b,up,-1);
		for (register int i=0;i<n;i++) b[i]=cheng(b[i],qsm(up,mo-2));
		for (register int i=n;i<=up;i++) b[i]=0;
		int kk=0;
	}
}

int main()
{
	scanf("%d",&n);
	jl=n;
	for (int i=0;i<n;i++)
	{
		scanf("%d",&a[i]);
	}
	solve1(n);
	for (int i=0;i<n;i++)
	{
		printf("%d ",b[i]);
	}
	printf("\n");
}