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

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

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

多项式求逆

介绍

多项式求逆即为求一个多项式的乘法逆元。
它长这样：

\[A(x)*B(x)\equiv 1(mod\ x^n) \]

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

过程

引理一： 考虑一个乘法逆元满足：\(A(x)*B(x)\equiv 1(mod\ x^n)\)，那么我们可以得到另一个柿子：\(A(x)*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}\)意义下的逆元。
由引理得：

\[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})\]

开方（这一步我也不能严谨证明，但是打表发现这是对的）

\[B(x)^2+G(x)^2-2*B(x)*G(x)\equiv 0(mod\ x^n) \]

乘个\(A(x)\)进去，然后由于是模\(x^n\)意义下，所以可以消去一些\(B(x)\):

\[B(x)+A(x)*G(x)^2-2*G(x)\equiv 0(mod\ x^n) \\B(x)\equiv G(x)(2-A(x)*G(x))(mod\ x^n)\]

这时，就可以用\(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)*D(x)+R(x) \]

其中\(A(x)\)和\(B(x)\)已经给出，且\(A(x)\)次数为n，\(B(x)\)次数为m，求\(D(x)\)和\(R(x)\)
也就是除法运算中求商和余数，话说如果上过初中就应该明白长除法。

过程

引理二：
考虑一个反转操作为将一个多项式的系数调换，用数学形式表达即为：

\[Pr(x)=x^n*P(\frac 1 x)=\sum_{i=0}^n a_{n-i}*x^i \]

证明显然。

那么现在就推柿子：

\[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^nA(\frac 1 x)=x^nB(\frac 1 x)D(\frac 1 x)+x^nR(x)\\x^nA(\frac 1 x)=x^mB(\frac 1 x)x^{n-m}D(\frac 1 x)+x^{n-m+1}x^{m-1}R(x) \]

然后再引入引理：

\[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(mod\ p)则称n为p的二次剩余 \]

再定义一个东东：

\[(\frac n p)=\left\{\begin{matrix}1(n为p二次剩余) & & \\ -1(n不为p的二次剩余) & & \\0(p|n) \end{matrix}\right.\]

这个东东叫做勒让德符号（听说还有个叫做高斯记号的东东）。当p为奇质数的时候（\(p=2\)的时候每个整数就都是其二次剩余了），还有个性质：

\[(\frac n p)=n^{\frac {p-1} 2} \]

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

定理一： 对于\(x^2\equiv n(mod\ p)\)中有\(\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}^pa^i\sqrt w^{p-i}*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\)时的情况。

\[B0^2\equiv A0(mod\ x) \]

其实有时候并不需要什么二次剩余，因为可能良心出题人给的\(A0\)是个比较好看的数，比如说1。

n>1

现在才是正轨。
考虑原式：

\[B(x)^2\equiv A(x)(mod\ x^n) \]

用同样的套路设\(G(x)\)表示\(A(x)\)在模\(x^{\lceil \frac n 2 \rceil}\)意义下的开方。
因此得到：

\[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)^2A(x)+A(x)^2\equiv 0(mod\ x^n)\\G(x)^4+2G(x)^2A(x)+A(x)^2\equiv 4G(x)^2A(x)(mod\ x^n)\\(G(x)^2+A(x))^2\equiv4G(x)^2A(x)(mod\ x^n)\\ \frac{(G(x)^2+A(x))^2}{4G(x)^2}\equiv A(x)(mod\ x^n)\]

再看到原式，可以得到：

\[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");
}