「算法笔记」多项式操作

注意下面的板子可能常数较大。

2021 年写的多项式求逆（已折叠）
一、基本思路

先给出多项式求逆的定义：

多项式求逆：给定一个 \(n-1\) 次多项式 \(A(x)\)，求一个多项式 \(B(x)\) 满足：\(A(x)B(x)\equiv 1\pmod {x^n}\)。模 \(x^n\) 即只考虑多项式的 \(x^0,x^1,\cdots,x^{n-1}\) 项。

考虑倍增。假设已经求出了多项式 \(A(x)\) 在模 \(x^{\lceil \tfrac{n}{2}\rceil}\) 意义下的逆元 \(B'(x)\)，那么有：

\[A(x)B'(x) \equiv 1\pmod{x^{\lceil \frac{n}{2}\rceil}} \]

因为 \(A(x)B(x)\equiv 1\pmod {x^n}\)，所以显然也有：\(A(x)B(x)\equiv 1\pmod{x^{\lceil \tfrac{n}{2}\rceil}}\)。

将两式相减得到：\(A(x)(B(x)-B'(x))\equiv 0\pmod{x^{\lceil\tfrac{n}{2}\rceil}}\)。

由于 \(A(x)\nmid x^{\lceil \tfrac{n}{2}\rceil}\)，那么有：\(B(x)-B'(x)\equiv 0\pmod{x^{\lceil\tfrac{n}{2}\rceil}}\)。

两边同时平方得：\(B(x)^2-2B(x)B'(x)+B'(x)^2\equiv 0\pmod {x^n}\)。

移项得：\(B(x)^2\equiv 2B(x)B'(x)-B'(x)^2\pmod {x^n}\)。

在两边同时乘上 \(A(x)\) 得到：

\[A(x)B(x)^2\equiv 2A(x)B(x)B'(x)-A(x)B'(x)^2\pmod {x^n} \]

通过逆元的定义 \(A(x)B(x)\equiv 1 \pmod {x^n}\) 可以化简为：

\[B(x)\equiv 2B'(x)-A(x)B'(x)^2\pmod {x^n} \]
二、具体实现

可以 递归 求解，递归边界为 \(n=1\)，此时答案为 常数项的逆元。

也可以 迭代 实现（枚举迭代长度），常数较小。

时间复杂度 \(T(n)=T(\tfrac{n}{2})+\mathcal O(n\log n)=\mathcal O(n\log n)\)。

//Luogu P4238
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=3e6+5,mod=998244353;
int n,f[N],a[N],b[N],res[N],tmp[N],len,r[N],inv;
int mul(int x,int n,int mod){
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans;
}
void NTT(int a[N],int n,int opt){    //opt=1/-1: DFT/IDFT
    for(int i=0;i<n;i++)
        if(i<r[i]) swap(a[i],a[r[i]]);
    for(int k=2;k<=n;k<<=1){
        int m=k>>1,x=mul(3,(mod-1)/k,mod),w=1,v; 
        if(opt==-1) x=mul(x,mod-2,mod);
        for(int i=0;i<n;i+=k,w=1)
            for(int j=i;j<i+m;j++) v=w*a[j+m]%mod,a[j+m]=(a[j]-v+mod)%mod,a[j]=(a[j]+v)%mod,w=w*x%mod;
    }
    if(opt==-1){
        inv=mul(len,mod-2,mod);
        for(int i=0;i<n;i++) a[i]=a[i]*inv%mod;
    }
} 
void solve(int A[N],int B[N],int n,int m,int res[N]){
    for(len=1;len<=n+m;len<<=1); 
    for(int i=0;i<len;i++)
        r[i]=(r[i>>1]>>1)|((i&1)?len>>1:0),a[i]=b[i]=0;
    for(int i=0;i<n;i++) a[i]=A[i];
    for(int i=0;i<m;i++) b[i]=B[i];
    NTT(a,len,1),NTT(b,len,1);
    for(int i=0;i<len;i++) a[i]=((2*b[i]%mod-a[i]*b[i]%mod*b[i]%mod)%mod+mod)%mod;
    NTT(a,len,-1);
    for(int i=0;i<n;i++) res[i]=a[i];
}
void invp(int a[N],int n,int res[N]){
    if(n==1){res[0]=mul(a[0],mod-2,mod);return ;}
    invp(a,(n+1)/2,res);
    solve(a,res,n,n,res);    //注意这里写 solve(a,res,n,(n+1)/2,res) 会出错。原因：B(x)≡2B′(x)−A(x)B′(x)^2 这个地方，最高次数可以达到 2n 次，如果 solve(a,res,n,(n+1)/2,res) 这样算出来的点值个数可能不到 2n 个，所以没有办法 IDFT 得到正确的多项式。
}
signed main(){
    scanf("%lld",&n);
    for(int i=0;i<n;i++)
        scanf("%lld",&f[i]);
    invp(f,n,res);
    for(int i=0;i<n;i++)
        printf("%lld%c",res[i],i==n-1?'\n':' ');
    return 0;
}

 2022.1.28 重学写了一遍迭代的板子：

//Luogu P4238
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,a[N],b[N],len,r[N],A[N],B[N],p[N],q[N];
int qpow(int x,int n){
    int ans=1;
    for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
    return ans;
}
void NTT(int *a,int n,int op){
    for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
    for(int k=2;k<=n;k<<=1){
        int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
        for(int i=0;i<n;i+=k,w=1)
            for(int j=i;j<i+m;j++)
                x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
    }
    if(op==-1){
        int inv=qpow(n,mod-2);
        for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
    }
} 
void mul(int *a,int *b,int len){
    for(int i=0;i<len;i++)
        p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?len>>1:0);
    NTT(p,len,1),NTT(q,len,1);
    for(int i=0;i<len;i++) a[i]=1ll*p[i]*q[i]%mod;
    NTT(a,len,-1);
}
void inv(int *a,int *b,int len){
    for(int i=0;i<len;i++) b[i]=0;
    b[0]=qpow(a[0],mod-2);
    for(int i=2;i<=len;i<<=1){
        for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
        //A(x)B'(x)^2 长度最大为 i*2
        mul(B,B,i),mul(A,B,i<<1);
        for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
    }
}
signed main(){
    scanf("%d",&n);
    for(int i=0;i<n;i++) scanf("%d",&a[i]);
    for(len=1;len<=n;len<<=1);
    inv(a,b,len);
    for(int i=0;i<n;i++) printf("%d ",b[i]);
    return 0;
}

多项式求逆

给出一个 \(n-1\) 次多项式 \(A(x)\)，求 \(B(x)\) 满足 \(A(x)B(x)\equiv 1\pmod {x^n}\)。

注意，这里的模 \(x^n\) 表示只考虑多项式的 \(x^0,x^1,\cdots,x^{n-1}\) 项（因为 \(x^n\) 及更高次模 \(x^n\) 就为 \(0\) 了）。

为了方便，设 \(n\) 为 \(2\) 的幂次。

考虑倍增。假设我们已知 \(A(x)B_0(x)\equiv 1\pmod{x^{n/2}}\)。因为 \(A(x)B(x)\equiv 1\pmod{x^n}\)，所以也有 \(A(x)B(x)\equiv 1\pmod{x^{n/2}}\)。

相减：\(A(x)(B(x)-B_0(x))\equiv 0\pmod{x^{n/2}}\)。

因为 \(A(x)\) 有逆元，故 \(B(x)-B_0(x)\equiv 0\pmod{x^{n/2}}\)。两边平方得 \(B(x)^2-2B(x)B_0(x)+B_0(x)^2\equiv 0\pmod {x^n}\)，即 \(B(x)^2\equiv 2B(x)B_0(x)-B_0(x)^2\pmod{x^n}\)。

两边同乘 \(A(x)\)，用 \(A(x)B(x)\equiv 1\pmod{x^n}\) 化简：

\[B(x)\equiv 2B_0(x)-A(x)B_0(x)^2\pmod{x^n} \]

边界为 \(n=1\) 时 \(B(x)=A(0)^{-1}\)。时间复杂度 \(T(n)=T(\frac n 2)+\mathcal O(n\log n)=\mathcal O(n\log n)\)。

//Luogu P4238
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,a[N],b[N],len,r[N],A[N],B[N],p[N],q[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void mul(int *a,int *b,int n){
	for(int i=0;i<n;i++)
		p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?n>>1:0);
	NTT(p,n,1),NTT(q,n,1);
	for(int i=0;i<n;i++) a[i]=1ll*p[i]*q[i]%mod;
	NTT(a,n,-1);
}
void inv(int *a,int *b,int n){
	for(int i=0;i<n;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=n;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		//A(x)B0(x)^2 长度最大为 i*2
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
signed main(){
	scanf("%d",&n);
	for(int i=0;i<n;i++) scanf("%d",&a[i]);
	for(len=1;len<=n;len<<=1);
	inv(a,b,len);
	for(int i=0;i<n;i++) printf("%d ",b[i]);
	return 0;
}

多项式开根

给出一个 \(n-1\) 次多项式 \(A(x)\)，求 \(B(x)\) 满足 \(B^2(x)\equiv A(x)\pmod {x^n}\)。若有多解，取 \(x^0\) 系数较小的作为答案。

保证 \(a_0=1\)。

倍增，假设已知 \(B_0^2(x)\equiv A(x)\pmod{x^{n/2}}\)。因为 \(B^2(x)\equiv A(x)\pmod{x^n}\)，所以也有 \(B^2(x)\equiv A(x)\pmod{x^{n/2}}\)。

\([B(x)-B_0(x)][B(x)+B_0(x)]\equiv 0\pmod{x^{n/2}}\)，发现这个式子中 \(B(x)\) 有两个解，取其中一个解得到 \(B(x)-B_0(x)\equiv 0\pmod{x^{n/2}}\)。平方，\(B^2(x)-2B(x)B_0(x)+B_0^2(x)\equiv 0\pmod {x^n}\)。

又由于 \(B^2(x)\equiv A(x)\pmod{x^n}\)，故 \(A(x)-2B(x)B_0(x)+B_0^2(x)\equiv 0\pmod{x^n}\)，即：

\[B(x)\equiv \frac{A(x)+B_0^2(x)}{2B_0(x)}\pmod{x^n} \]

边界为 \(n=1\) 时 \(G(x)=1\)。写个卷积、求逆求解，时间复杂度 \(\mathcal O(n\log n)\)。

//Luogu P5205
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,a[N],b[N],len,r[N],A[N],B[N],p[N],q[N],f[N],g[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void mul(int *a,int *b,int n){
	for(int i=0;i<n;i++)
		p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?n>>1:0);
	NTT(p,n,1),NTT(q,n,1);
	for(int i=0;i<n;i++) a[i]=1ll*p[i]*q[i]%mod;
	NTT(a,n,-1);
}
void inv(int *a,int *b,int n){
	for(int i=0;i<n;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=n;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
void sqrt(int *a,int *b,int n){
	for(int i=0;i<n;i++) b[i]=0;
	b[0]=1;
	for(int i=2;i<=n;i<<=1){
		for(int j=0;j<(i<<1);j++) f[j]=b[j],g[j]=0;
		mul(f,f,i);
		for(int j=0;j<i;j++) f[j]=(f[j]+a[j])%mod;
		inv(b,g,i),mul(f,g,i<<1);
		for(int j=0;j<i;j++) b[j]=1ll*f[j]*((mod+1)/2)%mod;
	}
}
signed main(){
	scanf("%d",&n);
	for(int i=0;i<n;i++) scanf("%d",&a[i]);
	for(len=1;len<=n;len<<=1);
	sqrt(a,b,len);
	for(int i=0;i<n;i++) printf("%d ",b[i]);
	return 0;
}

多项式求导/积分

主要是为了写一点前置知识。

求导：\((x^n)'=nx^{n-1}\)，\((c)'=0\)（\(c\) 为常数），\((\ln x)'=\large\frac 1 x\)，\((e^x)'=e^x\)。

\((fg)'=f'g+fg'\)，\((\frac{f}{g})'=\large\frac{f'g-fg'}{g^2}\)，链式法则 \([f(g(x))]'=f'(g(x))g'(x)\)。

  1. 若多项式 \(f(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n\)，则 \(f'(x)=a_1+2a_2x+3a_3x^2+\cdots+na_nx^{n-1}\)。
  2. 若 \(f(x)=\ln x\)，则 \(f'(x)=\large\frac 1 x\)；若 \(f(x)=e^x\)，则 \(f'(x)=e^x\)。

积分：\(\int f(x)\text dx=F(x)\)，\(f(x)\) 被称为“被积函数”，\(F(x)\) 被称为“ 原函数”。积分是导数的逆运算。

\(F(x)\) 是 \(f(x)\) 的一个原函数，\(F'(x)=f(x),\int f(x)\text dx=F(x)+C\)。

  1. 若多项式 \(f(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n\)，则 \(\int f(x)=C+a_0x+\large\frac{a_1}{2}\normalsize x^2+\large\frac{a_2}{3}\normalsize x^3\cdots+\large\frac{a_n}{n+1}\normalsize x^{n+1}\)，其实就是求导那个式子反着推。
  2. 若 \(f(x)=e^x\)，则 \(\int f(x)=e^x+C\)。
void diff(int *a,int *b,int n){
	for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%mod; b[n-1]=0;
} 
void integ(int *a,int *b,int n){
	for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*qpow(i,mod-2)%mod; b[0]=0;
}

多项式 ln

给出一个 \(n-1\) 次多项式 \(A(x)\)，求 \(B(x)\) 满足 \(B(x)\equiv \ln A(x)\pmod{x^n}\)。

保证 \(a_0=1\)。

根据链式法则，\([\ln(A(x))]'=\ln'(A(x))A'(x)=\large\frac{1}{A(x)}\normalsize A'(x)\)。故两边同时求导可得 \(B'(x)\equiv \large\frac{A'(x)}{A(x)}\normalsize\pmod{x^n}\)。

后面那东西可以多项式求逆 + 卷积求出，然后再对 \(B'(x)\) 积分就能解出 \(B(x)\)，\(B(x)=\int \large\frac{A'(x)}{A(x)}\)。

//Luogu P4725
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,a[N],b[N],len,r[N],A[N],B[N],p[N],q[N],f[N],g[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void mul(int *a,int *b,int n){
	for(int i=0;i<n;i++)
		p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?n>>1:0);
	NTT(p,n,1),NTT(q,n,1);
	for(int i=0;i<n;i++) a[i]=1ll*p[i]*q[i]%mod;
	NTT(a,n,-1);
}
void inv(int *a,int *b,int n){
	for(int i=0;i<n;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=n;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
void diff(int *a,int *b,int n){
	for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%mod; b[n-1]=0;
} 
void integ(int *a,int *b,int n){
	for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*qpow(i,mod-2)%mod; b[0]=0;
}
void ln(int *a,int *b,int n){
	for(int i=0;i<(n<<1);i++) f[i]=g[i]=0;
	diff(a,f,n),inv(a,g,n),mul(f,g,n<<1),integ(f,b,n);
}
signed main(){
	scanf("%d",&n);
	for(int i=0;i<n;i++) scanf("%d",&a[i]);
	for(len=1;len<=n;len<<=1);
	ln(a,b,len);
	for(int i=0;i<n;i++) printf("%d ",b[i]);
	return 0;
}

多项式除法/取模

给出一个 \(n\) 次多项式 \(A(x)\) 和一个 \(m\) 次多项式 \(B(x)\)，求一个 \(n-m\) 次多项式 \(C(x)\) 和 \(m-1\) 次多项式 \(D(x)\)，满足 \(A(x)=B(x)C(x)+D(x)\)。

对于一个 \(n\) 次多项式 \(A(x)\)，令 \(A^R(x)=x^nA(\large\frac 1 x\normalsize)\)。

发现，\(A(x)\) 里的每一项 \(a_ix^i\)，到 \(A^R(x)\) 里都会变成 \(a_ix^{n-i}\)，也就是说，\(A^R(x)\) 就是将 \(A(x)\) 的系数反转。

将 \(x\) 替换为 \(\large\frac 1 x\)，再同时乘 \(x^n\) 得到：

\[\begin{aligned} x^nA(\frac 1 x)&=x^m B(\frac 1 x)\times x^{n-m}C(\frac 1 x)+x^{n-m+1}\times x^{m-1}D(\frac 1x)\\ A^R(x)&=B^R(x)C^R(x)+x^{n-m+1}D^R(x)\\ A^R(x)&\equiv B^R(x)C^R(x)\pmod{x^{n-m+1}} \end{aligned} \]

用多项式求逆即可解出 \(C^R(x)\)，也就解出了 \(C(x)\)。最后把 \(C(x)\) 代回去就可以解出 \(D(x)\)。时间复杂度 \(\mathcal O(n\log n)\)。

//Luogu P4512
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,m,a[N],b[N],c[N],len,r[N],A[N],B[N],p[N],q[N],f[N],g[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void mul(int *a,int *b,int n){
	for(int i=0;i<n;i++)
		p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?n>>1:0);
	NTT(p,n,1),NTT(q,n,1);
	for(int i=0;i<n;i++) a[i]=1ll*p[i]*q[i]%mod;
	NTT(a,n,-1);
}
void inv(int *a,int *b,int n){
	for(int i=0;i<n;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=n;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%d",&a[i]);
	for(int i=0;i<=m;i++) scanf("%d",&b[i]);
	reverse(a,a+1+n),reverse(b,b+1+m);
	for(len=1;len<=n;len<<=1);
	inv(b,c,len),mul(c,a,len<<1),reverse(c,c+1+n-m);
	for(int i=n-m+1;i<(len<<1);i++) c[i]=0;
	for(int i=0;i<=n-m;i++) printf("%d ",c[i]); puts("");
	reverse(a,a+1+n),reverse(b,b+1+m),mul(c,b,len<<1);
	for(int i=0;i<m;i++) printf("%d ",(a[i]-c[i]+mod)%mod);
	return 0;
}

vector：效率 & 空间 相对较劣。

#include<bits/stdc++.h>
#define poly vector<int>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,m,len,r[N];
poly a,b,c;
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(poly &a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
poly operator*(poly a,poly b){
	for(len=1;len<(int)a.size()+(int)b.size();len<<=1);
	for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|(i&1?len>>1:0);
	a.resize(len),b.resize(len);
	NTT(a,len,1),NTT(b,len,1);
	for(int i=0;i<len;i++) a[i]=1ll*a[i]*b[i]%mod;
	return NTT(a,len,-1),a;
}
poly inv(poly a,int n){
	poly b(n,0);
	a.resize(n),b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=n;i<<=1){
		poly c(a.begin(),a.begin()+i),d(b.begin(),b.begin()+(i>>1));
		c=d*d*c;
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-c[j]+mod)%mod;
	}
	return b;
}
poly operator/(poly a,poly b){
	int n=a.size()-b.size()+1;
	reverse(a.begin(),a.end());
	reverse(b.begin(),b.end());
	for(len=1;len<=n;len<<=1);
	b=inv(b,len),b.resize(n);
	a=a*b,a.resize(n);
	return reverse(a.begin(),a.end()),a;
}
poly operator%(poly a,poly b){
	poly c=(a/b)*b;
	c.resize(b.size()-1);
	for(int i=0;i<(int)b.size()-1;i++) c[i]=(a[i]-c[i]+mod)%mod;
	return c;
}
signed main(){
	scanf("%d%d",&n,&m),a.resize(n+1),b.resize(m+1);
	for(int i=0;i<=n;i++) scanf("%d",&a[i]);
	for(int i=0;i<=m;i++) scanf("%d",&b[i]);
	c=a/b;
	for(int i:c) printf("%d ",i); puts("");
	c=a%b;
	for(int i:c) printf("%d ",i);
	return 0;
}

多项式 exp

泰勒展开

对于函数 \(f\)，设它在点 \(x_0\) 存在直到 \(n\) 阶的导数，则我们由这些导数构造一个 \(n\) 次多项式：

\[g(x)=f(x_0)+\frac{f'(x_0)}{1!}(x-x_0)+\frac{f''(x_0)}{2!}{(x-x_0)}^2+\cdots+\frac{f^{(n)}(x_0)}{n!}{(x-x_0)}^n \]

称为 \(f\) 在点 \(x_0\) 处的泰勒多项式。

主要思想是用多项式来逼近一个函数。也就是将一个不是形式幂级数的形式（如 \(\cos(x)\)），表示成形式幂级数，即 \(g(x)\sum_i a_ix^i\) 的形式。

以 \(x_0=0\) 为例，我们希望 \(g(x)\) 与 \(f(x)\) 尽可能重合，所以令 \(g(0)=f(0)\)。在此基础上，考虑曲线的变化趋势，即导数，保证在此处的导数相同，故我们再令 \(g'(0)=f'(0)\)。然后再令 \(g''(0)=f''(0),g^{(3)}(0)=f^{(3)}(0)\)，以此类推，来不断逼近函数。

显然求完 \(n\) 阶导，再代入 \(0\) 后就只剩下常数项了，而求完 \(n\) 阶导后的常数项实际上是 \(n!a_n\)，故我们有 \(a_n=\large\frac{f^{(n)}(0)}{n!}\)，于是 \(g(x)=f(0)+\dfrac{f^{(1)}(0)}{1!}x+\dfrac{f^{(2)}(0)}{2!}x^2+\dfrac{f^{(3)}}{3!}x^3+\cdots+\dfrac{f^{(r)}(0)}{r!}x^r+\cdots\)。

在上述栗子中我们是以 \(0\) 为关键点的，如果将 \(0\) 换为 \(x_0\)，重复上述步骤也能得到类似的式子。

\[g(x)=f(x_0)+\frac{f^{(1)}(x_0)}{1!}(x-x_0)+\frac{f^{(2)}(x_0)}{2!}{(x-x_0)}^2+\cdots+\frac{f^{(r)}(x_0)}{r!}{(x-x_0)}^r+\cdots \]

当然有时候我们不可能无止尽地递归下去，譬如我们只觉得只求 \(n\) 阶导就足够了。

\[g(x)=f(x_0)+\frac{f^{(1)}(x_0)}{1!}(x-x_0)+\frac{f^{(2)}(x_0)}{2!}{(x-x_0)}^2+\cdots+\frac{f^{(n)}(x_0)}{n!}{(x-x_0)}^n+r_n(x) \]

其中 \(r_n(x)\) 被称为泰勒公式的拉格朗日余项。

当 \(x_0=0\) 时，上述式子又被称为 \(f(x)\) 的麦克劳林展开式。

牛顿迭代

已知一个函数 \(g(x)\)，求一个多项式 \(f(x)\bmod x^n\)，满足 \(g(f(x))\equiv 0\pmod{x^n}\)。

考虑倍增，假设已知 \(g(f_0(x))\equiv 0\pmod{x^{n/2}}\)。

将 \(g(f(x))\) 在 \(f_0(x)\) 处泰勒展开，得到：

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

由于 \(g(f(x))\equiv 0\pmod{x^{n/2}}\)，故 \(f(x)\equiv f_0(x)\pmod{x^{n/2}}\)。所以 \(f(x)-f_0(x)\) 的最低次项至少为 \(x^{n/2}\)（因为 \(f(x),f_0(x)\) 的 \(x^0,x^1,\cdots,x^{n/2-1}\) 都是一样的，相减就是 \(0\)），则 \((f(x)-f_0(x))^2\) 最低次项至少为 \(x^n\)。因此在 \(\bmod x^n\) 意义下，整个式子从 \((f(x)-f_0(x))^2\) 这一项开始都为 \(0\) 了。即：

\[g(f(x))\equiv g(f_0(x))+g'(f_0(x))(f(x)-f_0(x))\pmod{x^n} \]

又因为 \(g(f(x))\equiv 0\pmod{x^n}\)，得到：

\[\begin{aligned} g'(f_0(x))f(x)&\equiv g'(f_0(x))f_0(x)-g(f_0(x))\pmod{x^n}\\ f(x)&\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\pmod{x^n} \end{aligned} \]

正题

给出一个 \(n-1\) 次多项式 \(A(x)\)，求 \(B(x)\) 满足 \(B(x)\equiv e^{A(x)}{\pmod{x^n}}\)。

保证 \(a_0=0\)。

两边同时取 \(\ln\) 得到 \(\ln B(x)\equiv A(x)\pmod{x^n}\)。

设 \(g(B(x))=\ln B(x)-A(x)\)，那么我们的目标是求满足 \(g(B(x))\equiv 0\pmod{x^n}\) 的 \(B(x)\)。

假设已知 \(g(B_0(x))\equiv 0\pmod{x^{n/2}}\)。牛顿迭代得 \(B(x)\equiv B_0(x)-\dfrac{g(B_0(x))}{g'(B_0(x))}\)。即：

\[\begin{aligned} B(x)&\equiv B_0(x)-\frac{\ln B_0(x)-A(x)}{\frac{1}{B_0(x)}}\pmod{x^n}\\ B(x)&\equiv B_0(x)(1-\ln B_0(x)+A(x))\pmod{x^n} \end{aligned} \]

//Luogu P4726
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,a[N],b[N],c[N],len,r[N],A[N],B[N],p[N],q[N],f[N],g[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void mul(int *a,int *b,int n){
	for(int i=0;i<n;i++)
		p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?n>>1:0);
	NTT(p,n,1),NTT(q,n,1);
	for(int i=0;i<n;i++) a[i]=1ll*p[i]*q[i]%mod;
	NTT(a,n,-1);
}
void inv(int *a,int *b,int n){
	for(int i=0;i<n;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=n;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
void diff(int *a,int *b,int n){
	for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%mod; b[n-1]=0;
} 
void integ(int *a,int *b,int n){
	for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*qpow(i,mod-2)%mod; b[0]=0;
}
void ln(int *a,int *b,int n){
	for(int i=0;i<(n<<1);i++) f[i]=g[i]=0;
	diff(a,f,n),inv(a,g,n),mul(f,g,n<<1),integ(f,b,n);
}
void exp(int *a,int *b,int n){
	for(int i=0;i<(n<<1);i++) b[i]=c[i]=0;
	b[0]=1;
	for(int i=2;i<=n;i<<=1){
		ln(b,c,i);
		for(int j=0;j<i;j++) c[j]=(-c[j]+a[j]+mod)%mod;
		c[0]=(c[0]+1)%mod,mul(b,c,i<<1);
	}
}
signed main(){
	scanf("%d",&n);
	for(int i=0;i<n;i++) scanf("%d",&a[i]);
	for(len=1;len<=n;len<<=1);
	exp(a,b,len);
	for(int i=0;i<n;i++) printf("%d ",b[i]);
	return 0;
}

vector 版，代码更短空间更小，效率相差不大：

#include<bits/stdc++.h>
#define poly vector<int>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,len,r[N];
poly a;
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(poly &a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
poly operator*(poly a,poly b){
	for(len=1;len<(int)a.size()+(int)b.size();len<<=1);
	for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|(i&1?len>>1:0);
	a.resize(len),b.resize(len);
	NTT(a,len,1),NTT(b,len,1);
	for(int i=0;i<len;i++) a[i]=1ll*a[i]*b[i]%mod;
	return NTT(a,len,-1),a;
}
poly inv(poly a,int n){
	poly b(n,0);
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=n;i<<=1){
		poly c(a.begin(),a.begin()+i),d(b.begin(),b.begin()+(i>>1));
		c=d*d*c;
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-c[j]+mod)%mod;
	}
	return b;
}
poly diff(poly a,int n){
	poly b(n,0);
	for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%mod;
	return b;
}
poly integ(poly a,int n){
	poly b(n,0);
	for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*qpow(i,mod-2)%mod;
	return b;
}
poly ln(poly a,int n){
	return integ(diff(a,n)*inv(a,n),n);
}
poly exp(poly a,int n){
	poly b(1,1),c;
	a.resize(n);
	for(int i=2;i<=n;i<<=1){
		c=ln(b,i);
		for(int j=0;j<i;j++) c[j]=(-c[j]+a[j]+mod)%mod;
		c[0]=(c[0]+1)%mod,b=b*c,b.resize(i<<1);
	}
	return b;
}
signed main(){
	scanf("%d",&n),a.resize(n);
	for(int i=0;i<n;i++) scanf("%d",&a[i]);
	for(len=1;len<=n;len<<=1);
	a=exp(a,len);
	for(int i=0;i<n;i++) printf("%d ",a[i]);
	return 0;
}

多项式快速幂

给出一个 \(n-1\) 次多项式 \(A(x)\)，求 \(B(x)\) 满足 \(B(x)\equiv A^k(x)\pmod{x^n}\)。

保证 \(a_0=1\)。

直接快速幂 + NTT，时间复杂度是 \(\mathcal O(n\log n\log k)\) 的，而且常数还比较大。

显然 \(A^k(x)=(e^{\ln A(x)})^k=e^{k\times \ln A(x)}\)。所以我们用多项式 ln,exp 就可以了。复杂度 \(\mathcal O(n\log n)\)。

//Luogu P5245
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,m,k,a[N],b[N],c[N],len,r[N],A[N],B[N],p[N],q[N],f[N],g[N],tmp[N];
char s[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void mul(int *a,int *b,int n){
	for(int i=0;i<n;i++)
		p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?n>>1:0);
	NTT(p,n,1),NTT(q,n,1);
	for(int i=0;i<n;i++) a[i]=1ll*p[i]*q[i]%mod;
	NTT(a,n,-1);
}
void inv(int *a,int *b,int n){
	for(int i=0;i<n;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=n;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
void diff(int *a,int *b,int n){
	for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%mod; b[n-1]=0;
} 
void integ(int *a,int *b,int n){
	for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*qpow(i,mod-2)%mod; b[0]=0;
}
void ln(int *a,int *b,int n){
	for(int i=0;i<(n<<1);i++) f[i]=g[i]=0;
	diff(a,f,n),inv(a,g,n),mul(f,g,n<<1),integ(f,b,n);
}
void exp(int *a,int *b,int n){
	for(int i=0;i<(n<<1);i++) b[i]=c[i]=0;
	b[0]=1;
	for(int i=2;i<=n;i<<=1){
		ln(b,c,i);
		for(int j=0;j<i;j++) c[j]=(-c[j]+a[j]+mod)%mod;
		c[0]=(c[0]+1)%mod,mul(b,c,i<<1);
	}
}
signed main(){
	scanf("%d%s",&n,s+1),m=strlen(s+1);
	for(int i=1;i<=m;i++) k=(10ll*k+s[i]-'0')%mod;
	for(int i=0;i<n;i++) scanf("%d",&a[i]);
	for(len=1;len<=n;len<<=1);
	ln(a,tmp,len);
	for(int i=0;i<len;i++) tmp[i]=1ll*tmp[i]*k%mod;
	exp(tmp,b,len);
	for(int i=0;i<n;i++) printf("%d ",b[i]);
	return 0;
}

vector 版：

#include<bits/stdc++.h>
#define poly vector<int>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,k,len,r[N];
char s[N];
poly a;
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(poly &a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
poly operator*(poly a,poly b){
	for(len=1;len<(int)a.size()+(int)b.size();len<<=1);
	for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|(i&1?len>>1:0);
	a.resize(len),b.resize(len);
	NTT(a,len,1),NTT(b,len,1);
	for(int i=0;i<len;i++) a[i]=1ll*a[i]*b[i]%mod;
	return NTT(a,len,-1),a;
}
poly inv(poly a,int n){
	poly b(n,0);
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=n;i<<=1){
		poly c(a.begin(),a.begin()+i),d(b.begin(),b.begin()+(i>>1));
		c=d*d*c;
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-c[j]+mod)%mod;
	}
	return b;
}
poly diff(poly a,int n){
	poly b(n,0);
	for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%mod;
	return b;
}
poly integ(poly a,int n){
	poly b(n,0);
	for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*qpow(i,mod-2)%mod;
	return b;
}
poly ln(poly a,int n){
	return integ(diff(a,n)*inv(a,n),n);
}
poly exp(poly a,int n){
	poly b(1,1),c;
	a.resize(n);
	for(int i=2;i<=n;i<<=1){
		c=ln(b,i);
		for(int j=0;j<i;j++) c[j]=(-c[j]+a[j]+mod)%mod;
		c[0]=(c[0]+1)%mod,b=b*c,b.resize(i<<1);
	}
	return b;
}
poly Pow(poly a,int k,int n){
	a.resize(n);
	poly b=ln(a,n);
	for(int i=0;i<n;i++) b[i]=1ll*b[i]*k%mod;
	return exp(b,n);
}
signed main(){
	scanf("%d%s",&n,s+1),a.resize(n);
	int m=strlen(s+1);
	for(int i=1;i<=m;i++) k=(10ll*k+s[i]-'0')%mod;
	for(int i=0;i<n;i++) scanf("%d",&a[i]);
	for(len=1;len<=n;len<<=1);
	a=Pow(a,k,len);
	for(int i=0;i<n;i++) printf("%d ",a[i]);
	return 0;
}

如果 \(a_0\neq 1\)：（以下代码适用于 \(k\) 小的情况）

poly Pow_(poly a,int k,int n){
	poly b=ln(a,n);
	for(int i=0;i<n;i++) b[i]=1ll*b[i]*k%mod;
	return exp(b,n);
}
poly Pow(poly a,int k,int n){
	a.resize(n);
	int pos=0;
	for(int i=0;i<n;i++) if(a[i]){pos=i;break;}
	poly tmp(n,0),b(n,0);
	for(int i=pos;i<n;i++) tmp[i-pos]=a[i];
	int x=tmp[0],ivx=qpow(x,mod-2);
	for(int i=0;i<n;i++) tmp[i]=1ll*tmp[i]*ivx%mod;
	tmp=Pow_(tmp,k,n),x=qpow(x,k);
	for(int i=0;i+pos*k<n;i++)
		b[i+pos*k]=1ll*tmp[i]*x%mod;
	return b.resize(n),b;
}

如果 \(k\) 很大：

  1. 记录 k 表示 \(k\bmod 998244353\)，用作 Pow_ 中的 b[i]=1ll*b[i]*k%mod 的 \(k\)。
  2. 记录 k2 表示 \(k\bmod (998244353-1)\)，用作费马小定理算 qpow(x,k)。
  3. 记录 k3 表示 \(\min(k,n)\)，用作 i+1ll*pos*k<n。注意要开 1ll。
//O2
#include<bits/stdc++.h>
#define poly vector<int>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,k,k2,k3,len,r[N];
char s[N];
poly a;
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(poly &a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
poly operator*(poly a,poly b){
	for(len=1;len<(int)a.size()+(int)b.size();len<<=1);
	for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|(i&1?len>>1:0);
	a.resize(len),b.resize(len);
	NTT(a,len,1),NTT(b,len,1);
	for(int i=0;i<len;i++) a[i]=1ll*a[i]*b[i]%mod;
	return NTT(a,len,-1),a;
}
poly inv(poly a,int n){
	poly b(n,0);
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=n;i<<=1){
		poly c(a.begin(),a.begin()+i),d(b.begin(),b.begin()+(i>>1));
		c=d*d*c;
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-c[j]+mod)%mod;
	}
	return b;
}
poly diff(poly a,int n){
	poly b(n,0);
	for(int i=1;i<n;i++) b[i-1]=1ll*a[i]*i%mod;
	return b;
}
poly integ(poly a,int n){
	poly b(n,0);
	for(int i=1;i<n;i++) b[i]=1ll*a[i-1]*qpow(i,mod-2)%mod;
	return b;
}
poly ln(poly a,int n){
	return integ(diff(a,n)*inv(a,n),n);
}
poly exp(poly a,int n){
	poly b(1,1),c;
	a.resize(n);
	for(int i=2;i<=n;i<<=1){
		c=ln(b,i);
		for(int j=0;j<i;j++) c[j]=(-c[j]+a[j]+mod)%mod;
		c[0]=(c[0]+1)%mod,b=b*c,b.resize(i<<1);
	}
	return b;
}
poly Pow_(poly a,int k,int n){
	poly b=ln(a,n);
	for(int i=0;i<n;i++) b[i]=1ll*b[i]*k%mod;
	return exp(b,n);
}
poly Pow(poly a,int k,int n){
	a.resize(n);
	int pos=0;
	for(int i=0;i<n;i++) if(a[i]){pos=i;break;}
	poly tmp(n,0),b(n,0);
	for(int i=pos;i<n;i++) tmp[i-pos]=a[i];
	int x=tmp[0],ivx=qpow(x,mod-2);
	for(int i=0;i<n;i++) tmp[i]=1ll*tmp[i]*ivx%mod;
	tmp=Pow_(tmp,k,n),x=qpow(x,k2);
	for(int i=0;i+1ll*pos*k3<n;i++)
		b[i+pos*k3]=1ll*tmp[i]*x%mod;
	return b.resize(n),b;
}
signed main(){
	scanf("%d%s",&n,s+1),a.resize(n);
	int m=strlen(s+1);
	for(int i=1;i<=m;i++){
		k=(10ll*k+s[i]-'0')%mod,k2=(10ll*k2+s[i]-'0')%(mod-1),
		k3=min(0ll+n,10ll*k3+s[i]-'0');
	}
	for(int i=0;i<n;i++) scanf("%d",&a[i]);
	for(len=1;len<=n;len<<=1);
	a=Pow(a,k,len);
	for(int i=0;i<n;i++) printf("%d ",a[i]);
	return 0;
}

这个板子非常慢，在 ZR#2618 被卡爆了。yny 模板。

各种 n^2 递推

更多

多项式求逆

设 \(B(x)\) 为 \(A(x)\) 的逆元，\(A(x)B(x)=1\)，\(\sum_{i=0}^n a_ib_{n-i}=[n=0]\)。

如果 \(n\geq 1\)，\(a_0b_n=-\sum_{i=1}^n a_ib_{n-i}\)，\(b_n=-\large\frac{1}{a_0}\normalsize\sum_{i=1}^n a_ib_{n-i}\)。

边界为 \(b_0=\large\frac{1}{a_0}\)。

多项式求 ln

设 \(B(x)=\ln A(x)\)，两边同时求导得 \(B'(x)=\large\frac{A'(x)}{A(x)}\)。然后就能 \(n^2\) 求逆，再对 \(B'(x)\) 积分解出 \(B(x)\) 了，\(B(x)=\int \large\frac{A'(x)}{A(x)}\)。

upd：不需要求逆。

\(A(x)B'(x)=A'(x)\)

\([x^{i-1}]A'(x)=ia_i\) 不方便，考虑转成 \([x^i]xA'(x)=ia_i\)

\(A(x)xB'(x)=xA'(x)\)

\(\sum_{i=0}^n ib_ia_{n-i}=na_n\)

\(nb_na_0=na_n-\sum_{i=0}^{n-1}ib_ia_{n-i}\)

由于 \(a_0=1\)，\(b_n=a_n-\frac 1 n\sum_{i=0}^{n-1}ib_ia_{n-i}\)。边界为 \(b_0=0\)。

void ln(){
	b[0]=0;
	for(int i=1;i<=n;i++){
		b[i]=0;
		for(int j=0;j<i;j++) b[i]=(b[i]+1ll*j*b[j]%mod*a[i-j])%mod;
		b[i]=(a[i]-1ll*b[i]*inv[i])%mod;
	}
}

多项式 exp

\([e^{A(x)}]'=e^{A(x)}A'(x)\)，设 \(B(x)=e^{A(x)}\)，两边求导得 \(B'(x)=B(x)A'(x)\)，\(b'_n=\sum_{i=0}^n b_ia'_{n-i}\)。

\(b_{n+1}(n+1)=\sum_{i=0}^n b_ia_{n-i+1}(n-i+1)\)，\(b_{n+1}=\large\frac{1}{n+1}\normalsize \sum_{i=0}^n b_ia_{n-i+1}(n-i+1)\)。

同样的，做 exp 要求 \(a_0=0\)，所以边界为 \(b_0=1\)。

upd：可以和 \(\ln\) 统一起来。

\(A(x)=e^{B(x)}\)，就是 \(B(x)=\ln A(x)\)，把上面求 \(\ln\) 的式子拎下来即可。

\(a_n=\frac 1 n\sum_{i=0}^n ib_ia_{n-i}\)

void exp(){
	b[0]=1;
	for(int i=1;i<=n;i++){
		b[i]=0;
		for(int j=0;j<=i;j++) b[i]=(b[i]+1ll*j*a[j]%mod*b[i-j])%mod;
		b[i]=1ll*b[i]*inv[i]%mod;
	}
}

例题

CF438E The Child and Binary Tree（*3100）

对所有 \(i\in[1,m]\)，求所有点权在集合 \(S\) 中，且点权之和为 \(i\) 的有根二叉树个数。对 \(998244353\) 取模。

\(1\leq n,m\leq 10^5\)，\(\forall x\in S,1\leq x\leq 10^5\)。

设 \(f_i\) 表示点权之和为 \(i\) 的符合要求的二叉树个数。记 \(g_i\) 表示 \(i\) 是否在 \(S\) 中出现过，则 \(f_i=\sum_{x+y+z=S}g_xf_yf_z\)，边界为 \(f_0=1\)。

设 \(F,G\) 分别为 \(f_i,g_i\) 的生成函数，那么 \(F=F^2G+1\)。

移项，\(F^2G-F+1=0\)。根据求根公式得，\(F=\large\frac{1\pm \sqrt{1-4G}}{2G}\)。

也就是 \(2FG=1\pm \sqrt{1-4G}\)，由于 \(g_0=0\)，\(2FG\) 的常数项也是 \(0\)，故排除正根。所以 \(F=\large\frac{1-\sqrt{1-4G}}{2G}\)。

而 \(g_0=0\) 没有逆元，不能直接多项式求逆，分母有理化得 \(F=\large\frac{2}{1+\sqrt{1-4G}}\)，多项式开根、多项式求逆即可。

#include<bits/stdc++.h>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,m,x,a[N],b[N],c[N],len,r[N],A[N],B[N],p[N],q[N],f[N],g[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void mul(int *a,int *b,int len){
	for(int i=0;i<len;i++)
		p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?len>>1:0);
	NTT(p,len,1),NTT(q,len,1);
	for(int i=0;i<len;i++) a[i]=1ll*p[i]*q[i]%mod;
	NTT(a,len,-1);
}
void inv(int *a,int *b,int len){
	for(int i=0;i<len;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=len;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
void sqrt(int *a,int *b,int len){
	for(int i=0;i<len;i++) b[i]=0;
	b[0]=1;
	for(int i=2;i<=len;i<<=1){
		for(int j=0;j<(i<<1);j++) f[j]=b[j],g[j]=0;
		mul(f,f,i);
		for(int j=0;j<i;j++) f[j]=(f[j]+a[j])%mod;
		inv(b,g,i),mul(f,g,i<<1);
		for(int j=0;j<i;j++) b[j]=1ll*f[j]*((mod+1)/2)%mod;
	}
}
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;i++) scanf("%d",&x),a[x]=1;
	for(int i=1;i<=m;i++)
		a[i]=(mod-4ll*a[i]%mod)%mod; a[0]++;
	for(len=1;len<=m;len<<=1);
	sqrt(a,b,len),b[0]++,inv(b,c,len);
	for(int i=1;i<=m;i++) printf("%lld\n",2ll*c[i]%mod);
	return 0;
}

P4389 付公主的背包

有 \(n\) 个物品，每个物品有无限多，第 \(i\) 个物品的体积为 \(v_i\)。

对于所有 \(i\in[1,m]\)，求用这些物品恰好装满容量为 \(i\) 的背包的方案数，对 \(998244353\) 取模。

\(1\leq n,m\leq 10^5\)，\(1\leq v_i\leq m\)。

设 \(f_i\) 表示选择的物品体积为 \(i\) 的方案数。

显然 \(F=(1+x^{v_1}+x^{2v_1}+\cdots)(1+x^{v_2}+x^{2v_2}+\cdots)\cdots(1+x^{v_n}+x^{2v_n}+\cdots)\)。

等比数列求和，\(1+x^v+x^{2v}+\cdots=\large\frac{1}{1-x^v}\)。故 \(F=\prod_{i=1}^n\large\frac{1}{1-x^{v_i}}\)。现在我们的目标就是求出 \((1-x^{v_1})(1-x^{v_2})\cdots(1-x^{v_n})\)。

根据泰勒展开，\(\ln(1+x)=x-\large\frac{x^2}{2}\normalsize+\large\frac{x^3}{3}\normalsize-\large\frac{x^4}{4}\normalsize+\cdots\)。

设 \(T=\ln((1-x^{v_1})(1-x^{v_2})\cdots(1-x^{v_n}))\)。易知 \(T=\ln(1-x^{v_1})+\ln(1-x^{v_2})+\cdots+\ln(1-x^{v_n})\)。

将 \(-x^v\) 代入 \(\ln(1+x)\) 可得 \(\ln(1-x^v)=-x^v-\large\frac{x^{2v}}{2}\normalsize-\large\frac{x^{3v}}{3}\normalsize-\large\frac{x^{4v}}{4}\normalsize-\cdots\)，由于我们要求 \(F\bmod x^{m+1}\)，只需枚举到 \(x^{\large\lfloor\frac{m+1}{v}\rfloor}\)。

这样多项式乘法就变成了多项式加法。如果 \(v_i\) 都很小，求一遍 \(\ln(1-x^v)\) 是 \(\mathcal O(m)\) 的怎么办？开个桶 \(c_i\) 记录 \(i\) 在 \(v\) 数组出现了多少次，然后将 \(c_i\ln(1-x^i)\) 累加进 \(T\)。这样复杂度就是 \(\large\frac m 1\normalsize+\large\frac m 2\normalsize+\cdots+\large\frac m m\normalsize=\mathcal O(m\ln m)\) 了。

求出 \(T\) 后做一遍 exp，再求个逆就能得到 \(F\)。

#include<bits/stdc++.h>
using namespace std;
const int N=(1<<18)+5,mod=998244353;
int n,m,x,cnt[N],a[N],b[N],c[N],len,r[N],A[N],B[N],p[N],q[N],f[N],g[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void mul(int *a,int *b,int len){
	for(int i=0;i<len;i++)
		p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?len>>1:0);
	NTT(p,len,1),NTT(q,len,1);
	for(int i=0;i<len;i++) a[i]=1ll*p[i]*q[i]%mod;
	NTT(a,len,-1);
}
void inv(int *a,int *b,int len){
	for(int i=0;i<len;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=len;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
void diff(int *a,int *b,int len){
	for(int i=1;i<len;i++) b[i-1]=1ll*a[i]*i%mod; b[len-1]=0;
} 
void integ(int *a,int *b,int len){
	for(int i=0;i<len;i++) b[i+1]=1ll*a[i]*qpow(i+1,mod-2)%mod; b[0]=0;
}
void ln(int *a,int *b,int len){
	for(int i=0;i<(len<<1);i++) f[i]=g[i]=0;
	diff(a,f,len),inv(a,g,len),mul(f,g,len<<1),integ(f,b,len);
}
void exp(int *a,int *b,int len){
	for(int i=0;i<(len<<1);i++) b[i]=c[i]=0;
	b[0]=1;
	for(int i=2;i<=len;i<<=1){
		ln(b,c,i);
		for(int j=0;j<i;j++) c[j]=(-c[j]+a[j]+mod)%mod;
		c[0]=(c[0]+1)%mod,mul(b,c,i<<1);
	}
}
signed main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;i++) scanf("%d",&x),cnt[x]++;
	for(int i=1;i<=m;i++) if(cnt[i])
		for(int j=1;j<=m/i;j++) a[i*j]=(a[i*j]+mod-1ll*cnt[i]*qpow(j,mod-2)%mod)%mod;
	for(len=1;len<=m;len<<=1);
	exp(a,b,len),inv(b,a,len);
	for(int i=1;i<=m;i++) printf("%d\n",a[i]);
	return 0;
}

P3784 [SDOI2017] 遗忘的集合

有一个集合 \(S\)，每个数有无限个，给出组成 \(x\,(1\leq x\leq n)\) 的方案数 \(f_x\bmod p\)，求 \(S\)。

\(1\leq n<2^{18}\)，\(10^6\leq p<2^{30}\)，\(0\leq f_i<p\)。

付公主的背包的逆操作。

\(f\) 的 OGF：\(F(x)=\prod_{v\geq 1}a_v\cdot \large\frac{1}{1-x^v}\)。

两边同时取 \(\ln\) 可得 \(\ln F(x)=-\sum_{v\geq 1}a_v\cdot \ln(1-x^v)\)。而 \(-\ln(1-x^v)=\sum_{i\geq 0}\large\frac{x^{iv}}{i}\)，故 \(\ln F(x)=\sum_{v\geq 1}a_v\sum_{i\geq 0}\large\frac{x^{iv}}{i}\)。

设 \(G=\ln F\)，\([x^t]G(x)=\sum_{v\mid t}a_v\large\frac v t\)，\(t\cdot [x^t]G(x)=\sum_{v\mid t}a_v\cdot v\)，\(g_i\gets g_i\cdot i,a_i\gets a_i\cdot i\)，那么求 \(a\) 直接莫反即可，\(G=a*1\)，\(a=G*\mu\)。要写 MTT、多项式求 \(\ln\)。

#include<bits/stdc++.h>
#define ll long long
#define db long double
using namespace std;
const int N=(1<<19)+5;
int n,mod,p,a[N],b[N],c[N],len,r[N],A[N],B[N],f[N],g[N],cnt,pr[N],mu[N],ans;
bool vis[N]; 
db pi=acos(-1);
complex<db>a0[N],a1[N],b0[N],b1[N],s1[N],s2[N],s3[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void FFT(complex<db>*a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1; 
		complex<db>v(cos(2*pi/k),sin(2*pi/k)*op),w(1,0),x,y;
		for(int i=0;i<n;i+=k,w=1) 
			for(int j=i;j<i+m;j++)
				x=a[j],y=w*a[j+m],a[j]=x+y,a[j+m]=x-y,w*=v;
	}
	if(op==-1) for(int i=0;i<n;i++) a[i]/=n;
} 
void mul(int *a,int *b,int len){
	for(int i=0;i<len;i++)
		a0[i]=a[i]/p,a1[i]=a[i]%p,b0[i]=b[i]/p,b1[i]=b[i]%p,r[i]=r[i>>1]>>1|(i&1?len>>1:0);
	FFT(a0,len,1),FFT(a1,len,1),FFT(b0,len,1),FFT(b1,len,1);
	for(int i=0;i<len;i++)
		s1[i]=a0[i]*b0[i],s2[i]=a0[i]*b1[i]+a1[i]*b0[i],s3[i]=a1[i]*b1[i];
	FFT(s1,len,-1),FFT(s2,len,-1),FFT(s3,len,-1);
	for(int i=0;i<len;i++)
		a[i]=((ll)(s1[i].real()+0.5)%mod*p%mod*p%mod+(ll)(s2[i].real()+0.5)%mod*p%mod+(ll)(s3[i].real()+0.5)%mod)%mod;
}
void inv(int *a,int *b,int len){
	for(int i=0;i<len;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=len;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
void diff(int *a,int *b,int len){
	for(int i=1;i<len;i++) b[i-1]=1ll*a[i]*i%mod; b[len-1]=0;
} 
void integ(int *a,int *b,int len){
	for(int i=0;i<len;i++) b[i+1]=1ll*a[i]*qpow(i+1,mod-2)%mod; b[0]=0;
}
void ln(int *a,int *b,int len){
	for(int i=0;i<(len<<1);i++) f[i]=g[i]=0;
	diff(a,f,len),inv(a,g,len),mul(f,g,len<<1),integ(f,b,len);
}
signed main(){
	scanf("%d%d",&n,&mod),p=sqrt(mod),a[0]=1;
	for(int i=1;i<=n;i++) scanf("%d",&a[i]);
	vis[0]=vis[1]=1,mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]) pr[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*pr[j]<=n;j++){
			vis[i*pr[j]]=1;
			if(i%pr[j]==0){mu[i*pr[j]]=0;break;} 
			mu[i*pr[j]]=-mu[i];
		}
	}
	for(len=1;len<=n;len<<=1);
	ln(a,b,len);
	for(int i=1;i<=n;i++) b[i]=1ll*b[i]*i%mod;
	for(int i=1;i<=n;i++)
		for(int j=1;i*j<=n;j++) c[i*j]=(c[i*j]+1ll*(mu[i]+mod)%mod*b[j]%mod)%mod;
	for(int i=1;i<=n;i++) ans+=c[i]>0;
	printf("%d\n",ans);
	for(int i=1;i<=n;i++)
		if(c[i]) printf("%d ",i);
	return 0;
}

P5860 「SWTR-03」Counting Trees

\(n\) 个点，每个点度数给定，求有多少个子集使得这些点能组成树，对 \(998244353\) 取模。

\(2\leq n\leq 5\times 10^5\)，\(1\leq d_i\leq n\)。

根据 Prufer 序列，度数分别为 \(d_1,d_2,\cdots,d_m\) 的点可以构成一棵树，当且仅当 \(\sum_{i=1}^m d_i=2m-2\)，即 \(\sum_{i=1}^m (d_i-2)=-2\)。

问题转化为一个 01 背包问题，构造生成函数 \(F(x)=\prod_{i=1}^n (1+x^{d_i-2})\)，\(ans=[x^{-2}]F(x)\)。

和上一题差不多，只是 \((1+x^{v_1}+x^{2v_1}+\cdots)\cdots(1+x^{v_n}+x^{2v_n}+\cdots)\) 变成了 \((1+x^{v_1})(1+x^{v_2})\cdots(1+x^{v_n})\)。\(\ln\) 之后加在一起再 \(\exp\) 回来。

非正数下标？\(d_i=1\) 时，\(\prod_{d_i=1}(1+x^{-1})=\large\frac{\prod_{d_i=1}(x+1)}{x}\)，恰好 \(1=|d_i-2|\)。\(d_i=2\) 时，选 \(1\) 和选 \(x^0\) 都一样。设 \(k\) 表示有多少 \(d_i=1\)，那么 \(x^kF(x)=\prod_{i=1}^n (1+x^{|d_i-2|})\)，\(ans=[x^{k-2}]x^kF(x)=[x^{k-2}]\prod_{i=1}^n(1+x^{|d_i-2|})\)。

在实现时，\(d_i=2\) 的另算，最后答案乘上 \(2^{cnt_2}\)。

#include<bits/stdc++.h>
using namespace std;
const int N=(1<<20)+5,mod=998244353;
int n,k,x,c2,cnt[N],a[N],b[N],c[N],len,r[N],A[N],B[N],p[N],q[N],f[N],g[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void mul(int *a,int *b,int len){
	for(int i=0;i<len;i++)
		p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?len>>1:0);
	NTT(p,len,1),NTT(q,len,1);
	for(int i=0;i<len;i++) a[i]=1ll*p[i]*q[i]%mod;
	NTT(a,len,-1);
}
void inv(int *a,int *b,int len){
	for(int i=0;i<len;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=len;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
void diff(int *a,int *b,int len){
	for(int i=1;i<len;i++) b[i-1]=1ll*a[i]*i%mod; b[len-1]=0;
} 
void integ(int *a,int *b,int len){
	for(int i=0;i<len;i++) b[i+1]=1ll*a[i]*qpow(i+1,mod-2)%mod; b[0]=0;
}
void ln(int *a,int *b,int len){
	for(int i=0;i<(len<<1);i++) f[i]=g[i]=0;
	diff(a,f,len),inv(a,g,len),mul(f,g,len<<1),integ(f,b,len);
}
void exp(int *a,int *b,int len){
	for(int i=0;i<(len<<1);i++) b[i]=c[i]=0;
	b[0]=1;
	for(int i=2;i<=len;i<<=1){
		ln(b,c,i);
		for(int j=0;j<i;j++) c[j]=(-c[j]+a[j]+mod)%mod;
		c[0]=(c[0]+1)%mod,mul(b,c,i<<1);
	}
}
signed main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
		scanf("%d",&x),k+=(x==1),c2+=(x==2),cnt[abs(x-2)]++;
	for(int i=1;i<=k;i++) if(cnt[i])
		for(int j=1;j<=k/i;j++)
			a[i*j]=(a[i*j]+1ll*(j&1?1:mod-1)*cnt[i]%mod*qpow(j,mod-2)%mod)%mod;
	for(len=1;len<=k;len<<=1);
	exp(a,b,len),printf("%lld\n",1ll*qpow(2,c2)*b[k-2]%mod);
	return 0;
}

P4841 [集训队作业2013]城市规划

求 \(n\) 个点的简单有标号无向连通图个数，对 \(1004535809\) 取模。

\(1\leq n\leq 1.3\times 10^5\)。

设 \(g_i\) 表示 \(i\) 个点的无向图的个数（显然 \(g_i=2^{\large\binom n 2}\)），\(f_i\) 表示 \(i\) 个点的无向连通图个数。

枚举 \(1\) 所在的连通块的大小，\(g_n=\sum_{i=1}^n\large\binom{n-1}{i-1}\normalsize f_ig_{n-i}\)。

拆组合数，\(g_n=\sum_{i=1}^n\large\frac{(n-1)!}{(i-1)!(n-i)!}\normalsize f_ig_{n-i}\)，变形 \(\large\frac{g_n}{(n-1)!}\normalsize=\sum_{i=1}^n \large\frac{f_i}{(i-1)!}\normalsize \times\large\frac{g_{n-i}}{(n-i)!}\)。

设 \(a_i=\large\frac{g_i}{(i-1)!}\)，\(b_i=\large\frac{f_i}{(i-1)!}\)，\(c_i=\large\frac{g_i}{i!}\)，则 \(a=b\times c\)，\(b=a\times c^{-1}\)，\(ans=b_n\times (n-1)!\)。多项式求逆即可。

#include<bits/stdc++.h>
using namespace std;
const int N=(1<<18)+5,mod=1004535809;
int n,a[N],b[N],c[N],fac[N],iv[N],len,r[N],A[N],B[N],p[N],q[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void mul(int *a,int *b,int len){
	for(int i=0;i<len;i++)
		p[i]=a[i],q[i]=b[i],r[i]=r[i>>1]>>1|(i&1?len>>1:0);
	NTT(p,len,1),NTT(q,len,1);
	for(int i=0;i<len;i++) a[i]=1ll*p[i]*q[i]%mod;
	NTT(a,len,-1);
}
void inv(int *a,int *b,int len){
	for(int i=0;i<len;i++) b[i]=0;
	b[0]=qpow(a[0],mod-2);
	for(int i=2;i<=len;i<<=1){
		for(int j=0;j<(i<<1);j++) A[j]=j<i?a[j]:0,B[j]=j<(i>>1)?b[j]:0; 
		mul(B,B,i),mul(A,B,i<<1);
		for(int j=0;j<i;j++) b[j]=(2ll*b[j]%mod-A[j]+mod)%mod;
	}
}
signed main(){
	scanf("%d",&n),fac[0]=iv[0]=iv[1]=c[0]=1;
	for(int i=2;i<=n;i++) iv[i]=1ll*iv[mod%i]*(mod-mod/i)%mod;
	for(int i=1;i<=n;i++){ 
		fac[i]=1ll*fac[i-1]*i%mod,iv[i]=1ll*iv[i-1]*iv[i]%mod; 
		a[i]=1ll*qpow(2,1ll*i*(i-1)/2%(mod-1))*iv[i-1]%mod,c[i]=1ll*qpow(2,1ll*i*(i-1)/2%(mod-1))*iv[i]%mod;
	} 
	for(len=1;len<=n;len<<=1);
	inv(c,b,len),mul(b,a,len),printf("%lld\n",1ll*b[n]*fac[n-1]%mod);
	return 0;
}

CF623E Transforming Sequence（*3300）

给出 \(n,k\)，求有多少个长度为 \(n\) 的序列 \(a\) 满足，\(0<a_i<2^k\)，且 \(a_i\) 的前缀或数组 \(b\) 严格单调递增。对 \(10^9+7\) 取模。

\(1\leq n\leq 10^{18}\)，\(1\leq k\leq 3\times 10^4\)。

显然 \(n>k\) 时答案为 \(0\)。

设 \(f_{i,j}\) 表示 \(b_i\) 有 \(j\) 位为 \(1\) 的方案数。则 \(f_{i,j}=\sum_{k=0}^{j-1}f_{i-1,k}\times \large\binom j k\normalsize\times 2^k\)，\(ans=\sum_{i=1}^k f_{n,i}\times\large\binom k i\)。

\(f_i\to f_{i+1}\)，要 \(n\) 次才能推到 \(f_n\)。假如改为 \(f_i\to f_{i+j}\) 呢？

\(f_{x+y,j}=\sum_{k=0}^{j-1}f_{x,k}\times f_{y,j-k}\times \large\binom j k\normalsize\times2^{y\times k}\)，这样就可以倍增 \(f_i\to f_{2i}\) 了。

转移可以用卷积优化。\(f_{x+y,j}=j!\sum_{k=0}^{j-1}\large\frac{f_{x,k}2^{y\times k}}{k!}\normalsize\times \large\frac{f_{y,j-k}}{(j-k)!}\)。

要写任意模数 NTT。时间复杂度 \(\mathcal O(k\log n\log k)\)。

#include<bits/stdc++.h>
#define ll long long
#define db long double
using namespace std;
const int N=(1<<18)+5,mod=1e9+7;
int k,p=sqrt(mod),len,r[N],fac[N],inv[N],a[N],b[N],f[N],g[N],ans;
ll n;
db pi=acos(-1);
complex<db>a0[N],a1[N],b0[N],b1[N],s1[N],s2[N],s3[N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void FFT(complex<db>*a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1; 
		complex<db>v(cos(2*pi/k),sin(2*pi/k)*op),w(1,0),x,y;
		for(int i=0;i<n;i+=k,w=1) 
			for(int j=i;j<i+m;j++)
				x=a[j],y=w*a[j+m],a[j]=x+y,a[j+m]=x-y,w*=v;
	}
	if(op==-1) for(int i=0;i<n;i++) a[i]/=n;
} 
void mul(int *a,int *b,int len){
	for(int i=0;i<len;i++)
		a0[i]=a[i]/p,a1[i]=a[i]%p,b0[i]=b[i]/p,b1[i]=b[i]%p,r[i]=r[i>>1]>>1|(i&1?len>>1:0);
	FFT(a0,len,1),FFT(a1,len,1),FFT(b0,len,1),FFT(b1,len,1);
	for(int i=0;i<len;i++)
		s1[i]=a0[i]*b0[i],s2[i]=a0[i]*b1[i]+a1[i]*b0[i],s3[i]=a1[i]*b1[i];
	FFT(s1,len,-1),FFT(s2,len,-1),FFT(s3,len,-1);
	for(int i=0;i<len;i++)
		a[i]=((ll)(s1[i].real()+0.5)%mod*p%mod*p%mod+(ll)(s2[i].real()+0.5)%mod*p%mod+(ll)(s3[i].real()+0.5)%mod)%mod;
}
void merge(int *a,int *b,ll x){
	for(int i=0;i<len;i++) f[i]=g[i]=0;
	for(int i=0;i<=k;i++)
		f[i]=1ll*a[i]*qpow(2,x*i)%mod*inv[i]%mod,g[i]=1ll*b[i]*inv[i]%mod;
	mul(f,g,len);
	for(int i=0;i<=k;i++) a[i]=1ll*f[i]*fac[i]%mod;
}
int C(int n,int m){return n<m?0:1ll*fac[n]*inv[m]%mod*inv[n-m]%mod;}
signed main(){
	scanf("%lld%d",&n,&k);
	if(n>k) puts("0"),exit(0);
	fac[0]=inv[0]=inv[1]=1;
	for(int i=2;i<=k;i++) inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
	for(int i=1;i<=k;i++)
		fac[i]=1ll*fac[i-1]*i%mod,inv[i]=1ll*inv[i-1]*inv[i]%mod;
	for(len=1;len<=2*k;len<<=1);
	a[0]=1,fill(b+1,b+1+k,1);	//初始时 a,b 分别表示 f[0],f[1]
	for(int x=1;n;n>>=1,merge(b,b,x),x<<=1)
		if(n&1) merge(a,b,x);
	for(int i=1;i<=k;i++) ans=(ans+1ll*a[i]*C(k,i)%mod)%mod;
	printf("%d\n",ans);
	return 0;
}

CF755G PolandBall and Many Other Balls（*3200）

对所有 \(i\in[1,k]\)，将排成一行的 \(n\) 个球分成 \(i\) 个组，一组可以是一个球或两个相邻的球，允许有些球不在任何一组里面，求方案数。对 \(998244353\) 取模。

\(1\leq n\leq 10^9\)，\(1\leq k\leq 2^{15}\)。

设 \(f_{i,j}\) 表示前 \(i\) 个球分成 \(j\) 组的方案数，\(f_{i,j}=f_{i-1,j}+f_{i-1,j-1}+f_{i-2,j-1}\)。

借鉴上一题的思路，分接缝处组成一组和接缝处不组成一组讨论，\(f_{x+y,i}=\sum_{k=0}^i f_{x,k}f_{y,i-k}+\sum_{k=0}^{i-1}f_{x-1,k}\times f_{y-1,i-k-1}\)。

记 \(F_i\) 表示 \(f_i\) 的生成函数，则 \(F_{x+y}=F_xF_y+xF_{x-1}F_{y-1}\)。

考虑在求 \(F_x\) 的同时把 \(F_{x-1}\) 也求出来。

求 \(F_{x+y}\) 的同时要求 \(F_{x+y-1}\)。\(F_{x+y-1}=F_xF_{y-1}+xF_{x-1}F_{y-2}\)，所以再记一个 \(F_{x-2}\) 即可。而 \(F_{x+y-2}=F_{x-1}F_{y-1}+xF_{x-2}F_{y-2}\)，于是可以 \(\mathcal O(n\log n)\) 地用 \((F_x,F_{x-1},F_{x-2})\) 和 \((F_y,F_{y-1},F_{y-2})\) 求出 \((F_{x+y},F_{x+y-1},F_{x+y-2})\) 了。再套个倍增即可。

时间复杂度 \(\mathcal O(n\log^2 n)\)。可能需要卡常。

#include<bits/stdc++.h>
using namespace std;
const int N=(1<<16)+5,mod=998244353;
int n,k,a[3][N],b[3][N],A[3][N],B[3][N],len,r[N],f[5][N];
int qpow(int x,int n){
	int ans=1;
	for(;n;n>>=1,x=1ll*x*x%mod) if(n&1) ans=1ll*ans*x%mod;
	return ans;
}
void NTT(int *a,int n,int op){
	for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]);
	for(int k=2;k<=n;k<<=1){
		int m=k>>1,v=qpow(~op?3:qpow(3,mod-2),(mod-1)/k),w=1,x,y; 
		for(int i=0;i<n;i+=k,w=1)
			for(int j=i;j<i+m;j++)
				x=a[j],y=1ll*w*a[j+m]%mod,a[j]=(x+y)%mod,a[j+m]=(x-y+mod)%mod,w=1ll*w*v%mod;
	}
	if(op==-1){
		int inv=qpow(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%mod;
	}
} 
void merge(int a[3][N],int b[3][N]){
	for(int i=0;i<3;i++){ 
		for(int j=0;j<len;j++) A[i][j]=a[i][j],B[i][j]=b[i][j];
		NTT(A[i],len,1),NTT(B[i],len,1);
	}
	for(int i=0;i<len;i++){ 
		f[0][i]=1ll*A[0][i]*B[0][i]%mod,f[1][i]=1ll*A[1][i]*B[1][i]%mod;
		f[2][i]=1ll*A[0][i]*B[1][i]%mod,f[3][i]=1ll*A[1][i]*B[2][i]%mod; 
		f[4][i]=1ll*A[2][i]*B[2][i]%mod;
	} 
	for(int i=0;i<5;i++) NTT(f[i],len,-1);
	for(int i=0;i<=k;i++){ 
		a[0][i]=(f[0][i]+(i?f[1][i-1]:0))%mod;
		a[1][i]=(f[2][i]+(i?f[3][i-1]:0))%mod;
		a[2][i]=(f[1][i]+(i?f[4][i-1]:0))%mod;
	} 
}
signed main(){
	scanf("%d%d",&n,&k);
	for(len=1;len<=2*k;len<<=1);
	for(int i=0;i<len;i++) r[i]=r[i>>1]>>1|(i&1?len>>1:0);
	a[0][0]=1,b[0][0]=b[0][1]=b[1][0]=1;
	for(;n;n>>=1,merge(b,b)) if(n&1) merge(a,b);
	for(int i=1;i<=k;i++) printf("%d ",a[0][i]);
	return 0;
}

卡常之前：

//TLE
void merge(int a[3][N],int b[3][N]){
	mul(a[0],b[0],f),mul(a[1],b[1],g);
	for(int i=0;i<=k;i++) c[0][i]=(f[i]+(i?g[i-1]:0))%mod;
	mul(a[0],b[1],f),mul(a[1],b[2],g);
	for(int i=0;i<=k;i++) c[1][i]=(f[i]+(i?g[i-1]:0))%mod;
	mul(a[1],b[1],f),mul(a[2],b[2],g);
	for(int i=0;i<=k;i++) c[2][i]=(f[i]+(i?g[i-1]:0))%mod;
	for(int i=0;i<3;i++)
		for(int j=0;j<=k;j++) a[i][j]=c[i][j];
}
posted @ 2021-03-13 09:36  maoyiting  阅读(419)  评论(0)    收藏  举报