多项式Polynomial

Posted 2021-02-21 chasedeath

多项式Polynomial

前置知识( ext{NTT})

所有操作均在对( ext{998244353})取模下进行

所有的操作均没有经过任何卡常

所有的操作均用$ ext{vector} $来实现，主要是为了理清思路，并且清零问题上会比较容易解决，同时对于每次计算完多项式的长度的要求会显得更加严格

相比经过卡常的代码来说，可读性应该强很多吧

代码总览

ll qpow(ll x,ll k) { // 朴素的快速乘
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res;
}

namespace Polynomial{
	typedef vector <int> Poly;
	#define Mod1(x) ((x>=P)&&(x-=P))
	#define Mod2(x) ((x<0)&&(x+=P))
	void Show(Poly a){ for(int i:a) printf("%d ",i); puts(""); }

	int rev[N];
	int PreMake(int n){
		int R=1,cc=-1;
		while(R<n) R<<=1,cc++;
		rep(i,1,R-1) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
		return R;
	}

	void NTT(int n,Poly &a,int f){ //朴素的NTT模板
		rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
		for(int i=1;i<n;i<<=1) {
			int w=qpow(f==1?3:(P+1)/3,(P-1)/i/2);
			for(int l=0;l<n;l+=i*2) {
				ll e=1;
				for(int j=l;j<l+i;++j,e=e*w%P) {
					int t=a[j+i]*e%P;
					a[j+i]=a[j]-t,Mod2(a[j+i]);
					a[j]+=t,Mod1(a[j]);
				}
			}
		}
		if(f==-1) {
			ll base=qpow(n,P-2);
			rep(i,0,n-1) a[i]=a[i]*base%P;
		}
	}

	Poly operator * (Poly a,Poly b){
		int n=a.size()+b.size()-1,R=PreMake(n);
		a.resize(R),b.resize(R);
		NTT(R,a,1),NTT(R,b,1);
		rep(i,0,R-1) a[i]=1ll*a[i]*b[i]%P;
		NTT(R,a,-1);
		a.resize(n);
		return a;
	}

	Poly operator + (Poly a,Poly b) { rep(i,0,a.size()-1) a[i]+=b[i],Mod1(a[i]); return a; }
	Poly operator - (Poly a,Poly b) { rep(i,0,a.size()-1) a[i]-=b[i],Mod2(a[i]); return a; } 

	Poly Inv(Poly a) { // 多项式求逆
		int n=a.size();
		if(n==1) { Poly tmp; tmp.pb(qpow(a[0],P-2)); return tmp; }
		Poly b=a; b.resize((n+1)/2); b=Inv(b);
		int R=PreMake(n<<1);
		a.resize(R),b.resize(R);
		NTT(R,a,1),NTT(R,b,1);
		rep(i,0,R-1) a[i]=(2-1ll*a[i]*b[i]%P+P)*b[i]%P;
		NTT(R,a,-1);
		a.resize(n);
		return a;
	}

	Poly operator / (vector <int> a,vector <int> b){
		reverse(a.begin(),a.end()),reverse(b.begin(),b.end());
		int n=a.size(),m=b.size();
		b.resize(n-m+1),b=Inv(b);
		a=a*b,a.resize(n-m+1);
		reverse(a.begin(),a.end());
		return a;
	}

	Poly operator % (vector <int> a,vector <int> b) { 
		a=a-a/b*b;
		a.resize(b.size()-1);
		return a;
	}

	Poly Sqrt(vector <int> a){ // 开根号
		int n=a.size();
		if(n==1) { Poly tmp; tmp.pb(1); return tmp; }
		Poly b=a; b.resize((n+1)/2),b=Sqrt(b),b.resize(n);
		Poly c=Inv(b);
		int R=PreMake(n*2);
		a.resize(R),c.resize(R);
		NTT(R,a,1),NTT(R,c,1);
		rep(i,0,R-1) a[i]=1ll*a[i]*c[i]%P;
		NTT(R,a,-1);
		a.resize(n);
		rep(i,0,n-1) a[i]=1ll*(P+1)/2*(a[i]+b[i])%P;
		return a;
	}


	Poly Deri(Poly a){ // 导数
		rep(i,1,a.size()-1) a[i-1]=1ll*i*a[i]%P;
		a.pop_back();
		return a;
	}

	int Mod_Inv[N];
	Poly IDeri(Poly a) { // 原函数
		Mod_Inv[0]=Mod_Inv[1]=1;
		rep(i,2,a.size()+1) Mod_Inv[i]=1ll*(P-P/i)*Mod_Inv[P%i]%P;
		a.pb(0);
		drep(i,a.size()-1,1) a[i]=1ll*a[i-1]*Mod_Inv[i]%P;
		a[0]=0;
		return a;
	}

	Poly ln(Poly a){
		int n=a.size();
		a=Inv(a)*Deri(a);
		a.resize(n-1);
		return IDeri(a);
	}

	Poly Exp(Poly a){ 
		int n=a.size();
		if(n==1) { Poly tmp; tmp.pb(1); return tmp; }
		Poly b=a; b.resize((n+1)/2),b=Exp(b),b.resize(n);
		Poly c=ln(b);
		int R=PreMake(n<<1);
		a.resize(R),b.resize(R),c.resize(R);
		NTT(R,b,1),NTT(R,a,1),NTT(R,c,1);
		rep(i,0,R-1) a[i]=1ll*b[i]*(1-c[i]+a[i]+P)%P;
		NTT(R,a,-1),a.resize(n);
		return a;
	}

	Poly Pow(Poly a,int k){
		a=ln(a);
		rep(i,0,a.size()-1) a[i]=1ll*a[i]*k%P;
		return exp(a);
	}

	#undef Mod1
	#undef Mod2
}

[ ]

[ ]

[ ]

1.多项式求逆

求 (G(x)equiv frac{1}{F(x)} (mod x^n))

考虑递归求解，设已经求出了

[H(x)equiv frac{1}{F(x)},(mod x^{lceil frac{n}{2} ceil}) ]

其中递归边界是(n=1)时,(G(x)=frac{1}{F(0)} (mod P))

则

[H(x)equiv G(x)(mod x^{lceil frac{n}{2} ceil}) ]

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

[(H(x)-G(x))^2equiv 0(mod x^n) ]

注意通过平方可以扩大模数，这很常用

[H(x)^2-2G(x)H(x)+G(x)^2equiv 0(mod x^n) ]

两边乘上(F(x))

[H(x)^2F(x)-2H(x)+G(x)equiv 0(mod x^n) ]

[G(x)equiv 2H(x)-H(x)^2F(x)(mod x^n) ]

[ ]

2.多项式开根号

求(G(x)^2equiv F(x) (mod x^n))

同样的，递归求解，设已经求出了，递归边界是(n=1)时,(G(x)=1)

[H(x)^2equiv F(x) (mod x^{lceil frac{n}{2} ceil}) ]

(H(x)equiv G(x) (mod x^{lceil frac{n}{2} ceil}))

[H(x)^2-2G(x)H(x)+G(x)^2equiv 0(mod x^n) ]

[H(x)^2-2G(x)H(x)+F(x)equiv 0 (mod x^n) ]

[G(x)equiv frac{H(x)^2+F(x)}{2H(x)} (mod x^n) ]

[ ]

3.多项式求(ln)

(G(x)equiv ln F(x) (mod x^n))

(G‘(x)equiv F‘(x)frac{1}{F(x)} (mod x^n))

求出(G‘(x))，然后求原函数即可

[ ]

4.多项式求exp

这个推导其实非常有趣的

不说了，先上一下泰勒公式

[f(x)=sum _{i=0}^{infty}frac{f^{(i)}(x_0)}{i!}(x-x_0)^i ]

其中(f^{(i)}(x))表示(f)的(i)阶导数

即通过(f(x))的导数在(x_0)处的函数值累加得到近似的多项式函数，累加到无穷之后认为相等

题目的问题是求

(G(x)=e^{ F(x)},F(x)=ln G(x))

(ln G(x)-F(x)equiv 0 (mod x^n))

那么我们把元(x)推广为多项式，把题目转化为

对于函数(f(G)=ln G-F)

求出在(mod x^n)意义下的零点

其中(f(x)=ln x-c)

考虑迭代求解，设已经求出(H(x)=e^{F(x)}( mod x^{frac{n}{2}}))

带入(f(G))在(H)上的泰勒展开式

[f(G)=sum_{i=0}^{infty} frac{(f^{(i)}(H))}{i!}(G-H)^i ]

和上面一样可以得到的是

(H(x)equiv G(x) (mod x^{frac{n}{2}}))

((H(x)-G(x))^2 equiv 0 (mod x^n))

所以只用考虑(i=0,1)的两项，后面的((G-H)^i)均为(0)

那么(f(G)=f(H)+f‘(H) cdot (G-H))

求导得到(f‘(x)=frac{1}{x})

所以(f(G)=ln H-F+frac{G-H}{H}equiv 0 (mod x^n))

(G=H(F-ln H+1))

所以就可以好好得递归写代码啦

[ ]

[ ]

5.多项式(k)次幂

(G(x)equiv F(x)^k(mod x^n))

(ln G(x)=k ln F(x) (mod x^n))

求出(ln G(x))之后，(exp)回来即可

很显然这个方法对于开根号也是适用的

[ ]