• 多项式全家桶——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 ceil})意义下的逆元。
    由引理得:

    [B(x)equiv G(x)(mod x^{lceil frac n 2 ceil})\B(x)-G(x)equiv 0(mod x^{lceil frac n 2 ceil})]

    开方(这一步我也不能严谨证明,但是打表发现这是对的)

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

    多项式除法

    介绍

    这玩意是求一个长成这样的柿子的:

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

    多项式开根

    介绍

    这玩意是求一个这样的东东:

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

    其中(A(x))次数为(n-1),要求一个(B(x))。注意!由于是模意义下的多项式,所以我们的(B(x))的次数不一定为(frac{n-1}2)。(废话)
    同样考虑用什么递归分治的思想来搞。

    前置芝士

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

    那么来康康怎么搞。
    首先,二次剩余的定义为:

    [exists 整数x,使得x^2equiv n(mod p)则称n为p的二次剩余 ]

    再定义一个东东:

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

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

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

    那么我们再引入两个东东:(现在所有的p都为奇质数)

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

    证明?我不费!

    定理二: 首先我们根据定理一,可以考虑random一个数(ain[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^isqrt w^{p-i}*C_p^i ]

    在模意义下这玩意只有当(i=0)(i=p)时有值。
    那么得到:

    [(a+sqrt w)^pequiv a^p+sqrt w^p ]

    根据费马小定理:(a^{p-1}equiv 1)
    根据条件:(w^{frac{p+1} 2}equiv 1)
    所以可以得到:

    [a^p+sqrt w^pequiv a-sqrt w ]

    那么把原式乘上这玩意即可得到:

    [(a+sqrt w)^{p+1}equiv(a+sqrt w)(a-sqrt w)equiv a^2-wequiv n ]

    证毕。

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

    代码:待填。QWQ

    过程

    n=1

    既然我们知道了二次剩余,那么我们可以快速得到(n=1)时的情况。

    [B0^2equiv A0(mod x) ]

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

    n>1

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

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

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

    [G(x)^2equiv A(x)(mod x^{lceil frac n 2 ceil})\G(x)^2-A(x)equiv 0(mod x^{lceil frac n 2 ceil})\(G(x)^2-A(x))^2equiv 0(mod x^n)\G(x)^4-2G(x)^2A(x)+A(x)^2equiv 0(mod x^n)\G(x)^4+2G(x)^2A(x)+A(x)^2equiv 4G(x)^2A(x)(mod x^n)\(G(x)^2+A(x))^2equiv4G(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("
    ");
    }
    
  • 相关阅读:
    css去掉iPhone、iPad默认按钮样式
    STL~Deque简介
    Centos 7 ssh登录速度慢
    C++ delete 两次
    编译gdb 报错 No module named gdb.frames
    gdb 脚本
    转载: CentOS/Linux 解决 SSH 连接慢
    百度经验:Win10查看已存储WiFi密码的两种方法
    git 操作
    Avoiding memory leaks in POSIX thread programming, 多线程避免内存泄漏
  • 原文地址:https://www.cnblogs.com/RainbowCrown/p/13454910.html
Copyright © 2020-2023  润新知