[营业日志]萌新的多项式入门

目录

• 多项式工业入门
  • 一些定义
  • 基础快速变幻
  • 多项式加、减、乘
  • 多项式求逆
  • 求导数、原函数
  • 多项式$ln$
  • 多项式$exp$
  • 多项式开根
    • $a_0=1$
    • $a_0≠1$留坑
  • 多项式快速幂
    • $a_0=1$
    • $a_0≠1$
  • 多项式三角函数
  • 多项式反三角函数
• 多项式黑科技
  • 任意模数NTT(MTT)

同步更新于caijiのBlog

菜鸡不会多项式石锤了

多项式工业入门

一些定义

template<const int mod>
struct modint{
    int x;
    modint<mod>(int o=0){x=o;}
    modint<mod> &operator = (int o){return x=o,*this;}
    modint<mod> &operator +=(modint<mod> o){return x=x+o.x>=mod?x+o.x-mod:x+o.x,*this;}
    modint<mod> &operator -=(modint<mod> o){return x=x-o.x<0?x-o.x+mod:x-o.x,*this;}
    modint<mod> &operator *=(modint<mod> o){return x=1ll*x*o.x%mod,*this;}
    modint<mod> &operator ^=(int b){
        modint<mod> a=*this,c=1;
        for(;b;b>>=1,a*=a)if(b&1)c*=a;
        return x=c.x,*this;
    }
    modint<mod> &operator /=(modint<mod> o){return *this *=o^=mod-2;}
    modint<mod> &operator +=(int o){return x=x+o>=mod?x+o-mod:x+o,*this;}
    modint<mod> &operator -=(int o){return x=x-o<0?x-o+mod:x-o,*this;}
    modint<mod> &operator *=(int o){return x=1ll*x*o%mod,*this;}
    modint<mod> &operator /=(int o){return *this *= ((modint<mod>(o))^=mod-2);}
	template<class I>friend modint<mod> operator +(modint<mod> a,I b){return a+=b;}
    template<class I>friend modint<mod> operator -(modint<mod> a,I b){return a-=b;}
    template<class I>friend modint<mod> operator *(modint<mod> a,I b){return a*=b;}
    template<class I>friend modint<mod> operator /(modint<mod> a,I b){return a/=b;}
    friend modint<mod> operator ^(modint<mod> a,int b){return a^=b;}
    friend bool operator ==(modint<mod> a,int b){return a.x==b;}
    friend bool operator !=(modint<mod> a,int b){return a.x!=b;}
    bool operator ! () {return !x;}
    modint<mod> operator - () {return x?mod-x:0;}
	modint<mod> &operator++(int){return *this+=1;}
};
const int N=4e6+5;

const int mod=998244353;
const modint<mod> GG=3,Ginv=modint<mod>(1)/3,I=86583718;
struct poly{
	vector<modint<mod>>a;
	modint<mod>&operator[](int i){return a[i];}
	int size(){return a.size();}
	void resize(int n){a.resize(n);}
};
int rev[N];
inline poly one(){poly a;a.a.push_back(1);return a;}

基础快速变幻

inline int ext(int n){int k=0;while((1<<k)<n)k++;return k;}
inline void init(int k){int n=1<<k;for(int i=0;i<n;i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));}
inline void ntt(poly&a,int k,int typ){
	int n=1<<k;
	for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int mid=1;mid<n;mid<<=1){
		modint<mod> wn=(typ>0?GG:Ginv)^((mod-1)/(mid<<1));
		for(int r=mid<<1,j=0;j<n;j+=r){
			modint<mod> w=1;
			for(int k=0;k<mid;k++,w=w*wn){
				modint<mod> x=a[j+k],y=w*a[j+k+mid];
				a[j+k]=x+y,a[j+k+mid]=x-y;
			}
		}
	}
	if(typ<0){
		modint<mod> inv=modint<mod>(1)/n;
		for(int i=0;i<n-1;i++)a[i]*=inv;
	}
}

多项式加、减、乘

poly operator +(poly a,poly b){
	int n=max(a.size(),b.size());a.resize(n),b.resize(n);
	for(int i=0;i<n;i++)a[i]+=b[i];return a;
}
poly operator -(poly a,poly b){
	int n=max(a.size(),b.size());a.resize(n),b.resize(n);
	for(int i=0;i<n;i++)a[i]-=b[i];return a;
}
inline poly operator*(poly a,poly b){
	int n=a.size()+b.size()-1,k=ext(n);
	a.resize(1<<k),b.resize(1<<k),init(k);
	ntt(a,k,1);ntt(b,k,1);for(int i=0;i<(1<<k);i++)a[i]*=b[i];
	ntt(a,k,-1),a.resize(n);return a;
}
inline poly operator*(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]*=b;return a; }
inline poly operator/(poly a,modint<mod> b){for(int i=0;i<a.size();i++)a[i]/=b;return a; }
inline poly operator-(poly a){for(int i=0;i<a.size();i++)a[i]=-a[i];return a; }

多项式求逆

如果多项式(F)只有一项，那么显然(G_0)就是(F_0)的逆元。

若有(n)项，递归求解。

假如我们已知(F(x)H(x) equiv 1 pmod{x^{lceil frac{n}{2} ceil}})

又显然(F(x)G(x) equiv 1 pmod{x^{lceil frac{n}{2} ceil}})

那么(F(x)(G(x)-H(x))equiv 0 pmod{x^{lceil frac{n}{2} ceil}})

即(G(x)-H(x)equiv 0 pmod{x^{lceil frac{n}{2} ceil}})

两边同时平方。

(G(x)^2+H(x)^2-2G(x)H(x) equiv 0 pmod{x^n})

两边同时乘(F(x))，再由(F(x)G(x) equiv 1 pmod{x^n})，得出：

(G(x) equiv 2H(x)-F(x)H(x)^2 pmod{x^n})

poly inv(poly F,int k){
	int n=1<<k;F.resize(n);
	if(n==1){F[0]=modint<mod>(1)/F[0];return F;}
	poly G,H=inv(F,k-1);
	G.resize(n),H.resize(n<<1),F.resize(n<<1);
	for(int i=0;i<n/2;i++)G[i]=H[i]*2;
	init(k+1),ntt(H,k+1,1),ntt(F,k+1,1);
	for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i]*F[i];
	ntt(H,k+1,-1),H.resize(n);
	for(int i=0;i<n;i++)G[i]-=H[i];return G;
}
inline poly inv(poly a){
	int n=a.size();
	a=inv(a,ext(n)),a.resize(n);return a;;
}

求导数、原函数

((x^a)^prime=ax^{a-1})

(int x^adx=frac{1}{a+1}x^{a+1})

inline poly deriv(poly a){//求导 
	int n=a.size()-1;
	for(int i=0;i<n;i++)a[i]=a[i+1]*(i+1);
	a.resize(n);return a;
}
inline poly inter(poly a){//求原 
	int n=a.size()+1;a.resize(n);
	for(int i=n;i>=1;i--)a[i]=a[i-1]/i;
	a[0]=0;return a;
}

多项式(ln)

(G(x)=F(A(x)),F(x)=ln(x))

两边同时求导，得到：

(G(x)^prime=F^prime(A(x))A^prime(x)=frac{A^prime(x)}{A(x)})

在求一遍逆即可。

inline poly ln(poly a){
	int n=a.size();
	a=inter(deriv(a)*inv(a));
	a.resize(n);return a;
}

多项式(exp)

计算(F(x)=e^{A(x)})

变形得(ln F(x)-A(x)=0)

设(G(F(x))=lnF(x)-A(x))，那么就是要求这一个函数的零点。

那么我们把(F(x))看做变量，(A(x))看做常数，对这个进行求导，得(G^prime(F(x))=frac{1}{F(x)})

那么代入牛顿迭代的公式得(F(x)equiv F_0(x)-frac{G(F_0(x))}{G'(F_0(x))}pmod{x^n})

(F(x)equiv F_0(x)(1-ln F_0(x)+A(x))pmod{x^n})

(这里的(F_0(x))是指在(mod x^{frac n 2})下的答案)

然后因为(A(0)=0)，所以(F(x))的常数项为(1)

poly exp(poly a,int k){
	int n=1<<k;a.resize(n);
	if(n==1)return one();
	poly f0=exp(a,k-1);f0.resize(n);
	return f0*(one()+a-ln(f0)); 
}
poly exp(poly a){
	int n=a.size();
	a=exp(a,ext(n));a.resize(n);return a;
}

多项式开根

设(H^2(x)equiv F(x)(pmod x^{frac n2}))

(G(x)equiv H(x)(pmod x^{frac n2}))

(G(x)-H(x)equiv 0(pmod x^{frac n2}))

((G(x)-H(x))^2equiv 0(pmod x^{frac n2}))

(G^2(x)-2H(x)*G(x)+H^2(x)equiv 0(pmod x^n))

(F(x)-2H(x)*G(x)+H^2(x)equiv 0(pmod x^n))

(G(x)=frac {F(x)+H^2(x)}{2H(x)})

(a_0=1)

poly sqrt(poly F,int k){
	int n=1<<k;F.resize(n);
	if(n==1){F[0]=1;return F;}
	poly G,H=sqrt(F,k-1);
	H.resize(n);init(k);poly invH=inv(H);
	G.resize(n),H.resize(n<<1),F.resize(n<<1);
	init(k+1),ntt(H,k+1,1);
	for(int i=0;i<(n<<1);i++)H[i]=H[i]*H[i];
	ntt(H,k+1,-1),H.resize(n);
	for(int i=0;i<n;i++)G[i]=H[i]+F[i];
	G=G*invH;G.resize(n);for(int i=0;i<n;i++)G[i]/=2;
	return G;
}
inline poly sqrt(poly a){
	int n=a.size();
	a=sqrt(a,ext(n)),a.resize(n);return a;;
}

(a_0≠1)留坑

多项式快速幂

(a_0=1)

(G(x)=(F(x))^kpmod{x^n})

(ln G(x)=kF(x)pmod{x^n})

(G(x)=e^{kF(x)}pmod{x^n})

inline poly pow(poly a,modint<mod> k){//保证a[0]=1 
	a=ln(a);for(int i=0;i<a.size();i++)a[i]*=k.x;
	return exp(a);
}

(a_0≠1)

把最低项改为(1)

(F(x)^k=left(frac{F(x)}{ax^t} ight)^k(ax)^{tk})

实现了一个modpair，因为(k)既要对(mod)取模，也要对(varphi(mod))取模

struct modpair{
	//用于快速幂中的次数
	modint<mod>k1;modint<mod-1>k2;
	struct trueint{
		double lg;int x;
 		trueint &operator=(int o){return x=o,lg=log10(o),*this;}
 		trueint &operator+=(int o){return lg<=8&&(x+=o,lg=log10(x)),*this;}
		trueint &operator*=(int o){return x*=o,lg+=log10(o),*this;}
	}k3;
	modpair(modint<mod>_1,modint<mod-1>_2,trueint _3):k1(_1),k2(_2),k3(_3){}
	modpair(int o=0){k1=o,k2=o,k3=o;}
    modpair &operator = (int o){return k1=o,k2=o,k3=o,*this;}
    modpair &operator +=(int o){return k1+=o,k2+=o,k3+=o,*this;}
    modpair &operator *=(int o){return k1*=o,k2*=o,k3*=o,*this;}
    friend modpair operator +(modpair a,int o){return a+=o;}
    friend modpair operator *(modpair a,int o){return a*=o;}
    modpair operator-(){return modpair(k1,k2,k3);}
};
inline poly pow2(poly a,modpair m){//不保证a[0]=1 
	int k=0;modint<mod> val;
	while(a[k]==0&&k<a.size())k++;
	if(k==a.size()||k!=0&&m.k3.lg>8||1ll*m.k3.x*k>=a.size()){
	for(int i=0;i<a.size();i++)a[i]=0;return a;}//bye~
	val=a[k];poly b;b.resize(a.size()-k);
	for(int i=0;i<b.size();i++)b[i]=a[i+k]/val;
	b=pow(b,m.k1);for(int i=0;i<a.size();i++)a[i]=0;
	for(int i=0;i<b.size()&&i+k*m.k1.x<a.size();i++)
		a[i+k*m.k1.x]=b[i]*(val^m.k2.x);
	return a;
}

多项式三角函数

由欧拉公式得到：

[egin{cases}e^{iF(x)}=cos(F(x))+isin(F(x))\ e^{-iF(x)}=cos(F(x))-isin(F(x)) end{cases}]

解方程得到

[cos F(x)=frac{e^{iF(x)}+e^{-iF(x)}}{2} ]

[sin F(x)=frac{e^{iF(x)}-e^{-iF(x)}}{2i} ]

在取模条件下，(i=sqrt{-1}=sqrt{mod-1})，即(mod-1)的二次剩余

inline poly sin(poly a){return (exp(a*I)-exp(-a*I))/(I*2);}
inline poly cos(poly a){return (exp(a*I)+exp(-a*I))/2;}

多项式反三角函数

(G(x)equivsin^{-1}F(x)space ( ext{mod }x^n))

(G'(x)equiv frac{F'(x)}{sqrt{1-F(x)^2}}space ({ ext{mod }x^n}))

(G(x)equiv int frac{F'(x)}{sqrt{1-F(x)^2}} ext dxspace ({ ext{mod }x^n}))

(H(x)equiv an^{-1}F(x)space ( ext{mod }x^n))

(H'(x)equivfrac{F'(x)}{1+F(x)^2}space ( ext{mod }x^n))

(H(x)equivintfrac{F'(x)}{1+F(x)^2} ext dxspace ( ext{mod }x^n))

inline poly asin(poly a){
	poly G=-a*a;G[0]++;
	return inter(deriv(a)*inv(sqrt(G)));
}
inline poly atan(poly a){
	poly G=a*a;G[0]++;
	return inter(deriv(a)*inv(G));
}

多项式黑科技

任意模数NTT(MTT)

把两个多项式 (A(x),B(x)) 分别拆成 (A_0(x),A_1(x),B_0(x),B_1(x))后，求出它们的点值表示。

构造两个多项式：

[P(x)=A(x)+iB(x) ]

[Q(x)=A(x)+iB(x) ]

由于(A,B)的虚部都为(0)，(P,Q)的每一项系数都互为共轭。对(P)做一遍DFT,可以std::conj在(O(n))求出(Q)的点值表示。然后(A_0,A_1,B_0,B_1)的点值就是有手就行了

接下求(A,B)之间的两两乘积，乘出来对(A_0B_0,A_1B_0,A_0B_0,A_0B_1)做IDFT。

由于虚部不为0，但仍然可以借用上述的方法。构造两个多项式：

[P(x)=A_0(x)B_0(x)+iA_1(x)B_0(x) ]

[Q(X)=A_0(x)B_1(x)+iA_1(x)B_1(x) ]

分别做IDFT，由于 (A_0B_0,A_0B_1,A_1B_0,A_1B_1)这四个多项式卷起来后的系数表示中虚部一定为(0)，那么(P)的实部与虚部就分别为(A_0(x)B_0(x))和(A_1(x)B_0(x))，而(Q)就位(A_0(x)B_1(x))与(A_1(x)B_1(x))

给出部分代码：

typedef std::complex<double>complex;
const int N=4e6+10;const double PI=acos(-1);const complex I=complex(0,1);
int rev[N];complex Wn[N];int M,mod;
int ksm(int x,int y) {
	int re=1;
	for(;y;y>>=1,x=1LL*x*x%mod)if(y&1)re=1LL*re*x%mod;
	return re;
}
inline long long num(complex x){double d=x.real();return d<0?(long long)(d-0.5)%mod:(long long)(d+0.5)%mod;}
struct poly{
	std::vector<complex>a0,a1;
	int size(){return a0.size();}
    void resize(int n){a0.resize(n);a1.resize(n);}
	void set(int x,long long y){
		y%=mod;
		a0[x]=y/M;
		a1[x]=y%M;
	}
	long long get(int x){return (M*M*num(a0[x].real())%mod +
				M*(num(a0[x].imag())+num(a1[x].real()))%mod+num(a1[x].imag()))%mod;}
	long long val(int x){return (long long)(M*a0[x].real()+a1[x].real()+mod)%mod;}
};
poly operator+(poly a,poly b){
    int n=std::max(a.size(),b.size());a.resize(n),b.resize(n);
    for(int i=0;i<n;i++)a.set(i,a.val(i)+b.val(i));return a;
}
poly operator-(poly a,poly b){
    int n=std::max(a.size(),b.size());a.resize(n),b.resize(n);
    for(int i=0;i<n;i++)a.set(i,a.val(i)-b.val(i));return a;
}
inline poly one(){poly a;a.resize(1);a.set(0,1);return a;}
inline int ext(int n){int k=0;while((1<<k)<n)k++;return k;}
inline void init(int k){
	int n=1<<k;
	for(int i=0;i<n;i++)
		rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
	for(int i=0;i<n;i++)
		Wn[i]={cos(PI/n*i),sin(PI/n*i)};
}
void FFT(std::vector<complex>&A,int n,int t){
	if(t<0)for(int i=1;i<n;i++)if(i<(n-i))std::swap(A[i],A[n-i]);
	for(int i=0;i<n;i++)
		if(i<rev[i])std::swap(A[i],A[rev[i]]);
	for(int m=1;m<n;m<<=1)
		for(int i=0;i<n;i+=m<<1)
			for(int k=i;k<i+m;k++){
				complex W=Wn[1ll*(k-i)*n/m];
				complex a0=A[k],a1=A[k+m]*W;
				A[k]=a0+a1;A[k+m]=a0-a1; 
			}
	if(t<0)for(int i=0;i<n;i++)A[i]/=n;
}
void MTT(poly &A,int n,int t){
	for(int i=0;i<n;i++)A.a0[i]=A.a0[i]+I*A.a1[i];
	FFT(A.a0,n,t);
	for(int i=0;i<n;i++)A.a1[i]=std::conj(A.a0[i?n-i:0]);
	for(int i=0;i<n;i++){
		complex p=A.a0[i],q=A.a1[i];
		A.a0[i]=(p+q)*0.5;A.a1[i]=(q-p)*0.5*I;
	}
}
inline poly operator*(poly a,poly b){
    int n=a.size()+b.size()-1,k=ext(n);
    a.resize(1<<k),b.resize(1<<k),init(k);
    MTT(a,1<<k,1);MTT(b,1<<k,1);
	for(int i=0;i<(1<<k);i++){
		complex p=a.a0[i]*b.a0[i]+I*a.a1[i]*b.a0[i];
		complex q=a.a0[i]*b.a1[i]+I*a.a1[i]*b.a1[i];
		a.a0[i]=p,a.a1[i]=q;
	}
    FFT(a.a0,1<<k,-1);FFT(a.a1,1<<k,-1);a.resize(n);
	long long tmp;for(int i=0;i<n;i++)
		tmp=a.get(i),a.set(i,tmp);
	return a;
}
inline poly deriv(poly a){//求导 
    int n=a.size()-1;
    for(int i=0;i<n;i++)a.set(i,a.val(i+1)*(i+1));
    a.resize(n);return a;
}
inline poly inter(poly a){//求原 
    int n=a.size()+1;a.resize(n);
    for(int i=n;i>=1;i--)a.set(i,a.val(i-1)*ksm(i,mod-2));
    a.set(0,0);return a;
}
poly inv(poly F,int k){
    int n=1<<k;F.resize(n);
    if(n==1){F.set(0,ksm(F.val(0),mod-2));return F;}
    poly G,H=inv(F,k-1);
    G.resize(n),H.resize(n<<1),F.resize(n<<1);
    for(int i=0;i<n/2;i++)G.set(i,H.val(i)*2);
    H=H*H;H.resize(n);H=H*F;H.resize(n);
    G=G-H;return G;
}
inline poly inv(poly a){
    int n=a.size();
    a=inv(a,ext(n)),a.resize(n);return a;;
}
inline poly ln(poly a){
    int n=a.size();
    a=inter(deriv(a)*inv(a));
    a.resize(n);return a;
}

别问我vector太慢怎么办下次一定改