hdu5896Running King

链接：http://acm.hdu.edu.cn/showproblem.php?pid=5896

题意：给定多个n，求有多少n个点的无向无重边无自环的图满足有环，$mod 1004535809$。

分析：其实看到这个$mod$是一个可ntt的费马素数就应该有点方向。直接构图我们不好求，我们反向求有多少n个点的森林设为$f(n)$，那么我们就可以用n个点的图的总数$2^{n*(n-1)/2}-f(n)$。那么怎么快速求出$f(n)$呢？我们选择增量如果之前已经求出$f(1) \ldots f(n-1)$当我们添加标号为n的这个点时我们枚举有多少个点和这个点n构成一棵树，比如组成一棵i个节点的书时我们要在前面的$n-1$个点中选出$i-1$个点，然后根据prufer定理i个点的树有$i^{i-2}$种，这样我们就得到了公式$f(n)=\sum_{i=1}^{n}f(n-i)*C(n-1,i-1)*i^{i-2}$，分解一下我们可以得到$\frac {f(n)}{n!}=\frac{1}{n}\sum_{i=1}^{n}\frac {f(n-i)}{(n-i)!} *\frac{i^{i-2}}{(i-1)!}$，其实看到这里我们就能看出这是一个分治ntt就行了，令$g(n)=\frac{f(n)}{n!}$，$h(i)=\frac{i^{i-2}}{(i-1)!}$就会有$g(n)=\frac{1}{n}\sum_{i=1}^{n}g(n-i)*h(i)$。如果是不懂分治fft的可以转hdu5730。

#include<map>
#include<set>
#include<cmath>
#include<queue>
#include<bitset>
#include<math.h>
#include<vector>
#include<string>
#include<stdio.h>
#include<cstring>
#include<iostream>
#include<algorithm>
#pragma comment(linker, "/STACK:102400000,102400000")
using namespace std;
const int N=200010;
const int M=50010;
const int mod=1000000007;
const int MOD1=1000000007;
const int MOD2=1000000009;
const double EPS=0.00000001;
typedef long long ll;
const ll MOD=1004535809;
const int INF=1000000010;
const ll MAX=1ll<<55;
const double eps=1e-5;
const double inf=~0u>>1;
const double pi=acos(-1.0);
typedef long double db;
typedef unsigned int uint;
typedef unsigned long long ull;
ll q[15],g[N],h[N],f[N],p[N],inv[N];
ll G=3,P=1004535809,a[3*N],b[3*N],wn[25];
ll q_pow(ll a,ll b,ll M) {
    ll ret=1;a%=M;
    while (b) {
        if (b&1) ret=ret*a%M;
        a=a*a%M;b>>=1;
    }
    return ret;
}
void getwn() {
    for (int i=0;i<25;i++) wn[i]=q_pow(G,(P-1)/(1<<i),P);
}
void NTT(ll x[],int n,int rev) {
    int i,j,k,t,ds,id=0;
    ll w,u,v;
    for (i=1,j=n>>1,k=n>>1;i<n-1;i++,k=n>>1) {
        if (i<j) swap(x[i],x[j]);
        while (j>=k) { j-=k;k>>=1; }
        if (j<k) j+=k;
    }
    for (i=2,ds=1;i<=n;i<<=1,ds++)
        for (j=0;j<n;j+=i) {
            w=1;
            for (k=j;k<j+i/2;k++) {
                u=x[k]%P;v=w*x[k+i/2]%P;
                x[k]=(u+v)%P;
                x[k+i/2]=(u-v+P)%P;
                w=w*wn[ds]%P;
            }
        }
    if (rev==-1) {
        for (i=1;i<n/2;i++) swap(x[i],x[n-i]);
        w=q_pow(n,P-2,P);
        for (i=0;i<n;i++) x[i]=x[i]*w%P;
    }
}
void cdq(int l,int r) {
    if (l==r) { (g[l]+=h[l]*inv[l]%MOD)%=MOD;return ; }
    int i,len=1,mid=(l+r)>>1,len1=mid-l+1,len2=r-l;
    cdq(l,mid);
    while (len<len1+len2) len<<=1;
    for (i=0;i<len1;i++) a[i]=g[l+i];
    for (i=len1;i<len;i++) a[i]=0;
    for (i=0;i<len2;i++) b[i]=h[i+1];
    for (i=len2;i<len;i++) b[i]=0;
    NTT(a,len,1);NTT(b,len,1);
    for (i=0;i<len;i++) (a[i]*=b[i])%=P;
    NTT(a,len,-1);
    for (i=mid+1;i<=r;i++) (g[i]+=inv[i]*a[i-l-1]%MOD)%=MOD;
    cdq(mid+1,r);
}
int main()
{
    int i,T;ll n=0;
    scanf("%d", &T);
    for (i=1;i<=T;i++) scanf("%I64d", &q[i]),n=max(n,q[i]);
    inv[1]=1;f[1]=1;p[1]=1;g[0]=1;h[1]=1;getwn();
    for (i=2;i<=n;i++) {
        inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;f[i]=(f[i-1]*inv[i])%MOD;
        p[i]=p[i-1]*i%MOD;h[i]=q_pow(i,i-2,MOD)*f[i-1]%MOD;
    }
    cdq(1,n);
    for (i=1;i<=T;i++) printf("%I64d\n", (q_pow(2ll,q[i]*(q[i]-1)/2,MOD)-g[q[i]]*p[q[i]]%MOD+MOD)%MOD);
    return 0;
}