多项式

正睿集训三道题考两道多项式，于是我决定补一补这个巨大的坑。

一、NTT

前置知识

阶

如果\(\gcd(a,p)=1\)，那么对于方程\(a^r\equiv1\pmod p\)，使它成立的最小的 \(r\) 称为 \(a\) 关于 \(p\) 的阶，记作 \(ord_p(a)\)

性质：对于 \(i\in[0,ord_p(a))\) ，所有的 \(a^i\space mod\space p\) 的结果互不相同。

原根

若 \(g\) 模 \(n\) 的阶为 \(\phi(n)\) ，且 \(gcd(g,n)=1\) ，则称 \(g\) 为 \(n\) 的一个原根。

令 \(n\) 为大于 \(1\) 的 \(2\) 的幂，\(p\) 为素数，且 \(p\) 能表示成 \(kn+1\) ，\(g\)为\(p\)的一个原根。

设

\[g_n=g^{\frac{p-1}{n}} \]

所以

\[g_n^n=g^{n\times \frac{p-1}{n} }=g^{p-1} \]

\[g_n^{\frac{n}{2}}=g^{\frac{p-1}{2}} \]

\[g_{an}^{ak}=g^{\frac{ak(p-1)}{an} }=g^{\frac{k(p-1)}{n}}=g_n^k \]

显然

\[g_n^n\equiv1(mod\space p) \]

\[g_n^{\frac{n}{2}}\equiv-1(mod\space p) \]

\[(g_n^{k+\frac{n}{2}})^2=g_n^{2k+n}\equiv g_n^{2k}(mod\space p) \]

由此可以发现，单位根有的性质，原根都有，因此可以用原根代替单位根。

\(code:\)

#include<iostream>
#include<cstdio>
#define int long long
using namespace std;
const int p=998244353,g=3,gi=332748118;
int n,m,lim,len,a[2100005],b[2100005],rev[2100005];
int f(int x,int y){
    int s=x;x=1;
    while(y){
        if(y&1)
            x=x*s%p;
        s=s*s%p;y>>=1;
    }
    return x%p;
}
void ntt(int *a,int lim,int tmp){
    for(int i=0;i<lim;++i)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        int wn;
        if(tmp==1)
            wn=f(g,(p-1)/(mid<<1));
        else
            wn=f(gi,(p-1)/(mid<<1));
        for(int j=0;j<lim;j+=(mid<<1)){
            int w=1;
            for(int k=0;k<mid;++k,w=w*wn%p){
                int x=a[j+k],y=w*a[j+k+mid]%p;
                a[j+k]=(x+y)%p;
                a[j+k+mid]=(x-y+p)%p;
            }
        }
    }
    if(tmp==-1){
        int inv=f(lim,p-2);
        for(int i=0;i<lim;++i)
            a[i]=a[i]*inv%p;
    }
}
signed main(){
    scanf("%lld%lld",&n,&m);//这里的n和m是次数
    for(int i=0;i<=n;++i)
        scanf("%lld",&a[i]),a[i]=(a[i]+p)%p;
    for(int i=0;i<=m;++i)
        scanf("%lld",&b[i]),b[i]=(b[i]+p)%p;
    lim=1;
    while(lim<=n+m)
        lim<<=1,++len;
    for(int i=0;i<lim;++i)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));//rev[i]表示i的二进制翻转，如rev[11001]=10011
    ntt(a,lim,1);ntt(b,lim,1);
    for(int i=0;i<lim;++i)
        a[i]=a[i]*b[i]%p;
    ntt(a,lim,-1);
    for(int i=0;i<=n+m;++i)
        printf("%lld ",a[i]);
    return 0;
}

二、牛顿迭代法

已知多项式 \(g(x)\) ，求一个 \(f(x)\) 使得方程 \(g(f(x))\equiv 0(mod\space x^n)\)成立。

结论：

当 \(n=1\) 时，方程的解需要单独求出。

假设已知 \(f_0(x)\) 使得方程 \(g(f_0(x))\equiv 0(mod\space x^{\lceil \frac{n}{2}\rceil} )\) 成立。那么

\[f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\space(mod\space x^n) \]

三、多项式求逆

已知多项式 \(g(x)\) ，求一个 \(f(x)\) 使得方程 \(g(x)-\frac{1}{f(x)}\equiv 0(mod\space x^n)\)成立。

假设我们已知 \(f_0(x)\) 使得方程 \(g(x)-\frac{1}{f_0(x)}\equiv 0(mod\space x^{\lceil \frac{n}{2}\rceil} )\) 成立。

根据牛顿迭代法，有

\[f(x)=f_0(x)-\frac{g(x)-\frac{1}{f_0(x)}}{\frac{1}{f_0^2(x)}} \]

\[f(x)=2f_0(x)-g(x)f_0^2(x) \]

\(code:\)

#include<iostream>
#include<cstdio>
using namespace std;
const long long maxn=3e6+5,p=998244353,g=3,gi=332748118;//gi为原根的逆元 
long long n,inv,a[2100005],r[2100005],b[2100005],c[2100005];
long long f(long long a,long long b){
	long long s=a;a=1;
	while(b){
		if(b&1) a=a*s%p;
		s=s*s%p;b>>=1;
	}
	return a%p;
}
void ntt(long long *a,long long lim,long long tmp){
	for(int i=0;i<lim;++i)
		if(i<r[i])
			swap(a[i],a[r[i]]);
	for(int mid=1;mid<lim;mid<<=1){
		long long wn;
		if(tmp==1)
			wn=f(g,(p-1)/(mid<<1));
		else
			wn=f(gi,(p-1)/(mid<<1));
		for(int j=0;j<lim;j+=(mid<<1)){
			long long w=1;
			for(int k=0;k<mid;++k,w=(w*wn)%p){
				long long x=a[j+k],y=w*a[j+k+mid]%p;
				a[j+k]=(x+y)%p;
				a[j+k+mid]=(x-y+p)%p;
			}
		}
	}
	if(tmp==-1){
		inv=f(lim,p-2);
		for(int i=0;i<lim;++i)
			a[i]=a[i]*inv%p;
	}
}
void work(long long n,long long *a,long long *b){
	if(n==1){
		b[0]=f(a[0],p-2);
		return ;
	}
	work((n+1)>>1,a,b);
	long long lim=1,l=0;
	while(lim<(n<<1))
		lim<<=1,++l;
	for(int i=0;i<lim;++i)
		r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(int i=0;i<n;++i)
		c[i]=a[i];
	for(int i=n;i<lim;++i)
		c[i]=0;
	ntt(c,lim,1);ntt(b,lim,1);
	for(int i=0;i<lim;++i)
		b[i]=(2-c[i]*b[i]%p+p)%p*b[i]%p;//B=2B'-AB'^2
	ntt(b,lim,-1);
	for(int i=n;i<lim;++i)
		b[i]=0;
}
int main(){
	scanf("%lld",&n);
	for(int i=0;i<n;++i)
		scanf("%lld",&a[i]),a[i]=(a[i]+p)%p;
	work(n,a,b);
	for(int i=0;i<n;++i)
		printf("%lld ",b[i]); 
	return 0;
}

四、多项式开方

已知多项式 \(g(x)\) ，求一个 \(f(x)\) 使得方程 \(g(x)-f^2(x)\equiv 0(mod\space x^n)\)成立。

假设我们已知 \(f_0(x)\) 使得方程 \(g(x)-f_0^2(x)\equiv 0(mod\space x^{\lceil \frac{n}{2}\rceil} )\) 成立。

根据牛顿迭代法，有

\[f(x)=f_0(x)-\frac{g(x)-f_0^2(x)}{-2f_0(x)} \]

\[f(x)=\frac{f_0^2(x)+g(x)}{2f_0(x)} \]

这里有分式，套用多项式求逆即可。

\(code:\)

#include<iostream>
#include<cstdio>
using namespace std;
const long long maxn=3e6+5,P=998244353,G=3,Gi=332748118,inv2=499122177;//gi为原根的逆元
long long n,m,inv,r[500005],c[500005],f[500005],g[500005],d[500005];
long long ginv[500005],gr[500005],fr[500005];
long long _pow(long long a,long long b){
	long long s=a;a=1;
	while(b){
		if(b&1) a=a*s%P;
		s=s*s%P;b>>=1;
	}
	return a%P;
}
void ntt(long long *a,long long lim,long long tmp){
	for(int i=0;i<lim;++i)
		if(i<r[i])
			swap(a[i],a[r[i]]);
	for(int mid=1;mid<lim;mid<<=1){
		long long wn;
		if(tmp==1)
			wn=_pow(G,(P-1)/(mid<<1));
		else
			wn=_pow(Gi,(P-1)/(mid<<1));
		for(int j=0;j<lim;j+=(mid<<1)){
			long long w=1;
			for(int k=0;k<mid;++k,w=(w*wn)%P){
				long long x=a[j+k],y=w*a[j+k+mid]%P;
				a[j+k]=(x+y)%P;
				a[j+k+mid]=(x-y+P)%P;
			}
		}
	}
	if(tmp==-1){
		inv=_pow(lim,P-2);
		for(int i=0;i<lim;++i)
			a[i]=a[i]*inv%P;
	}
}
void finv(long long n,long long *a,long long *b){
	if(n==1){
		b[0]=_pow(a[0],P-2);
		return ;
	}
	finv((n+1)>>1,a,b);
	long long lim=1,l=0;
	while(lim<(n<<1))
		lim<<=1,++l;
	for(int i=0;i<lim;++i)
		r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(int i=0;i<n;++i)
		c[i]=a[i];
	for(int i=n;i<lim;++i)
		c[i]=0;
	ntt(c,lim,1);ntt(b,lim,1);
	for(int i=0;i<lim;++i)
		b[i]=(2-c[i]*b[i]%P+P)%P*b[i]%P;//B=2B'-AB'^2
	ntt(b,lim,-1);
	for(int i=n;i<lim;++i)
		b[i]=0;
}
void ssqrt(long long n,long long *a,long long *b){
	if(n==1){
		b[0]=1;return ;
	}
	ssqrt((n+1)>>1,a,b);
	long long lim=1,l=0;
	while(lim<(n<<1))
		lim<<=1,++l;
	for(int i=0;i<lim;++i)
		d[i]=0;
	finv(n,b,d);
	for(int i=0;i<n;++i)
		c[i]=a[i];
	for(int i=n;i<lim;++i)
		c[i]=0;
	for(int i=0;i<lim;++i)
		r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	ntt(c,lim,1);ntt(d,lim,1);
	for(int i=0;i<lim;++i)
		c[i]=(c[i]*d[i])%P;//B=2B'-AB'^2
	ntt(c,lim,-1);
	for(int i=0;i<n;++i)
		b[i]=(b[i]+c[i])%P*inv2%P;
	for(int i=n;i<lim;++i)
		b[i]=0;
}
int main(){
	scanf("%lld",&n);
	for(int i=0;i<n;++i)
		scanf("%lld",&f[i]);
	ssqrt(n,f,g);
	for(int i=0;i<n;++i)
		printf("%lld ",g[i]);
	printf("\n");
	return 0;
}

五、多项式ln（多项式对数函数）

已知多项式 \(g(x)\) ，求一个 \(f(x)\) 使得方程 \(ln\space g(x)\equiv\space f(x)(mod\space x^n)\)成立。

可以先求出\(f(x)\)的导数\(f'(x)\)，再求\(f'(x)\)的积分，得到\(f(x)\)。根据复合函数的求导法则，有\(f'(x)=g'(x)g^{-1}(x)\) 。因此只需要求\(g(x)\)的逆，和\(g(x)\)的导数即可。然后再求\(f'(x)\)的积分即可。

\(code:\)

#include<iostream>
#include<cstdio>
#define int long long
using namespace std;
const int mod=998244353,g=3,gi=332748118;
int n,lim,len,a[400005],b[400005],c[400005],r[400005];
int fpow(int x,int y){
    int s=x;x=1;
    while(y){
        if(y&1)
            x=x*s%mod;
        y>>=1;s=s*s%mod;
    }
    return x%mod;
}
void ntt(int *a,int lim,int tmp){
    for(int i=0;i<lim;++i)
        if(i<r[i])
            swap(a[i],a[r[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        int wn;
        if(tmp==1)
            wn=fpow(g,(mod-1)/(mid<<1));
        else
            wn=fpow(gi,(mod-1)/(mid<<1));
        for(int j=0;j<lim;j+=mid<<1){
            int w=1;
            for(int k=0;k<mid;++k,w=w*wn%mod){
                int x=a[j+k],y=a[j+k+mid]*w%mod;
                a[j+k]=(x+y)%mod;
                a[j+k+mid]=(x-y+mod)%mod;
            }
        }
    }
    if(tmp==-1){
        int inv=fpow(lim,mod-2);
        for(int i=0;i<lim;++i)
            a[i]=a[i]*inv%mod;
    }
}
void work(int n,int *a,int *b){
    if(n==1){
        b[0]=fpow(a[0],mod-2);
        return ;
    }
    work((n+1)>>1,a,b);
    for(int i=0;i<n;++i)
        c[i]=a[i];
    for(int i=n;i<lim;++i)
        c[i]=0;
    ntt(c,lim,1);ntt(b,lim,1);
    for(int i=0;i<lim;++i)
        b[i]=(2-c[i]*b[i]%mod+mod)%mod*b[i]%mod;
    ntt(b,lim,-1);
    for(int i=n;i<lim;++i)
        b[i]=0;
}
void mul(int n,int *a,int *b){
    ntt(a,lim,1);ntt(b,lim,1);
    for(int i=0;i<lim;++i)
        a[i]=a[i]*b[i]%mod;
    ntt(a,lim,-1);
    for(int i=n;i<lim;++i)
        a[i]=0;
}
void dx(int n,int *a,int *b){
    b[n-1]=0;
    for(int i=n-1;i>=1;--i)
        b[i-1]=a[i]*i%mod;
}
void sum(int n,int *a,int *b){
    b[0]=0;
    for(int i=1;i<n;++i)
        b[i]=fpow(i,mod-2)*a[i-1]%mod;
}
signed main(){
    scanf("%lld",&n);
    for(int i=0;i<n;++i)
        scanf("%lld",&a[i]),a[i]=(a[i]+mod)%mod;
    lim=1;len=0;
    while(lim<(n<<1))
        lim<<=1,++len;
    for(int i=0;i<lim;++i)
        r[i]=(r[i>>1]>>1)|((i&1)<<(len-1));
    work(n,a,b);
    dx(lim,a,c);
    mul(n,b,c);
    sum(n,b,c);
    for(int i=0;i<n;++i)
        printf("%lld ",c[i]);
    printf("\n");
    return 0;
}

六、多项式exp(多项式指数函数)

已知多项式 \(g(x)\) ，求一个 \(f(x)\) 使得方程 \(g(x)-ln\space f(x)\equiv 0(mod\space x^n)\)成立。

假设我们已知 \(f_0(x)\) 使得方程 \(ln\space f_0(x)-g(x)\equiv 0(mod\space x^n)\)成立。

根据牛顿迭代法，有

\[f(x)=f_0(x)-\frac{ln\space f_0(x)-g(x)}{\frac{1}{f_0(x)}} \]

\[f(x)=f_0(x)(1-ln\space f_0(x)+g(x)) \]

七、多项式幂函数

已知多项式 \(f(x)\) ，求一个 \(g(x)\) 使得方程 \(f(x)^k\equiv g(x)(mod\space x^n)\)成立。

将等式两边取对数，得到：

\[k\space ln\space f(x)\equiv ln\space g(x) \]

因此只要先求一次对数，再将每一个系数乘以 \(k\) ，再求一次 \(exp\) 即可。

加强版：不保证 \(f_0=1\)

找到第一个不为 \(0\) 的项，然后提出它的公因式即可。

void power(int n,int k,int k2,int *a,int *b){
    int lim=1,len=0,inv1=1,pos=-1;
    while(lim<(n<<1))
        lim<<=1,++len;
    for(int i=0;i<lim;++i)
        r[i]=(r[i>>1]>>1)|((i&1)<<(len-1));
    for(int i=0;i<n;++i){
        if(a[i]!=0&&pos==-1)
            inv1=fpow(a[i],mod-2),pos=i;
        a[i]=a[i]*inv1%mod;
    }
    if(ok&&!a[0]){
        for(int i=0;i<n;++i)
            b[i]=0;
        return ;
    }
    for(int i=0;i<n;++i){
        if(i<n-pos) a[i]=a[i+pos];
        else a[i]=0;
    }
    n-=pos;
    ln(n,a,b);
    for(int i=0;i<n;++i)
        a[i]=b[i]*k%mod;
    for(int i=0;i<lim;++i){
        if(i>=n) a[i]=0;
        b[i]=0;
    }
    exp(n,a,b);
    n+=pos;
    inv1=fpow(inv1,mod-2);
    inv1=fpow(inv1,k2);
    for(int i=n-1;i>=k*pos;--i)
        b[i]=b[i-k*pos]*inv1%mod;
    for(int i=0;i<min(k*pos,n);++i)
        b[i]=0;
    return ;
}