Learning：多项式（三）一些NTT的扩展

任意模数FFT

如果模数不是一些NTT的模数，那我们又如何解决呢？（当然用拆位FFT也可以做，但是我没写过）
首先我们知道有个东西叫中国剩余定理。
我们可以选取多个模数先用NTT求一下，然后再用中国剩余定理合并。一般取三个模数。
但是要注意的一点，我们在合并的时候可能会炸long long，所以我们需要使用快速乘来解决。
模板：

#include<cstdio>
#include<cctype>
#include<algorithm>
using namespace std;
const long long p[3]={998244353,1004535809,469762049},g=3,gi[3]={332748118,334845270,156587350},M=p[0]*p[1];
int n,m,mod,limit=1,l,r[4000010];
long long a[3][400010],b[3][400010];
int rd(){
    int x=0;
    char c;
    do c=getchar();
    while(!isdigit(c));
    do{
        x=(x<<1)+(x<<3)+(c^48);
        c=getchar();
    }while(isdigit(c));
    return x;
}
long long qpow(long long x,int n,long long mod){
    long long ans=1;
    while(n){
        if(n&1)
            ans=(ans*x)%mod;
        x=(x*x)%mod;
        n>>=1;
    }
    return ans;
}
void NTT(int limit,long long *a,int type,int idx){
    for(int i=0;i<limit;i++)
        if(i<r[i])
            swap(a[i],a[r[i]]);
    for(int mid=1;mid<limit;mid<<=1){
        long long wn;
        if(type==1)
            wn=qpow(g,(p[idx]-1)/(mid<<1),p[idx]);
        else wn=qpow(gi[idx],(p[idx]-1)/(mid<<1),p[idx]);
        for(int i=0;i<limit;i+=(mid<<1)){
            long long w=1;
            for(int j=0;j<mid;j++,w=w*wn%p[idx]){
                int x=a[i+j],y=w*a[i+j+mid]%p[idx];
                a[i+j]=(x+y)%p[idx];
                a[i+j+mid]=(x-y+p[idx])%p[idx];
            }
        }
    }
    if(~type)
        return;
    long long inv=qpow(limit,p[idx]-2,p[idx]);
    for(int i=0;i<=n+m;i++)
        a[i]=(a[i]*inv)%p[idx];
    return;
}
long long mul(long long a,long long b,long long mod){
    a%=mod;
    b%=mod;
    long long d=(long long)(a*(double)b/mod+0.5);
    long long ret=a*b-d*mod;
    if(ret<0)
        return ret+mod;
    return ret;
}
int main(){
    n=rd();
    m=rd();
    mod=rd();
    for(int i=0;i<=n;i++){
        a[0][i]=a[1][i]=a[2][i]=rd();
        a[0][i]%=p[0];
        a[1][i]%=p[1];
        a[2][i]%=p[2];
    }
    for(int i=0;i<=m;i++){
        b[0][i]=b[1][i]=b[2][i]=rd();
        b[0][i]%=p[0];
        b[1][i]%=p[1];
        b[2][i]%=p[2];
    }
    while(limit<=n+m){
        limit<<=1;
        l++;
    }
    for(int i=0;i<limit;i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<3;i++){
        NTT(limit,a[i],1,i);
        NTT(limit,b[i],1,i);
        for(int j=0;j<limit;j++)
            a[i][j]=(a[i][j]*b[i][j])%p[i];
        NTT(limit,a[i],-1,i);
    }
    long long ni1=qpow(p[1]%p[0],p[0]-2,p[0]),ni2=qpow(p[0]%p[1],p[1]-2,p[1]),ni3=qpow(M%p[2],p[2]-2,p[2]);
    for(int i=0;i<=n+m;i++){
        long long A=(mul(a[0][i]*p[1],ni1,M)+mul(a[1][i]*p[0],ni2,M))%M,k=((a[2][i]-A)%p[2]+p[2])%p[2]*ni3%p[2];
        printf("%lld ",(k%mod*(M%mod)%mod+A%mod)%mod);
    }
    return 0;
}

注意：以下部分代码中的$n$均已补成$2$的幂了。

多项式求逆

已知多项式$A(x)$，求一个多项式$B(x)$，满足：$A(x)B(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$
我们先考虑$n=1$的情况，此时直接将多项式$A(x)$的常数项求个逆元即可。
假设我们已经求出了满足$A(x)B(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\frac{n}{2}})$的$B(x)$。
移项得：$0\equiv 1-A(x)B(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\frac{n}{2}})$
两边同时平方得：$0\equiv 1-2A(x)B(x)+A^2(x)B^2(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$
移项得：$2A(x)B(x)-A^2(x)B^2(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$
提公因式得：$A(x)(2B(x)-A(x)B^2(x))\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$
我们用NTT求出$2B(x)-A(x)B^2(x)$，就可以得到满足$A(x)B(x)\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$的了$B(x)$。
递归求解即可。
代码：

void Inverse(int *a,int *b,int n){
	if(n==1){
		b[0]=qpow(a[0],mod-2);
		return;
	}
	Inverse(a,b,n>>1);
	for(int i=0;i<(n<<1);i++)
		r[i]=(r[i>>1]>>1)|(i&1?n:0);
	for(int i=0;i<n;i++)
		c[i]=a[i];
	for(int i=n;i<(n<<1);i++)
		c[i]=0;
	NTT(b,n<<1,1);
	NTT(c,n<<1,1);
	for(int i=0;i<(n<<1);i++)
		b[i]=(2ll*b[i]-1ll*b[i]*b[i]%mod*c[i]%mod+mod)%mod;
	NTT(b,n<<1,-1);
	for(int i=n;i<(n<<1);i++)
		b[i]=0;
	return;
}

多项式除法

已知$n$次多项式$A(x)$和$m$次多项式$B(x)$，求满足$A(x)\equiv B(x)D(x)+R(x)$的多项式$D(x)$和$R(x)$。其中多项式$D(x)$的次数为$n-m$，多项式$R(x)$的次数严格小于$m$。
对于一个多项式$A(x)$，我们记$A^r(x)$表示$A(x)$的系数反转后的多项式（即如果$A(x)={\sum }_{i=0}^n{a}_i{x}^i$，则$A^r(x)={\sum }_{i=0}^n{a}_{n-i}{x}^i$。那么我们会发现：$x^{n}A(\frac{1}{x})=A^r(x)$。
我们把$\frac{1}{x}$代入原式就可以得到：$A(\frac{1}{x})=B(\frac{1}{x})D(\frac{1}{x})+R(\frac{1}{x})$
然后就是：$$x^{n}A(\frac{1}{x})=x^{n}B(\frac{1}{x})D(\frac{1}{x})+x^{n}R(\frac{1}{x})$$
$$x^{n}A(\frac{1}{x})=x^{m}B(\frac{1}{x})x^{n-m}D(\frac{1}{x})+x^{n}R(\frac{1}{x})$$
$$A^r(x)=B^r(x)D^r(x)+x^{n-m+1}{R}^r(x)$$
因为$D(x)$的次数为$n-m$，所以我们可以在模$x^{n-m+1}$的意义下来搞，那么：
$$A^r(x)\equiv B^r(x)D^r(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{n-m+1})$$
我们会发现这个时候我们只需要一个求逆就可以求出$D^r(x)$了，然后我们反转一下代入原式就可以求出$R(x)$了。
代码：

void Divide(int *a,int N,int *b,int M,int *d){
	for(int i=0;i<=N;i++)
		tmp2[i]=a[i];
	for(int i=0;i<=M;i++)
		tmp3[i]=b[i];
	reverse(tmp2,tmp2+N+1);
	reverse(tmp3,tmp3+M+1);
	int limit=1;
	while(limit<2*(N-M+1))
		limit<<=1;
	for(int i=N-M+1;i<limit;i++)
		tmp2[i]=0;
	for(int i=N-M+1;i<limit;i++)
		tmp3[i]=0;
	Inverse(tmp3,d,limit);
	for(int i=N-M+1;i<limit;i++)
		d[i]=0;
	for(int i=0;i<limit;i++)
		r[i]=(r[i>>1]>>1)|(i&1?(limit>>1):0);
	NTT(tmp2,limit,1);
	NTT(d,limit,1);
	for(int i=0;i<limit;i++)
		d[i]=1ll*d[i]*tmp2[i]%mod;
	NTT(d,limit,-1);
	for(int i=N-M+1;i<limit;i++)
		d[i]=0;
	reverse(d,d+N-M+1);
	return;
}
void Mod(int *a,int N,int *b,int M,int *d,int &n){
	if(N<M){
		n=N;
		for(int i=0;i<=n;i++)
			d[i]=a[i];
		return;
	}
	int limit=1;
	while(limit<=N)
		limit<<=1;
	Divide(a,N,b,M,tmp4);
	for(int i=0;i<=M;i++)
		tmp6[i]=b[i];
	for(int i=M+1;i<limit;i++)
		tmp6[i]=0;
	for(int i=N-M+1;i<limit;i++)
		tmp4[i]=0;
	for(int i=0;i<limit;i++)
		r[i]=(r[i>>1]>>1)|(i&1?(limit>>1):0);
	NTT(tmp4,limit,1);
	NTT(tmp6,limit,1);
	for(int i=0;i<limit;i++)
		tmp6[i]=1ll*tmp6[i]*tmp4[i]%mod;
	NTT(tmp6,limit,-1);
	n=M-1;
	for(int i=0;i<=n;i++)
		d[i]=Minus(a[i],tmp6[i]);
	while(!d[n]&&n)
		n--;
	return;
}

多项式开根

已知多项式$A(x)$，求一个多项式$B(x)$，满足$A(x)\equiv B^2(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$
我们先考虑n=1的情况，此时直接将多项式的常数项开根即可。（这里需要用到二次剩余的知识，我也不知道怎么做。但有时候常数项确定的时候，我们可以直接赋值）
假设我们已经求出了满足$A(x)\equiv B^2(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\frac{n}{2}})$的B(x)。
移项得：$A(x)-B^2(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^{\frac{n}{2}})$
两边同时平方得：$A^2(x)-2A(x)B^2(x)+B^4(x)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$
两边同时加$4A(x)B^2(x)$得：$A^2(x)+2A(x)B^2(x)+B^4(x)\equiv 4A(x)B^2(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$
两边同时除以$4B^2(x)$，左边配方得：$(\frac{B^2(x)+A(x)}{2B(x)}{)}^2\equiv A(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$
再变形可以得到：$(\frac{1}{2}B(x)+\frac{A(x)}{B(x)}{)}^2\equiv A(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$
这个时候，通过NTT，多项式求逆求出$(\frac{1}{2}B(x)+\frac{A(x)}{B(x)}{)}^2$，就可以得到满足$A(x)\equiv B^2(x)(\mathrm{m}\mathrm{o}\mathrm{d}{x}^n)$的$B(x)$了。
代码：

void Sqrt(int *a,int *b,int n){
    if(n==1){
        b[0]=1;
        return;
    }
    Sqrt(a,b,n>>1);
    for(register int i=0;i<n;i++)
        d[i]=0;
    inv(b,d,n);
    for(register int i=0;i<(n<<1);i++)
        r[i]=(r[i>>1]>>1)|(i&1?n:0);
    for(register int i=0;i<n;i++)
        c[i]=a[i];
    for(register int i=n;i<(n<<1);i++)
        c[i]=0;
    NTT(d,n<<1,1);
    NTT(c,n<<1,1);
    for(register int i=0;i<(n<<1);i++)
        d[i]=1ll*c[i]*d[i]%mod;
    NTT(d,limit,-1);
    for(register int i=0;i<n;i++)
        b[i]=1ll*(b[i]+d[i])*inv2%mod;
    for(register int i=n;i<(n<<1);i++)
        b[i]=0;
    return;
}

多项式求导

如果多项式$A(x)={\sum }_{i=0}^n{a}_i{x}^i$，那么多项式$A(x)$的导数就是$A^{\prime }(x)={\sum }_{i=0}^n(i+1)a_{i+1}{x}_i$
怎么算出来的？
把多项式的每一项分别求导就好了。
代码：

void Derivative(int *a,int *b,int n){
    for(register int i=0;i<n;i++)
        b[i]=1ll*a[i+1]*(i+1)%mod;
    b[n]=0;
    return;
}

多项式积分

如果多项式$A(x)={\sum }_{i=1}^n{a}_i{x}^i$，那么这个多项式的积分$\int A(x)dx={\sum }_{i=1}^n\frac{1}{i}{a}_{i-1}{x}^i$。
大概就是把求导的过程反过来就行了。
代码：

void Integrate(int *a,int *b,int n){
    for(register int i=n-1;i;i--)
        b[i]=1ll*a[i-1]*inv[i]%mod;
    b[0]=0;
    return;
}

多项式求ln

我们已知多项式$A(x)$，要求一个多项式$B(x)$满足$lnA(x)=B(x)$。
我们把$lnA(x)$求导得到$\frac{A^{\prime }(x)}{A(x)}$，又因为$lnA(x)=B(x)$，所以$B(x)=\int \frac{A^{\prime }(x)}{A(x)}$，于是就可以用多项式求逆+多项式求导+多项式积分搞定啦。
代码：

void Ln(int *a,int *b,int n){
    Derivative(a,b,n);
    for(int i=0;i<(n<<1);i++)
        tmp2[i]=0;
    Inverse(a,tmp2,n);
    for(int i=0;i<(n<<1);i++)
        r[i]=(r[i>>1]>>1)|(i&1?n:0);
    NTT(b,n<<1,1);
    NTT(tmp2,n<<1,1);
    for(int i=0;i<(n<<1);i++)
        b[i]=1ll*b[i]*tmp2[i]%mod;
    NTT(b,n<<1,-1);
    Integrate(b,b,n<<1);
    return;
}

留个坑，剩下的暂时没有心情写，等心情好了再填坑。

多项式exp

多项式求幂

多点求值

快速插值

这些模板先放上来：
$exp$&求幂

void Exp(int *a,int *b,int n){
    if(n==1){
        b[0]=1;
        return;
    }
    Exp(a,b,n>>1);
    for(register int i=0;i<(n<<1);i++)
        tmp3[i]=0;
    Ln(b,tmp3,n);
    for(register int i=0;i<n;i++)
        tmp3[i]=Minus(a[i],tmp3[i]);
    tmp3[0]=Add(tmp3[0],1);
    for(register int i=n;i<(n<<1);i++)
        b[i]=tmp3[i]=0;
    for(register int i=0;i<(n<<1);i++)
        r[i]=(r[i>>1]>>1)|(i&1?n:0);
    NTT(b,n<<1,1);
    NTT(tmp3,n<<1,1);
    for(register int i=0;i<(n<<1);i++)
        b[i]=1ll*b[i]*tmp3[i]%mod;
    NTT(b,n<<1,-1);
    for(register int i=n;i<(n<<1);i++)
        b[i]=0;
    return;
}
void Pow(int *a,int *b,int n,int k){
    Ln(a,tmp4,n);
    for(register int i=0;i<n;i++)
        tmp4[i]=1ll*tmp4[i]*k%mod;
    Exp(tmp4,b,n);
    return;
}

多点求值&快速插值

#include<cstdio>
#include<algorithm>
using namespace std;
const int mod=998244353,g=3,gi=332748118;
int n,x[50010],y[50010],*ans[200010],*a[200010],*b[200010],*c[200010],r[200010],val[200010],a0[200010],a1[200010],b0[200010],b1[200010],tmp1[200010],tmp2[200010],tmp3[200010],tmp4[200010],tmp5[200010],tmp6[200010];
char buffer[1000010],*head,*tail;
char Getchar(){
	if(head==tail){
		int len=fread(buffer,1,1000000,stdin);
		head=buffer,tail=head+len;
		if(head==tail)
			return EOF;
	}
	return *head++;
}
int rd(){
	int x=0;
	char c;
	do c=Getchar();
	while(!isdigit(c));
	do{
		x=(x<<1)+(x<<3)+(c^48);
		c=Getchar();
	}while(isdigit(c));
	return x;
}
int Add(int a,int b){
	return a+b>=mod?a+b-mod:a+b;
}
int Minus(int a,int b){
	return a<b?a-b+mod:a-b;
}
int qpow(int x,int n){
	int ret=1;
	while(n){
		if(n&1)
			ret=1ll*ret*x%mod;
		x=1ll*x*x%mod;
		n>>=1;
	}
	return ret;
}
void NTT(int *a,int limit,int type){
	for(int i=0;i<limit;i++)
		if(i<r[i])
			swap(a[i],a[r[i]]);
	for(int mid=1;mid<limit;mid<<=1){
		int wn=qpow(type==1?g:gi,(mod-1)/(mid<<1));
		for(int i=0;i<limit;i+=(mid<<1)){
			int w=1;
			for(int j=0;j<mid;j++,w=1ll*w*wn%mod){
				int x=a[i+j],y=1ll*w*a[i+j+mid]%mod;
				a[i+j]=Add(x,y);
				a[i+j+mid]=Minus(x,y);
			}
		}
	}
	if(type==-1){
		int inv=qpow(limit,mod-2);
		for(int i=0;i<limit;i++)
			a[i]=1ll*a[i]*inv%mod;
	}
	return;
}
void Inverse(int *a,int *b,int N){
	if(N==1){
		b[0]=qpow(a[0],mod-2);
		return;
	}
	Inverse(a,b,N>>1);
	for(register int i=0;i<N;i++)
		tmp1[i]=a[i];
	for(register int i=N;i<(N<<1);i++)
		tmp1[i]=b[i]=0;
	for(register int i=0;i<(N<<1);i++)
		r[i]=(r[i>>1]>>1)|(i&1?N:0);
	NTT(b,N<<1,1);
	NTT(tmp1,N<<1,1);
	for(register int i=0;i<(N<<1);i++)
		b[i]=Minus(Add(b[i],b[i]),1ll*tmp1[i]*b[i]%mod*b[i]%mod);
	NTT(b,N<<1,-1);
	for(register int i=N;i<(N<<1);i++)
		b[i]=0;
	return;
}
void Divide(int *a,int N,int *b,int M,int *d){
	for(int i=0;i<=N;i++)
		tmp2[i]=a[i];
	for(int i=0;i<=M;i++)
		tmp3[i]=b[i];
	reverse(tmp2,tmp2+N+1);
	reverse(tmp3,tmp3+M+1);
	int limit=1;
	while(limit<2*(N-M+1))
		limit<<=1;
	for(int i=N-M+1;i<limit;i++)
		tmp2[i]=0;
	for(int i=N-M+1;i<limit;i++)
		tmp3[i]=0;
	Inverse(tmp3,d,limit);
	for(int i=N-M+1;i<limit;i++)
		d[i]=0;
	for(int i=0;i<limit;i++)
		r[i]=(r[i>>1]>>1)|(i&1?(limit>>1):0);
	NTT(tmp2,limit,1);
	NTT(d,limit,1);
	for(int i=0;i<limit;i++)
		d[i]=1ll*d[i]*tmp2[i]%mod;
	NTT(d,limit,-1);
	for(int i=N-M+1;i<limit;i++)
		d[i]=0;
	reverse(d,d+N-M+1);
	return;
}
void Mod(int *a,int N,int *b,int M,int *d,int &n){
	if(N<M){
		n=N;
		for(int i=0;i<=n;i++)
			d[i]=a[i];
		return;
	}
	int limit=1;
	while(limit<=N)
		limit<<=1;
	Divide(a,N,b,M,tmp4);
	for(int i=0;i<=M;i++)
		tmp6[i]=b[i];
	for(int i=M+1;i<limit;i++)
		tmp6[i]=0;
	for(int i=N-M+1;i<limit;i++)
		tmp4[i]=0;
	for(int i=0;i<limit;i++)
		r[i]=(r[i>>1]>>1)|(i&1?(limit>>1):0);
	NTT(tmp4,limit,1);
	NTT(tmp6,limit,1);
	for(int i=0;i<limit;i++)
		tmp6[i]=1ll*tmp6[i]*tmp4[i]%mod;
	NTT(tmp6,limit,-1);
	n=M-1;
	for(int i=0;i<=n;i++)
		d[i]=Minus(a[i],tmp6[i]);
	while(!d[n]&&n)
		n--;
	return;
}
void Mul(int *a,int lena,int *b,int lenb,int *c){
	int limit=1;
	while(limit<=lena+lenb)
		limit<<=1;
	for(int i=0;i<=lena;i++)
		tmp1[i]=a[i];
	for(int i=lena+1;i<limit;i++)
		tmp1[i]=0;
	for(int i=0;i<=lenb;i++)
		tmp2[i]=b[i];
	for(int i=lenb+1;i<limit;i++)
		tmp2[i]=0;
	for(int i=0;i<limit;i++)
		r[i]=(r[i>>1]>>1)|(i&1?(limit>>1):0);
	NTT(tmp1,limit,1);
	NTT(tmp2,limit,1);
	for(int i=0;i<limit;i++)
		tmp2[i]=1ll*tmp2[i]*tmp1[i]%mod;
	NTT(tmp2,limit,-1);
	for(int i=0;i<=lena+lenb;i++)
		c[i]=tmp2[i];
	return;
}
void Derivative(int *a,int *b,int N){
	for(int i=0;i<N-1;i++)
		b[i]=1ll*a[i+1]*(i+1)%mod;
	b[N-1]=0;
	return;
}
void solve1(int l,int r,int o){
	if(l==r){
		a[o]=new int[2];
		a[o][0]=x[l]?mod-x[l]:x[l];
		a[o][1]=1;
		return;
	}
	int mid=(l+r)>>1;
	solve1(l,mid,o<<1);
	solve1(mid+1,r,o<<1|1);
	a[o]=new int[r-l+2];
	Mul(a[o<<1],mid-l+1,a[o<<1|1],r-mid,a[o]);
	return;
}
void solve2(int *b,int n,int l,int r,int o){
	int *now=new int[r-l+1];
	Mod(b,n,a[o],r-l+1,now,n);
	if(l==r){
		val[l]=now[0];
		return;
	}
	int mid=(l+r)>>1;
	solve2(now,n,l,mid,o<<1);
	solve2(now,n,mid+1,r,o<<1|1);
	return;
}
void solve3(int l,int r,int o){
	ans[o]=new int[r-l+1];
	if(l==r){
		ans[o][0]=y[l];
		return;
	}
	int mid=(l+r)>>1;
	solve3(l,mid,o<<1);
	solve3(mid+1,r,o<<1|1);
	Mul(ans[o<<1],mid-l,a[o<<1|1],r-mid,tmp4);
	Mul(a[o<<1],mid-l+1,ans[o<<1|1],r-mid-1,tmp3);
	for(int i=0;i<r-l+1;i++)
		ans[o][i]=Add(tmp4[i],tmp3[i]);
	return;
}
int main(){
	n=rd();
	for(int i=0;i<=n;i++)
		x[i]=rd(),y[i]=rd();
	solve1(0,n,1);
	Derivative(a[1],tmp5,n+2);
	solve2(tmp5,n,0,n,1);
	for(int i=0;i<=n;i++)
		y[i]=1ll*y[i]*qpow(val[i],mod-2)%mod;
	solve3(0,n,1);
	for(int i=0;i<=n;i++)
		printf("%d ",ans[1][i]);
	return 0;
}