[n点无向图个数 分治FFT || 多项式求逆] BZOJ 3456 城市规划

然后就分治FFT了

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

inline char nc(){
  static char buf[100000],*p1=buf,*p2=buf;
  if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
  return *p1++;
}

inline void read(int &x){
  char c=nc(),b=1;
  for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
  for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

inline void read(ll &x){
  char c=nc(),b=1;
  for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
  for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

const int MAXN=300005;
const ll P=1004535809;
const ll G=3;

ll inv[MAXN],fac[MAXN],Inv[MAXN],C[MAXN];

inline ll Pow(ll a,ll b){
  ll ret=1;
  for (;b;b>>=1,a=a*a%P)
    if (b&1)
      ret=ret*a%P;
  return ret;
}

inline ll Pre(){
  fac[0]=1;
  for (int i=1;i<MAXN;i++)
    fac[i]=fac[i-1]*i%P;
  Inv[1]=inv[1]=1;
  for (int i=2;i<MAXN;i++)
    Inv[i]=inv[i]=(P-P/i)*inv[P%i]%P;
  inv[0]=1;
  for (int i=1;i<MAXN;i++)
    inv[i]=inv[i-1]*inv[i]%P;
  for (int i=1;i<MAXN;i++)
    C[i]=Pow(2,(ll)i*(i-1)/2)*inv[i]%P;
}

ll N,n;  
ll L,w[2][MAXN],R[MAXN];  
  
inline void Init(int m){  
    L=0; for (n=1;n<=(m<<1);n<<=1) L++;    
    for (int i=0;i<n;i++) R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));    
    w[0][0]=w[0][n]=1;
    ll t=Pow(G,(P-1)/n);    
    for (int i=1;i<n;i++) w[0][i]=w[0][i-1]*t%P;    
    for (int i=0;i<n;i++) w[1][i]=w[0][n-i];  
}

inline void NTT(ll *a,int r){  
  for (int i=0;i<n;i++) if (i<R[i]) swap(a[i],a[R[i]]);   
  for (int i=1;i<n;i<<=1)      
    for (int j=0;j<n;j+=(i<<1))      
      for (int k=0;k<i;k++){     
	ll x=a[j+k],y=w[r][n/(i<<1)*k]*a[j+k+i]%P;      
	a[j+k]=x+y; a[j+k]=a[j+k]>P?a[j+k]-P:a[j+k]; 
	a[j+k+i]=x-y+P; a[j+k+i]=a[j+k+i]>P?a[j+k+i]-P:a[j+k+i];   
      }     
  if (r==1) for (int i=0,inv=Inv[n];i<n;i++) a[i]=a[i]*inv%P;    
}

ll f[MAXN];
ll A[MAXN<<1],B[MAXN<<1];  

inline void solve(int l,int r){
  if (l==r) return; int mid=(l+r)>>1,len=r-l;
  solve(l,mid);
  Init(len-1);
  for (int i=0;i<n;i++) A[i]=B[i]=0;  
  for (int i=l;i<=mid;i++) A[i-l]=f[i]*inv[i-1]%P;
  for (int i=0;i<len;i++) B[i]=C[i+1];
  NTT(A,0); NTT(B,0);
  for (int i=0;i<n;i++) A[i]=A[i]*B[i]%P;
  NTT(A,1);
  for (int i=mid+1;i<=r;i++)
    f[i]+=P-fac[i-1]*A[i-l-1]%P,f[i]=f[i]>P?f[i]-P:f[i];
  solve(mid+1,r);
}

int main(){
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  ll tem;
  Pre();
  read(tem);
  for (N=1;N<tem;N<<=1);
  for (int i=1;i<=N;i++)
    f[i]=Pow(2,(ll)i*(i-1)/2);
  solve(1,N);
  printf("%lld\n",f[tem]);
  return 0;
}