多项式总结及模板

89 篇文章 0 订阅
72 篇文章 0 订阅

1.求逆

设要对多项式A求逆,逆为B:

求在模ceil(n/2)的逆,将原式与之相减后平方,发现可得在模n意义下也是0了,然后两边同乘A,移项即可倍增了(注意此时A也是在模当前长度意义下的)。

最后倍增形式为B'=2*B-B*B*A  (b’等于2b减b方a),边界为常数取逆元,也可发现多项式有无逆取决于常数有无逆元

2.开根

与求逆类似,先求ceil(n/2),然后将A移到左边使等式右边变成0,平方,加上4b^2a,再把右边b^2除到左边就发现ok了。

最终倍增形式:B'=(A+B*B)/2B(B'等于A加B方后除以2B),除法就上求逆,边界是常数开根,暴力枚举(并不会二次剩余...)?模板题保证是1...

(开根里面套求逆,虽说复杂度仍然是nlog,但已经比大部分两个log慢了。。。)

 3.Ln

两边求导得到所求的导数B’=A'/A,再积分即可。

4.Exp

B=e^(A)

根据牛顿迭代的那套理论,设一个函数C,使C(B)=0,那么C(B)=LnB-A

假设已知在mod x^n的B0,要求在mod x^2n的B1

有C(B)=C(B0)/0!+C'(B0)/1!*(B-B0)+...

发现...可以省略,因为在mod x^n下B是唯一确定的,而mod x^2n下一定也满足mod x^n下,所以在mod x^2n时后面低的n位与mod x^n相等,当求的导数大于等于二阶发现(B-B0)那一项的幂次大于等于2,也就是最小的x项幂次大于等于2n,省略即可。

最终推得

B1=B0-C(B0)/C'(B0)

而C求导时是在把B0这个多项式看成变量求导(不要把每一项的x看成变量),A是已知多项式要看成一个常数

故B1=B0(1-Ln(B0)+A)

倍增即可。(注意能Exp一定要满足0次项为0,不然无法赋初值)

 

先贴4个的代码:

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int mod=998244353;
const int N=1e5+100;

void Ad(int &x,int y)
{if((x+=y)>=mod)x-=mod;}
void Mul(int &x,int y)
{x=1LL*x*y%mod;}
int qpow(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)res=1LL*res*x%mod;
        x=1LL*x*x%mod,y>>=1;
    }
    return res;
}

namespace Poly{

vector<int>A,B,C;
int wh[N*3],cc,len,n,inv[N*3];
void Pre_Ntt(int l)
{
    cc=0,len=1;
    while(len<=l)++cc,len<<=1;
}
void Cal_wh()
{for(int i=1;i<len;i++)wh[i]=(wh[i>>1]>>1)|((i&1)<<(cc-1));}
void Ntt(vector<int>&a,bool inv)
{
    for(int i=0;i<len;i++)
        if(i<wh[i])swap(a[i],a[wh[i]]);
    int tp,mo,ha;
    for(int l=2,md;l<=len;l<<=1)
    {
        md=l>>1,tp=qpow(3,(mod-1)/l);
        for(int i=0;i<len;i+=l)
        {
            mo=1;
            for(int j=0;j<md;j++,mo=1LL*mo*tp%mod)
            {
                ha=1LL*mo*a[i+j+md]%mod;
                a[i+j+md]=(a[i+j]-ha+mod)%mod;
                a[i+j]=(a[i+j]+ha)%mod;
            }
        }
    }
    if(inv)
    {
        tp=qpow(len,mod-2);
        for(int i=1;i<len/2;i++)swap(a[i],a[len-i]);
        for(int i=0;i<len;i++)Mul(a[i],tp);
    }
}

void Get_Inv(vector<int>&b,int n)
{
    if(n==1){b[0]=qpow(A[0],mod-2);return;}
    int md=(n+1)/2;
    vector<int>a,c;
    Pre_Ntt(n+md+md),a.resize(len),c.resize(len);
    Get_Inv(c,md);
    Pre_Ntt(n+md+md),Cal_wh();
    for(int i=0;i<n;i++)a[i]=A[i];
    Ntt(a,0),Ntt(c,0);
    for(int i=0;i<len;i++)a[i]=(2*c[i]-1LL*c[i]*c[i]%mod*a[i]%mod+mod)%mod;
    Ntt(a,1);
    for(int i=0;i<n;i++)b[i]=a[i];
}

void Inv(vector<int> &a,vector<int> &res,int n)
{
    A=a,res.clear(),res.resize(n);
    Get_Inv(res,n);
}

void Get_Sqrt(vector<int>&b,int n)
{
    if(n==1){b[0]=1;return;}
    int md=(n+1)/2;
    vector<int>a,c;
    Pre_Ntt(n+n),a.resize(len),c.resize(len);
    Get_Sqrt(c,md);
    for(int i=0;i<n;i++)b[i]=2*c[i]%mod;
    Inv(b,b,n);
    Pre_Ntt(n+n),Cal_wh();
    Ntt(c,0);
    for(int i=0;i<len;i++)c[i]=1LL*c[i]*c[i]%mod;
    Ntt(c,1);
    for(int i=0;i<n;i++)a[i]=(B[i]+c[i])%mod;
    Ntt(a,0),Ntt(b,0);
    for(int i=0;i<len;i++)b[i]=1LL*b[i]*a[i]%mod;
    Ntt(b,1);
}

void Sqrt(vector<int> &a,vector<int> &res,int n)
{
    Pre_Ntt(n+n);
    B=a,res.clear(),res.resize(len);
    Get_Sqrt(res,n);
}

void Dao(vector<int> &a,vector<int> &res,int n)
{
	res.clear(),res.resize(n);
	for(int i=1;i<n;i++)
		res[i-1]=1LL*a[i]*i%mod;
	res[n-1]=0;
}

void Jf(vector<int> &a,vector<int> &res,int n)
{
	res.clear(),res.resize(n+1),inv[1]=1;
	for(int i=2;i<=n;i++)
		inv[i]=1LL*(mod-mod/i)*inv[mod%i]%mod;
	for(int i=n-1;i>=0;i--)
		res[i+1]=1LL*a[i]*inv[i+1]%mod;
	res[0]=0;
}

void Ln(vector<int>&a,vector<int> &res,int n)
{
	static vector<int>mo,ha;
	Dao(a,mo,n),Inv(a,ha,n);
	Pre_Ntt(n+n),Cal_wh();
	mo.resize(len),Ntt(mo,0);
	ha.resize(len),Ntt(ha,0);
	for(int i=0;i<len;i++)mo[i]=1LL*mo[i]*ha[i]%mod;
	Ntt(mo,1),Jf(mo,res,n);
}

void Get_Exp(vector<int>&b1,int n)
{
	if(n==1){b1[0]=1;return;}
	int md=(n+1)/2;
	vector<int>b0;
	b0.resize(n),b1.resize(n);
	Get_Exp(b0,md),Ln(b0,b1,n);
	for(int i=0;i<n;i++)b1[i]=((i==0)-b1[i]+C[i]+mod)%mod;
	Pre_Ntt(n+n),b0.resize(len),b1.resize(len);
	Cal_wh(),Ntt(b0,0),Ntt(b1,0);
	for(int i=0;i<len;i++)
		b1[i]=1LL*b0[i]*b1[i]%mod;
	Ntt(b1,1);
}

void Exp(vector<int>&a,vector<int> &res,int n)
{C=a,res.resize(n),Get_Exp(res,n);}

}

int main()
{
    int n;
    scanf("%d",&n);
    vector<int>a,b;
    a.resize(n);
    for(int i=0;i<n;i++)
        scanf("%d",&a[i]);
    Poly::Exp(a,b,n);
    for(int i=0;i<n;i++)
        printf("%d ",b[i]);
    puts("");
}

刷到一道题要多项式取模/除法,补上:

https://www.luogu.org/problemnew/show/P4512

A(x)=B(x)D(x)+R(x)

翻转,即最高n次的F(x) -> x^n*F(1/x)

整个等式操作一下:

x^n*A(1/x)=x^n(B(1/x)D(1/x)+R(1/x))

拆开:

x^nA(1/x)=x^mB(1/x)x^(n-m)D(1/x)+x^(m-1)R(1/x)x^(n-m+1)

换个元:

A'=B'D' mod(x n-m+1)

即可求D反转后的东西,求完转回来即可,R也随之解决

代码:

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const ll mod=998244353;
const int N=1e5+100;

ll qpow(ll x,ll y=mod-2)
{
	ll res=1;
	while(y)
	{
		if(y&1)res=res*x%mod;
		x=x*x%mod,y>>=1;
	}
	return res;
}
inline ll Add(ll x,ll y)
{return x+y>=mod?x+y-mod:x+y;}
inline ll Sub(ll x,ll y)
{return x-y<0?x-y+mod:x-y;}
inline ll Mul(ll x,ll y)
{return x*y%mod;}

namespace Poly{
	int wh[N<<2],len,cc;
	void pre_ntt(int nn)
	{
		len=1,cc=0;
		while(len<nn)len<<=1,++cc;
		for(int i=1;i<len;i++)
			wh[i]=(wh[i>>1]>>1)|((i&1)<<(cc-1));
	}
	void ntt(vector<ll> &a,bool inv)
	{
		for(int i=0;i<len;i++)
			if(i<wh[i])swap(a[i],a[wh[i]]);
		ll tp,mo,ha;
		for(int l=2,md;l<=len;l<<=1)
		{
			md=l>>1,tp=qpow(3,(mod-1)/l);
			for(int i=0;i<len;i+=l)
			{
				mo=1;
				for(int j=0;j<md;j++,mo=Mul(mo,tp))
				{
					ha=Mul(a[i+j+md],mo);
					a[i+j+md]=Sub(a[i+j],ha);
					a[i+j]=Add(a[i+j],ha);
				}
			}
		}
		if(inv)
		{
			tp=qpow(len);
			for(int i=1;i<len/2;i++)swap(a[i],a[len-i]);
			for(int i=0;i<len;i++)a[i]=Mul(a[i],tp);
		}
	}
	vector<ll> Poly_Mul(vector<ll> x,vector<ll> y)
	{
		pre_ntt(x.size()+y.size());
		x.resize(len),y.resize(len);
		ntt(x,0),ntt(y,0);
		for(int i=0;i<len;i++)
			x[i]=Mul(x[i],y[i]);
		ntt(x,1);
		return x;
	}
	vector<ll> Inv(vector<ll> x)
	{
		int lim=1;vector<ll>A,B;
		while(lim<x.size())lim<<=1;
		for(int i=1;i<=lim;i<<=1)
		{
			if(i==1){B.resize(1),B[0]=qpow(x[0]);continue;}
			A.resize(i);
			for(int j=0;j<i&&j<x.size();j++)A[j]=x[j];
			pre_ntt(i+B.size());
			A.resize(len),B.resize(len);
			ntt(A,0),ntt(B,0);
			for(int i=0;i<len;i++)
				B[i]=(2*B[i]-A[i]*B[i]%mod*B[i]%mod+mod)%mod;
			ntt(B,1),B.resize(i);
		}
		B.resize(x.size());
		return B;
	}
	vector<ll> Div(vector<ll> A,vector<ll> B)
	{
		int n=A.size()-1,m=B.size()-1;
		reverse(A.begin(),A.end());
		reverse(B.begin(),B.end());
		A.resize(n-m+1),B.resize(n-m+1);
		B=Inv(B),A=Poly_Mul(A,B),A.resize(n-m+1);
		reverse(A.begin(),A.end());
		return A;
	}
	vector<ll> Poly_Sub(vector<ll> x,vector<ll> y)
	{
		if(x.size()<y.size())x.resize(y.size());
		for(int i=0;i<y.size();i++)
			x[i]=Sub(x[i],y[i]);
		return x;
	}
}

int main()
{
	vector<ll>x,y,res;
	int n,m;
	scanf("%d%d",&n,&m),++n,++m;
	x.resize(n),y.resize(m);
	for(int i=0;i<n;i++)
		scanf("%lld",&x[i]);
	for(int i=0;i<m;i++)
		scanf("%lld",&y[i]);
	res=Poly::Div(x,y);
	for(int i=0;i<res.size();i++)
		printf("%lld ",res[i]);
	puts("");
	res=Poly::Poly_Sub(x,Poly::Poly_Mul(y,res));
	for(int i=0;i<m-1;i++)
		printf("%lld ",res[i]);
	puts("");
}

 

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值