多项式

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

一、NTT

前置知识

阶

如果$gcd(a,p)=1$，那么对于方程$a^r\equiv 1(\mathrm{mod}p)$，使它成立的最小的 $r$ 称为 $a$ 关于 $p$ 的阶，记作 $or{d}_p(a)$

性质：对于 $i\in [0,or{d}_p(a))$ ，所有的 $a^{i}modp$ 的结果互不相同。

原根

若 $g$ 模 $n$ 的阶为 $\varphi (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\equiv 1(modp)$$

$$g_n^{\frac{n}{2}}\equiv -1(modp)$$

$$(g_n^{k+\frac{n}{2}}{)}^2=g_n^{2k+n}\equiv g_n^{2k}(modp)$$

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

$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{x}^n)$成立。

结论：

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

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

$$f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g^{\prime }(f_0(x))}(mod{x}^n)$$

三、多项式求逆

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

假设我们已知 $f_0(x)$ 使得方程 $g(x)-\frac{1}{f_0(x)}\equiv 0(mod{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{x}^n)$成立。

假设我们已知 $f_0(x)$ 使得方程 $g(x)-f_0^2(x)\equiv 0(mod{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)$ 使得方程 $lng(x)\equiv f(x)(mod{x}^n)$成立。

可以先求出$f(x)$的导数$f^{\prime }(x)$，再求$f^{\prime }(x)$的积分，得到$f(x)$。根据复合函数的求导法则，有$f^{\prime }(x)=g^{\prime }(x)g^{-1}(x)$ 。因此只需要求$g(x)$的逆，和$g(x)$的导数即可。然后再求$f^{\prime }(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)-lnf(x)\equiv 0(mod{x}^n)$成立。

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

根据牛顿迭代法，有

$$f(x)=f_0(x)-\frac{ln{f}_0(x)-g(x)}{\frac{1}{f_0(x)}}$$

$$f(x)=f_0(x)(1-ln{f}_0(x)+g(x))$$

七、多项式幂函数

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

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

$$klnf(x)\equiv lng(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 ;
}