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

令f_n表示n个点的方案数，则有递推公式: f_n=2^C(n,2)-∑_i(i>1&&i<n) C(n-1,i-1)*f_i*2^C(n-i,2) (用总方案数减去不连通的方案数)

两边都乘上(n-1)!^-1然后再移项整理得:2^C(n,2)*(n-1)^-1=∑_i(i>1&&i<=n)f_i*(i-1)!*2^C(n-i,2)*(n-i)!

令A_i=f_i*(i-1)!,B_i=2^C(i,2)*i!,C_i=2^C(n,2)*(n-1)!^-1,则有A×B=C,A=C*B^-1

这样就只用算出B在mod xⁿ⁺¹意义下的逆元，就可以计算A了。

Portal

#include <stdio.h>
#include <string.h>
#include <algorithm>
#define G 3
#define DFT(a,n) fft(a,n,0)
#define IDFT(a,n) fft(a,n,1)
#define INV(x,P) qpow(x,P-2,P)
#define rep(i,a,b) for(int i=a,s=b;i<s;++i)
using namespace std;

const int MAXN=(1<<18)+5;
typedef long long ll;
const int P=1004535809;

inline ll qpow(ll x,int k,int M){
	ll ans=1;
	for(;k;x=x*x%M,k>>=1)
		if(k&1)ans=ans*x%M;
	return ans;
}

int r[MAXN],w[2][MAXN],rn;

inline void init_fft(int n,int b){
	rn=INV(n,P);
	int g=qpow(G,(P-1)/n,P);
	w[0][0]=w[0][n]=1;
	rep(i,1,n)w[0][i]=1ll*w[0][i-1]*g%P;
	rep(i,0,n+1)w[1][i]=w[0][n-i];
	rep(i,1,n+1){
		r[i]=(r[i>>1]>>1)|((i&1)<<(b-1));
	}
}

inline void fft(int*a,int n,int t){
	rep(i,1,n)if(i<r[i])swap(a[i],a[r[i]]);
	for(int i=2;i<=n;i<<=1)
		for(int j=0;j<n;j+=i)
			for(int k=0;k<i>>1;++k){
				int v=1ll*a[(i>>1)+j+k]*w[t][n/i*k]%P;
				a[(i>>1)+j+k]=(a[j+k]+P-v)%P;
				a[j+k]=(a[j+k]+v)%P;
			}
	if(t)rep(i,0,n)a[i]=1ll*a[i]*rn%P;
}

inline void inv(int*a,int*b,int n,int bit){
	b[0]=INV(a[0],P),b[1]=0;
	int l=1;
	static int c[MAXN];
	rep(i,1,bit+1){
		l<<=1;
		rep(j,0,l)c[j]=a[j];
		rep(j,l,l<<1)b[j]=c[j]=0;
		init_fft(l<<1,i+1);
		DFT(b,l<<1),DFT(c,l<<1);
		rep(j,0,l<<1)
			b[j]=1ll*b[j]*(2+P-1ll*b[j]*c[j]%P)%P;
		IDFT(b,l<<1);
		rep(j,l,l<<1)b[j]=0;
	}
}

int n,l,b;

int fac[MAXN],inv_fac[MAXN];

inline void Init(){
	scanf("%d",&n);
	l=n+1,b=1;
	while(1<<b<l)++b;
	l=1<<b;
	fac[0]=1;
	rep(i,1,l)fac[i]=1ll*fac[i-1]*i%P;
	inv_fac[0]=inv_fac[1]=1;
	rep(i,2,l){
		int a=P/i,b=P%i;
		inv_fac[i]=1ll*a*(P-inv_fac[b])%P;
	}
	rep(i,2,l)
		inv_fac[i]=1ll*inv_fac[i]*inv_fac[i-1]%P;
}

inline ll C(ll n){return n*(n-1)>>1;}

int a[MAXN],c[MAXN];

inline void Solve(){
	rep(i,0,n+1)
		c[i]=qpow(2,C(i)%(P-1),P)*inv_fac[i]%P;
	inv(c,a,l,b);
	c[0]=0;
	rep(i,1,n+1)
		c[i]=qpow(2,C(i)%(P-1),P)*inv_fac[i-1]%P;
	rep(i,l,l<<1)a[i]=c[i]=0;
	init_fft(l<<1,b+1);
	DFT(a,l<<1),DFT(c,l<<1);
	rep(i,0,l<<1)a[i]=1ll*a[i]*c[i]%P;
	IDFT(a,l<<1);
	int ans=1ll*a[n]*fac[n-1]%P;
	printf("%d\n",ans);
}

int main(){
	Init(); Solve(); return 0;
}