多项式 求逆/除法/取模/开根/ln/exp/求幂 (待添加多点求值)

43 篇文章 0 订阅
39 篇文章 6 订阅


模板直接到最底下

多项式求逆

给定 A ( x ) A(x) A(x) B ( x ) = A − 1 ( x ) ( m o d x n ) B(x)=A^{-1}(x)(mod x^n) B(x)=A1(x)(modxn)
使用倍增构造。
若n=1,则直接求0次项的乘法逆元。
若n>1
B ( x ) B(x) B(x)满足
A ( x ) B ( x ) ≡ 1 m o d    x n / 2 A(x)B(x)\equiv 1 \mod x^{n/2} A(x)B(x)1modxn/2

A ( x ) B ( x ) − 1 ≡ 0 m o d    x n / 2 A(x)B(x)-1\equiv 0 \mod x^{n/2} A(x)B(x)10modxn/2

( A ( x ) B ( x ) − 1 ) 2 ≡ 0 m o d    x n (A(x)B(x)-1)^2\equiv 0 \mod x^{n} (A(x)B(x)1)20modxn

A 2 ( x ) B 2 ( x ) − 2 A ( x ) B ( x ) + 1 ≡ 0 m o d    x n A^2(x)B^2(x)-2A(x)B(x)+1\equiv 0 \mod x^{n} A2(x)B2(x)2A(x)B(x)+10modxn

A ( x ) ( 2 B ( x ) − A ( x ) B 2 ( x ) ) ≡ 1 m o d    x n A(x)(2B(x)-A(x)B^2(x))\equiv 1 \mod x^{n} A(x)(2B(x)A(x)B2(x))1modxn

2 B ( x ) − A ( x ) B 2 ( x ) ≡ A − 1 ( x ) m o d    x n 2B(x)-A(x)B^2(x)\equiv A^{-1}(x) \mod x^{n} 2B(x)A(x)B2(x)A1(x)modxn

递归即可,复杂度 O ( n l o g ( n ) ) O(nlog(n)) O(nlog(n))

多项式除法

已知 A ( x ) , B ( x ) A(x),B(x) A(x),B(x),它们的次数分别为 n , m ( n > m ) n,m(n>m) n,m(n>m)

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

C ( x ) C(x) C(x)是除法, D ( x ) D(x) D(x)是取模,其中 C ( x ) C(x) C(x)的系数不大于 n − m n-m nm D ( x ) D(x) D(x)的系数小于 m m m

首先可以发现 A ′ ( x ) = x n A ( 1 x ) A'(x)=x^nA({1\over x}) A(x)=xnA(x1)就是将A的系数翻转过来。

同乘 x n x^n xn
那么
x n A ( 1 x ) = x m B ( 1 x ) x n − m C ( 1 x ) + x n − m + 1 x m − 1 D ( 1 x ) x^nA({1\over x})=x^mB({1\over x})x^{n-m}C({1\over x})+x^{n-m+1}x^{m-1}D({1\over x}) xnA(x1)=xmB(x1)xnmC(x1)+xnm+1xm1D(x1)

将这个东西 m o d    x n − m + 1 \mod x^{n-m+1} modxnm+1,可以发现对 C ( x ) C(x) C(x)的系数不会有影响。
于是
A ′ ( x ) ≡ B ′ ( x ) C ′ ( x ) m o d    x n − m + 1 A'(x)\equiv B'(x)C'(x) \mod x^{n-m+1} A(x)B(x)C(x)modxnm+1

C ′ ( x ) ≡ A ′ ( x ) B ′ − 1 ( x ) m o d    x n − m + 1 C'(x)\equiv A'(x)B'^{-1}(x) \mod x^{n-m+1} C(x)A(x)B1(x)modxnm+1

于是多项式求逆,然后乘起来就没了。

多项式取模

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

先除,再乘,再减,就没了。

多项式开方

给定 A ( x ) A(x) A(x) B 2 ( x ) = A ( x ) ( m o d x n ) B^2(x)=A(x)(mod x^n) B2(x)=A(x)(modxn)
使用倍增构造。
当n=1时,相当于求0次项的二次剩余。
(然而我并不会二次剩余,而且很多题里的常数项是真的常数,所以这里忽略二次剩余)

补一下二次剩余这个坑吧
我这辈子应该都学不会那个C开头的算法了。
但是可以BSGS对原根求出离散对数,把指数除个2,再幂回去。
对于二次剩余这个做法是正确的,如果指数是奇数不能除2,则这个数是没有二次剩余的。

(好像用后面的多项式求幂也能求这个)
当n>1时,
B ( x ) , G ( x ) B(x),G(x) B(x),G(x)满足
B 2 ( x ) ≡ A ( x ) m o d    x n / 2 , G 2 ( x ) ≡ A ( x ) m o d    x n B^2(x)\equiv A(x) \mod x^{n/2},G^2(x)\equiv A(x) \mod x^{n} B2(x)A(x)modxn/2,G2(x)A(x)modxn

G 2 ( x ) − B 2 ( x ) ≡ 0 m o d    x n / 2 G^2(x)-B^2(x)\equiv 0 \mod x^{n/2} G2(x)B2(x)0modxn/2

两边平方
G 4 ( x ) + B 4 ( x ) − 2 G 2 ( x ) B 2 ( x ) ≡ 0 m o d    x n G^4(x)+B^4(x)-2G^2(x)B^2(x)\equiv 0 \mod x^{n} G4(x)+B4(x)2G2(x)B2(x)0modxn

( G 2 ( x ) + B 2 ( x ) ) 2 ≡ 4 G 2 ( x ) B 2 ( x ) m o d    x n / 2 (G^2(x)+B^2(x))^2\equiv 4G^2(x)B^2(x) \mod x^{n/2} (G2(x)+B2(x))24G2(x)B2(x)modxn/2

G 2 ( x ) + B 2 ( x ) ≡ 2 G ( x ) B ( x ) m o d    x n / 2 G^2(x)+B^2(x)\equiv 2G(x)B(x) \mod x^{n/2} G2(x)+B2(x)2G(x)B(x)modxn/2

因为 G 2 ( x ) ≡ A ( x ) m o d    x n G^2(x)\equiv A(x) \mod x^{n} G2(x)A(x)modxn
A ( x ) + B 2 ( x ) ≡ 2 G ( x ) B ( x ) m o d    x n / 2 A(x)+B^2(x)\equiv 2G(x)B(x) \mod x^{n/2} A(x)+B2(x)2G(x)B(x)modxn/2

A ( x ) 2 B ( x ) + B ( x ) 2 ≡ G ( x ) m o d    x n / 2 {A(x)\over 2B(x)}+{B(x)\over 2}\equiv G(x) \mod x^{n/2} 2B(x)A(x)+2B(x)G(x)modxn/2

求逆然后乘一乘,加一加就行了。

多项式求ln

l n ( A ( x ) ) = ∫ [ l n ( A ( x ) ) ] ′ d x ln(A(x))=\int [ln(A(x))]'dx ln(A(x))=[ln(A(x))]dx
= ∫ A ′ ( x ) A ( x ) d x =\int {A'(x)\over A(x)}dx =A(x)A(x)dx
求导,求逆,然后积分即可。

多项式exp

牛顿迭代

有函数 g ( ) g() g(),解多项式 f f f的前n项使 g ( f ) = 0 g(f)=0 g(f)=0
f f f有2n项,设前n项使 f 0 f_0 f0
g ( f ) g(f) g(f) f 0 f_0 f0处泰勒展开

g ( f ) ≡ g ( f 0 ) + g ′ ( f 0 ) ( f − f 0 ) + g ′ ′ ( f 0 ) ( f − f 0 ) 2 2 ! + . . . m o d    x 2 n g(f)\equiv g(f_0)+g'(f_0)(f-f_0)+{g''(f0)(f-f_0)^2\over 2!}+... \mod x^{2n} g(f)g(f0)+g(f0)(ff0)+2!g(f0)(ff0)2+...modx2n

注意到从第三项开始,在mod 2n意义下=0
所以
g ( f ) ≡ g ( f 0 ) + g ′ ( f 0 ) ( f − f 0 ) m o d    x 2 n g(f)\equiv g(f_0)+g'(f_0)(f-f_0) \mod x^{2n} g(f)g(f0)+g(f0)(ff0)modx2n

∵ g ( f ) = 0 \because g(f)=0 g(f)=0

∴ g ( f 0 ) + g ′ ( f 0 ) ( f − f 0 ) ≡ 0 m o d    x 2 n \therefore g(f_0)+g'(f_0)(f-f_0)\equiv 0 \mod x^{2n} g(f0)+g(f0)(ff0)0modx2n

f ≡ f 0 − g ( f 0 ) g ′ ( f 0 ) m o d    x 2 n f\equiv f_0-{g(f_0)\over g'(f_0)} \mod x^{2n} ff0g(f0)g(f0)modx2n

求exp

B ( x ) ≡ e x p ( A ( x ) ) m o d    x n B(x)\equiv exp(A(x)) \mod x^n B(x)exp(A(x))modxn
g ( x ) = l n ( x ) − A g(x)=ln(x)-A g(x)=ln(x)A,那么目标是 g ( B ) = 0 g(B)=0 g(B)=0
B ≡ B 0 − g ( B 0 ) g ′ ( B 0 ) m o d    x n B\equiv B_0-{g(B_0)\over g'(B_0)}\mod x^n BB0g(B0)g(B0)modxn

∵ g ( B 0 ) ′ = 1 B 0 \because g(B_0)'={1\over B_0} g(B0)=B01

∴ B ≡ B 0 − B 0 ( l n ( B 0 ) − A ( x ) ) m o d    x n \therefore B\equiv B_0-B_0(ln(B_0)-A(x))\mod x^n BB0B0(ln(B0)A(x))modxn

B ≡ B 0 ( 1 + A ( x ) − l n ( B 0 ) ) m o d    x n B\equiv B_0(1+A(x)-ln(B_0))\mod x^n BB0(1+A(x)ln(B0))modxn
当n=1时直接求(程序中默认常数项=0)
当n>1时倍增

多项式求幂

先ln,再乘上幂数,再exp。

代码

这个代码对应的是bzoj3625 CF438E

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
#define N 3201000
#define ll long long
#define mo 998244353
using namespace std;
int n,m;
ll len,a[N],b[N],c[N];
ll mi(ll a,ll b)
{
	a%=mo;ll jy=1;
	for(;b;b/=2,a=a*a%mo) if(b%2==1) jy=jy*a%mo;
	return jy;
}
struct polynomial{
	ll D[N],W[N],A[N],B[N],C[N],E[N],F[N];
	void DFT(ll *a,int sig,int deg)
	{
		int len,lg;
		for(len=1,lg=0;len<=deg;len*=2,lg++);
		W[0]=1;W[1]=mi(3,(mo-1)/len);fo(i,2,len) W[i]=W[i-1]*W[1]%mo;
		fo(i,0,len-1)
		{
			int p=0;
			for(int j=0,tp=i;j<lg;j++,tp/=2) p=(p<<1)+(tp%2);
			D[p]=a[i];
		}
		for(int m=2;m<=len;m*=2)
		{
			int half=m/2,tmp=len/m;
			for(int j=0;j<len;j+=m)
			{
				fo(i,j,j+half-1)
				{
					ll u=D[i],v=D[i+half]*((sig==1)?W[(i-j)*tmp]:W[len-(i-j)*tmp])%mo;
					D[i]=(u+v)%mo;
					D[i+half]=(u-v+mo)%mo;
				}
			}
		}
		if(sig==-1)
		{
			ll inv=mi(len,mo-2);
			fo(i,0,len-1) D[i]=D[i]*inv%mo;
		}
		fo(i,0,len-1) a[i]=D[i];
	}
	void mul(ll *a,ll *b,ll *c,int len) //C
	{
		fo(i,0,len-1) A[i]=a[i],B[i]=b[i];
		DFT(A,1,len);DFT(B,1,len);
		fo(i,0,len+len) A[i]=A[i]*B[i]%mo;
		DFT(a,-1,len);
		fo(i,0,len+len) c[i]=A[i];
	}
	void inv(ll *a,ll *b,int deg) ///B
	{
		if(deg==1)
		{
			b[0]=mi(a[0],mo-2);
			b[1]=0;
			return;
		} 
		inv(a,b,deg/2);
		fo(i,0,deg-1) A[i]=a[i],B[i]=b[i];
		fo(i,deg,deg+deg-1) A[i]=B[i]=0;
		DFT(A,1,deg);DFT(B,1,deg);
		fo(i,0,deg+deg-1) A[i]=(A[i]*B[i]%mo*B[i]%mo)%mo;
		DFT(A,-1,deg);
		fo(i,0,deg-1) b[i]=(b[i]*2ll-A[i]+mo)%mo;
		fo(i,deg,deg+deg) b[i]=0;
	}
	void divide(ll *a,int n,ll *b,int m,ll *c) //C
	{
		int len=1;
		for(;len<n-m+1;len*=2);
		fo(i,0,m-1) c[i]=b[m-i-1],C[i]=0;
		fo(i,min(m,n-m+1),len+len) c[i]=C[i]=0;
		inv(c,C,len);
		fo(i,0,n-1) A[i]=a[n-i-1];
		fo(i,n-m+1,len+len) A[i]=C[i]=0;
		DFT(A,1,len);DFT(C,1,len);
		fo(i,0,len+len) A[i]=A[i]*C[i]%mo;
		DFT(A,-1,len);
		fo(i,0,n-m) c[i]=A[n-m-i];	
	}
	void mod(ll *a,int n,ll *b,int m,ll *c) //C
	{
		int len=1;
		for(;len<=n;len*=2);
		divide(a,n,b,m,c);
		fo(i,n-m+1,len+len) c[i]=0;
		fo(i,0,m-1) A[i]=b[i];
		fo(i,m,len+len) A[i]=0;
		DFT(A,1,len);DFT(c,1,len);
		fo(i,0,len+len) A[i]=A[i]*c[i]%mo;
		DFT(A,-1,len);
		fo(i,0,m-2) c[i]=(a[i]-A[i]+mo)%mo;
	}
	void sqrt(ll *a,ll *b,int deg)  //C
	{
		if(deg==1)
		{
			b[0]=1;   //默认常数为1(不用二次剩余)
			return;
		}
		sqrt(a,b,deg/2);
		inv(b,C,deg);
		fo(i,0,deg-1) A[i]=a[i];
		fo(i,deg,deg+deg) A[i]=0,C[i]=0;
		DFT(A,1,deg);DFT(C,1,deg);
		fo(i,0,deg+deg-1) A[i]=A[i]*C[i]%mo;
		DFT(A,-1,deg);
		ll ny2=mi(2,mo-2);
		fo(i,0,deg-1) b[i]=(b[i]*ny2%mo+A[i]*ny2%mo)%mo;
	}
	void ln(ll *a,ll *b,int deg)  //C
	{
		inv(a,C,deg);
		fo(i,0,deg-1) b[i]=1ll*a[i+1]*(i+1)%mo;
		fo(i,deg,deg+deg) C[i]=b[i]=0;
		DFT(b,1,deg);DFT(C,1,deg);
		fo(i,0,deg+deg-1) b[i]=b[i]*C[i]%mo;
		DFT(b,-1,deg);
		fd(i,deg,1) b[i]=b[i-1]*mi(i,mo-2)%mo;
		fo(i,deg,deg+deg) b[i]=0;
		b[0]=0;
	}
	void exp(ll *a,ll *b,int deg)  //E
	{
		if(deg==1)
		{
			b[0]=1;  //默认常数为0
			b[1]=0;
			return;
		}
		exp(a,b,deg/2);
		ln(b,E,deg);
		fo(i,0,deg-1) E[i]=(a[i]-E[i]+mo)%mo;
		E[0]=(1+E[0]+mo)%mo;
		fo(i,deg,deg+deg) b[i]=E[i]=0;
		DFT(b,1,deg);DFT(E,1,deg);
		fo(i,0,deg+deg-1) b[i]=b[i]*E[i]%mo;
		DFT(b,-1,deg);
		fo(i,deg,deg+deg) b[i]=0;
	}
	void pmi(ll *a,ll *b,ll sig,ll deg) //F
	{
		ln(a,F,deg);
		fo(i,0,deg-1) F[i]=F[i]*sig%mo;
		exp(F,b,deg);
	}
}poly;
int main()
{
	int mx=0;
	scanf("%d%d",&n,&m);
	fo(i,1,n)
	{
		int x;scanf("%d",&x);
		a[x]=mo-4;
		mx+=x;
	}
	for(len=1;len<=m;len*=2);
	a[0]++;
	poly.pmi(a,b,mi(2,mo-2),len);
	b[0]++;
	poly.inv(b,a,len);
	fo(i,1,m) printf("%lld\n",2ll*a[i]%mo);
}
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值