题解 洛谷 P5434 【有标号荒漠计数】

因为是有标号计数，考虑用指数型生成函数来解决。

设有标号荒漠的生成函数为 (F(x))，有标号仙人掌的生成函数为 (f(x))，因为荒漠是仙人掌所组成的集合，得：

[F(x) = exp (f(x)) ]

解决有标号仙人掌的计数问题后，就能解决有标号荒漠的计数问题。

先在仙人掌上确定一点来考虑，和该点相邻的边有在环上和不在环上两种情况，即与该点相邻的元素为独立边和环。独立边可看作单独的一个仙人掌，得独立边的生成函数为 (f(x))。环的大小至少为 (3)，环上每个点也都对应单独的一个仙人掌，正方向反方向都会计算一遍，所以要除 (2)，得环的生成函数为 (sumlimits_{i geqslant 2 } frac{f^i(x)}{2})。

该点相邻的所有元素可以看作一个集合，同时还要考虑确定的那个点，得：

[f(x) = x exp(f(x) + sum_{i geqslant 2 } frac{f^i(x)}{2}) = x exp(frac{f^2(x) - 2f(x)}{2f(x)-2}) ]

通过牛顿迭代来求解。设已经求得 (f_0(x))，满足 (f_0(x) equiv f(x) pmod{x^n})，设 (g(f(x))= x exp(frac{f^2(x) - 2f(x)}{2f(x)-2}) - f(x),h(x) = x exp(frac{f_0^2(x) - 2f_0(x)}{2f_0(x)-2}))，得：

[egin{aligned} f(x) & equiv f_0(x) - frac{g(f_0(x))}{{g}'(f_0(x))} pmod{x^{2n}} \ & equiv f_0(x) - frac{h(x)-f_0(x)}{h(x){(frac{f_0^2(x) - 2f_0(x)}{2f_0(x)-2})}'-1} pmod{x^{2n}} \ & equiv f_0(x) - frac{2h(x)-2f_0(x)}{h(x)(frac{1}{(f_0(x)-1)^2}+1)-2} pmod{x^{2n}}end{aligned} ]

因为之前的推导是确定了仙人掌上的一个点，所以对于 (n) 个点的仙人掌，一个方案会计算 (n) 次，要除以 (n)。

(code:)

#include<bits/stdc++.h>
#define maxn 800010
#define P 998244353
#define G 3
using namespace std;
typedef long long ll;
template<typename T> inline void read(T &x)
{
    x=0;char c=getchar();bool flag=false;
    while(!isdigit(c)){if(c=='-')flag=true;c=getchar();}
    while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
    if(flag)x=-x;
}
ll n;
ll rev[maxn],f[maxn],g[maxn],fac[maxn],iv[maxn],w[22][maxn];
ll qp(ll x,ll y)
{
    ll v=1;
    while(y)
    {
        if(y&1) v=v*x%P;
        x=x*x%P,y>>=1;
    }
    return v;
}
int calc(int n)
{
    int lim=1;
    while(lim<=n) lim<<=1;
    for(int i=0;i<lim;++i)
        rev[i]=(rev[i>>1]>>1)|((i&1)?lim>>1:0);
    return lim;
}
void NTT(ll *a,int lim,int type)
{
    for(int i=0;i<lim;++i)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int len=1,t=1;len<lim;len<<=1,++t)
    {
        for(int i=0;i<lim;i+=len<<1)
        {
            for(int j=0;j<len;++j)
            {
                ll x=a[i+j],y=w[t][j]*a[i+j+len]%P;
                a[i+j]=(x+y)%P,a[i+j+len]=(x-y+P)%P;
            }
        }
    }
    if(type==1) return;
    ll inv=qp(lim,P-2);
    for(int i=0;i<lim;++i) a[i]=a[i]*inv%P;
    reverse(a+1,a+lim);
}
void Inv(int deg,ll *a,ll *b)
{
    static ll t[maxn];
    if(deg==1)
    {
        b[0]=qp(a[0],P-2);
        return;
    }
    Inv((deg+1)>>1,a,b);
    int lim=calc(deg<<1);
    for(int i=0;i<deg;++i) t[i]=a[i];
    for(int i=deg;i<lim;++i) t[i]=b[i]=0;
    NTT(t,lim,1),NTT(b,lim,1);
    for(int i=0;i<lim;++i) b[i]=b[i]*(2-t[i]*b[i]%P+P)%P;
    NTT(b,lim,-1);
    for(int i=deg;i<lim;++i) b[i]=0;
}
void Ln(int deg,ll *a,ll *b)
{
    static ll inva[maxn],dera[maxn];
    Inv(deg,a,inva);
    for(int i=0;i<deg-1;++i) dera[i]=a[i+1]*(i+1)%P;
    dera[deg-1]=0;
    int lim=calc(deg<<1);
    for(int i=deg;i<lim;++i) dera[i]=inva[i]=0;
    NTT(dera,lim,1),NTT(inva,lim,1);
    for(int i=0;i<lim;++i) b[i]=dera[i]*inva[i]%P;
    NTT(b,lim,-1);
    for(int i=deg-1;i>=1;--i) b[i]=b[i-1]*qp(i,P-2)%P;
    b[0]=0;
    for(int i=deg;i<lim;++i) b[i]=0;
}
void Exp(int deg,ll *a,ll *b)
{
    static ll t[maxn],lnb[maxn];
    if(deg==1)
    {
        b[0]=1;
        return;
    }
    Exp((deg+1)>>1,a,b),Ln(deg,b,lnb);
    int lim=calc(deg<<1);
    for(int i=0;i<deg;++i) t[i]=(a[i]-lnb[i]+P)%P;
    for(int i=deg;i<lim;++i) t[i]=b[i]=0;
    NTT(t,lim,1),NTT(b,lim,1);
    for(int i=0;i<lim;++i) b[i]=b[i]*(1+t[i])%P;
    NTT(b,lim,-1);
    for(int i=deg;i<lim;++i) b[i]=0;
}
void work(int deg,ll *f)
{
    static ll h[maxn],t1[maxn],t2[maxn],t3[maxn];
    if(deg==1)
    {
        f[0]=0;
        return;
    }
    work((deg+1)>>1,f);
    int lim=calc(deg<<1);
    for(int i=0;i<deg;++i) t1[i]=f[i],t2[i]=f[i]*2;
    NTT(t1,lim,1);
    for(int i=0;i<lim;++i) t1[i]=(t1[i]*t1[i]-2*t1[i]+P)%P;
    t2[0]-=2,Inv(deg,t2,t3),NTT(t3,lim,1);
    for(int i=0;i<lim;++i) t3[i]=t1[i]*t3[i]%P,t1[i]=t1[i]+1;
    NTT(t3,lim,-1),Exp(deg,t3,h),NTT(t1,lim,-1),Inv(deg,t1,t2);
    for(int i=deg-1;i>=1;--i) h[i]=h[i-1];
    h[0]=0;
    for(int i=0;i<deg;++i) t1[i]=(2*h[i]-2*f[i]+P)%P;
    t2[0]++,NTT(t2,lim,1),NTT(h,lim,1);
    for(int i=0;i<lim;++i) t2[i]=t2[i]*h[i]%P;
    NTT(t2,lim,-1),t2[0]-=2,Inv(deg,t2,t3),NTT(t1,lim,1),NTT(t3,lim,1);
    for(int i=0;i<lim;++i) t1[i]=t1[i]*t3[i]%P;
    NTT(t1,lim,-1);
    for(int i=0;i<deg;++i) f[i]=(f[i]-t1[i]+P)%P;
    for(int i=deg;i<lim;++i) f[i]=0;
}
void init()
{
    fac[0]=fac[1]=iv[1]=1;
    for(int i=2;i<=n;++i)
        fac[i]=fac[i-1]*i%P,iv[i]=(P-P/i)*iv[P%i]%P;
    for(int i=1;i<=19;++i) w[i][0]=1,w[i][1]=qp(G,(P-1)/(1<<i));
    for(int i=1;i<=19;++i)
        for(int j=2;j<1<<i;++j)
            w[i][j]=w[i][j-1]*w[i][1]%P;
}
int main()
{
    read(n),init(),work(n+1,f);
    for(int i=0;i<=n;++i) f[i]=f[i]*iv[i]%P;
    Exp(n+1,f,g),printf("%lld",g[n]*fac[n]%P);
	return 0;
}