BZOJ3456：城市规划——题解

#include<cstdio>
#include<cctype>
#include<cstring>
#include<vector>
#include<cmath>
#include<algorithm>
#include<iostream>
using namespace std;
typedef long long ll;
const ll P=1004535809;
const int G=3;
const int N=5e5+5;
ll qpow(ll a,ll n,ll p){
    ll res=1;
    while(n){
        if(n&1)res=res*a%p;
        a=a*a%p;n>>=1;
    }
    return res;
}
void MTT(ll a[],int n,int on){
    for(int i=1,j=n>>1;i<n-1;i++){
        if(i<j)swap(a[i],a[j]);
        int k=n>>1;
        while(j>=k){j-=k;k>>=1;}
        if(j<k)j+=k;
    }
    for(int i=2;i<=n;i<<=1){
        ll res=qpow(G,(P-1)/i,P);
        for(int j=0;j<n;j+=i){
            ll w=1;
            for(int k=j;k<j+i/2;k++){
                ll u=a[k],t=w*a[k+i/2]%P;
                a[k]=(u+t)%P;
                a[k+i/2]=(u-t+P)%P;
                w=w*res%P;
            }
        }
    }
    if(on==-1){
        ll inv=qpow(n,P-2,P);
        a[0]=a[0]*inv%P;
        for(int i=1;i<=n/2;i++){
            a[i]=a[i]*inv%P;
            if(i!=n-i)a[n-i]=a[n-i]*inv%P;
            swap(a[i],a[n-i]);
        }
    }
}
ll t[N];
void inv(int deg,ll a[],ll b[]){
    if(deg==1){
        b[0]=qpow(a[0],P-2,P);
        return;
    }
    inv((deg+1)>>1,a,b);
    int n=1;
    while(n<(deg<<1))n<<=1;
    for(int i=0;i<deg;i++)t[i]=a[i];
    for(int i=deg;i<n;i++)t[i]=0;
    MTT(t,n,1);MTT(b,n,1);
    for(int i=0;i<n;i++)
        b[i]=b[i]*(2-b[i]*t[i]%P+P)%P;
    MTT(b,n,-1);
    for(int i=deg;i<n;i++)b[i]=0;
}
int n;
ll f[N],g[N],tmp[N],c[N],jc[N],C[N][3];
inline void init(){
    jc[0]=jc[1]=1;C[0][2];
    for(int i=2;i<=n;i++)jc[i]=jc[i-1]*i%P;
    for(int i=1;i<=n;i++)C[i][2]=C[i-1][2]+i-1;
}
int main(){
    scanf("%d",&n);init();
    for(int i=0;i<=n;i++){
        g[i]=qpow(2,C[i][2],P)*qpow(jc[i],P-2,P)%P;
        if(i)c[i]=qpow(2,C[i][2],P)*qpow(jc[i-1],P-2,P)%P;
    }
    int nn=1;
    while(nn<=n)nn<<=1;
    
    inv(nn,g,tmp);memcpy(g,tmp,sizeof(tmp));
    
    MTT(g,nn,1);MTT(c,nn,1);
    for(int i=0;i<nn;i++)f[i]=g[i]*c[i]%P;
    MTT(f,nn,-1);
    
    f[n]=f[n]*jc[n-1]%P;
    printf("%lld\n",f[n]);
    return 0;
}