HDU 5896 Running King（dp+NTT+CDQ分治）

Description

求$n$个点的无向有环图个数

Input

第一行一整数$T$表示用例组数，每组用例输入一整数$n(1\le T\le 10,1\le n\le 2\cdot 10^5)$

Output

输出$n$个点的无向有环图个数，结果模$1004535809$

Sample Input

3
2
3
4

Sample Output

0
1
26

Solution

$n$个点的无向有环图个数$=2^{\frac{n(n-1)}{2}}-n$个点的无向森林个数

令$f(n)$为$n$个点的无向森林个数，$g(n)$为$n$个点的无向树的个数，则$g(n)=n^{n-2}$(每棵$n$个节点无向树唯一对应一个长度为$n-2$的$Purfer$序列)

枚举第一个点所在树的节点数有转移$f(n)=g(n)+\sum\limits_{i=1}^{n-1}C_{n-1}^{i-1}g(i)f(n-i)$，令$f(0)=1$，则有$\frac{f(n)}{(n-1)!}=\sum\limits_{i=1}^n\frac{i^{i-1}}{i!}\frac{f(n-i)}{(n-i)!}$，$CDQ$分治即可

假设当前求解$f(l),...,f(r)$，而$f(l),...,f(mid)$已经求出，考虑$f(l),...,f(mid)$对$f(mid+1),...,f(r)$的贡献，令$x(i)=\frac{(i+1)^{i}}{(i+1)!},y(i)=\frac{f(i+l)}{(i+l)!},i=0,...,mid-l$，对$mid+1\le t\le r$，$f(l),...,f(mid)$对$\frac{f(t)}{(t-1)!}$的贡献为$h(t)=\sum\limits_{i=l}^{mid}x(t-i-1)\frac{f(i)}{i!}=\sum\limits_{i=0}^{t-l-1}x(t-l-1-i)y(i)=(x*y)(t-l-1)$，其中$x*y$表示$x,y$两序列的卷积，$NTT$即可快速求卷积

Code

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
using namespace std;
typedef long long ll;
#define maxn 222222
typedef long long ll;
const ll mod=1004535809ll;
const ll g=3;
ll mod_pow(ll a,ll b,ll p)
{
    a%=p;
    ll ans=1;
    while(b)
    {
        if(b&1)ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans;
}
void change(ll *x,int len) 
{
    for(int i=1,j=len/2;i<len-1;i++)
    {
        if(i<j)swap(x[i],x[j]);
        int k=len/2;
        while(j>=k)
        {
            j-=k;
            k/=2;
        }
        if(j<k)j+=k;
    }
}
void ntt(ll *x,int len,int sta) 
{
    change(x,len);
    for(int m=2;m<=len;m<<=1) 
    {
        ll Wn=mod_pow(g,(mod-1)/m,mod);
        if(sta==-1)Wn=mod_pow(Wn,mod-2,mod);
        for(int i=0;i<len;i+=m) 
        {
            ll W=1;
            for(int j=i;j<i+m/2;j++) 
            {
                ll x1=x[j],x2=W*x[j+m/2]%mod;
                x[j]=(x1+x2)%mod,x[j+m/2]=(x1-x2+mod)%mod;
                W=W*Wn%mod;
            }
        }
    }
    if(sta==-1)
    {
        int temp=mod_pow(len,mod-2,mod);
        for(int i=0;i<len;i++)
            x[i]=x[i]*temp%mod;
    }
}
ll mod_pow(ll a,ll b)
{
    ll ans=1ll;
    while(b)
    {
        if(b&1)ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
} 
ll n,a[maxn],dp[maxn];
ll x[2*maxn],y[2*maxn];
ll f[maxn],h[maxn],inv[maxn];
void init()
{
    f[0]=1,inv[0]=1;
    for(int i=1;i<maxn;i++)f[i]=f[i-1]*i%mod,inv[i]=mod_pow(f[i],mod-2);
    for(int i=0;i<maxn;i++)h[i]=mod_pow(i+1,i);
}
void deal(int l,int r)
{
    if(l==r)return ;
    int mid=(l+r)>>1;
    deal(l,mid);
    int len=1;
    while(len<=r-l+1)len<<=1;
    for(int i=0;i<len;i++)
    {
        x[i]=a[i];
        if(i+l<=mid)y[i]=dp[i+l]*inv[i+l]%mod;
        else y[i]=0;
    }
    ntt(x,len,1),ntt(y,len,1);
    for(int i=0;i<len;i++)x[i]=x[i]*y[i]%mod;
    ntt(x,len,-1);
    for(int i=mid+1;i<=r;i++)
        dp[i]+=x[i-l-1]*f[i-1]%mod,dp[i]%=mod;
    deal(mid+1,r);
}
int main()
{
    int T;
    scanf("%d",&T);
    init();
    while(T--)
    {
        scanf("%d",&n);
        for(int i=0;i<=n;i++)a[i]=h[i]*inv[i+1]%mod;
        memset(dp,0,sizeof(dp));
        dp[0]=1;
        deal(0,n);
        ll ans=mod_pow(2,1ll*n*(n-1)/2)-dp[n];
        printf("%I64d\n",(ans%mod+mod)%mod);
    }
    return 0;
}