记录编号	372715	评测结果	AAAAAAAAAA
题目名称	[BZOJ 3456] 城市规划	最终得分	100
用户昵称	stdafx.h	是否通过	通过
代码语言	C++	运行时间	4.104 s
提交时间	2017-02-18 22:07:42	内存使用	10.59 MiB

/*
    G(x) = sigma ( 2^C(i,2) x ^ i / ( i! ) )
    F(x) = sigma ( fi * x ^ i / ( i! ) )
    G(x) = e ^ F(x)
    F(x) = ln ( G(x) )
*/
 
#define fastcall __attribute__((optimize("-Os")))
#define IL __inline__ __attribute__((always_inline))
 
#define maxn 270010ul
 
#define ll long long
 
#include <stdio.h>
#include <string.h>
#include <algorithm>
 
const int G=3,P=1004535809;

int R;

fastcall IL int mul_mod(int a,int b){
    __asm__ __volatile__ ("\tmull %%ebx\n\tdivl %%ecx\n" :"=d"(R):"a"(a),"b"(b),"c"(P));
    return R;
}

fastcall IL int sub(int x,int y){
    x-=y; return x<0?x+P:x;
}

namespace polynomial{
	int W[maxn];
    fastcall IL int ksm(int p,ll k){
        int ans=1;
        while(k){
            if(k&1) ans=mul_mod(ans,p);
            k>>=1ll,p=mul_mod(p,p);
        } return ans;
    }
    fastcall void ntt(int *a,int n,int flag){
        for(int i=0,j=0;i<n;i++){
            if(i>j) std::swap(a[i],a[j]);
            for(int t=n>>1;(j^=t)<t;t>>=1);
        } register int wn,u,t;
        for(register int m=2,i,k,p;m<=n;m<<=1){
            wn=ksm(G,flag==1?(P-1)/m:(P-1-(P-1)/m));
            W[0]=1; for(i=1,p=m>>1;i<p;i++) W[i]=mul_mod(W[i-1],wn);
            for(i=0;i<n;i+=m) for(k=0;k<p;k++)
                t=1ll*W[k]*a[i+k+p]%P,u=a[i+k],a[i+k]=(u+t)%P,a[i+k+p]=sub(u,t);
        }
        if(flag==-1){
            int inv=ksm(n,P-2);
            for(int i=0;i<n;i++) a[i]=1ll*a[i]*inv%P;
        } return;
    }
    void poly_derivative(int *poly,int *der,int n){
        poly[n]=0;
        for(int i=0;i<n;i++) der[i]=1ll*(i+1)*poly[i+1]%P;
        return;
    }
    void poly_integrate(int *poly,int n){
        static int tmp[maxn];
        for(int i=0;i<n;i++) tmp[i+1]=1ll*poly[i]*ksm(i+1,P-2)%P;
        for(int i=0;i<n;i++) poly[i]=tmp[i];
        return;
    }
    /*
 
    & polynomial inverse
 
    -> G(x) * F(x) = 1 ( mod x^n )
    -> H(x) * F(x) = 1 ( mod x^2n )
 
    -> we want to know -> H(x)
 
    -> G(x) * F(x) - 1 =0 ( mod x^n )
    -> [G(x) * F(x) - 1] ^ 2 = 0 ( mod x^2n )
    -> [G(x)F(x)] ^ 2 - 2 * G(x) * F(x) + 1 = 0 ( mod x^2n )
    -> G(x) ^ 2 * F(x) - 2 * G(x) + H(x) = 0 ( mod x^2n )
    -> H(x) = 2 * G(x) - G(x) ^ 2 * F(x)
 
    -> H(x) = G(x) * [ 2 - G(x) * F(x) ]
 
    */
    void poly_inverse(int *f,int *g,int n){
        static int tmp[maxn];
        if(n==1){g[0]=ksm(f[0],P-2);return;}
        poly_inverse(f,g,n>>1);
        memcpy(tmp,f,sizeof(f[0])*n);
        memset(tmp+n,0,sizeof(f[0])*n);
        ntt(tmp,n<<1,1),ntt(g,n<<1,1);
        for(int i=0;i<(n<<1);i++) tmp[i]=1ll*g[i]*(2ll-1ll*g[i]*tmp[i]%P)%P;
        ntt(tmp,n<<1,-1);
        for(int i=0;i<n;i++) g[i]=tmp[i];
        memset(g+n,0,sizeof(g[0])*n);
        return;
    }
    /*
 
    & polynomial sqrt
 
    ->   G(x) ^ 2 = F(x) ( mod x^n )
    ->   H(x) ^ 2 = F(x) ( mod x^2n )
 
    ->   we want to know -> H(x)
 
    ->   G(x) ^ 2 - F(x) = 0 ( mod x^n )
    ->   [ G(x) ^ 2 - F(x) ] ^ 2 = 0 ( mod x^2n )
    ->   [ G(x) ^ 2 + F(x) ] ^ 2 = 4 * ( G(x) ^ 2 ) * F(x) ( mod x^2n )
    ->   [ (G(x) ^ 2 + F(x)) / ( 2 * G(x) ) ] ^ 2 = F(x) ( mod x^2n )
 
    ->   H(x) = (G(x) ^ 2 + F(x)) / ( 2 * G(x) )
 
    */
    void poly_sqrt(int *f,int *g,int n){
        static int tmp[maxn],ginv[maxn];
        if(n==1){g[0]=1;return;}
        int inv2=ksm(2,P-2);
        poly_sqrt(f,g,n>>1);
        memset(ginv,0,sizeof(ginv[0])*n);
        poly_inverse(g,ginv,n);
        memcpy(tmp,f,sizeof(f[0])*n);
        memset(tmp+n,0,sizeof(f[0])*n);
        ntt(tmp,n<<1,1),ntt(g,n<<1,1),ntt(ginv,n<<1,1);
        for(int i=0;i<(n<<1);i++) tmp[i]=1ll*inv2*(g[i]+1ll*tmp[i]*ginv[i]%P)%P;
        ntt(tmp,n<<1,-1);
        for(int i=0;i<n;i++) g[i]=tmp[i];
        memset(g+n,0,sizeof(g[0])*n);
        return;
    }
    /*
 
    & polynomial ln
    ->   F(x) = ln ( G(x) )
    ->   F'(x) = G'(x) / G(x)
 
    ->   F(x) = ∫( G'(x) / G(x) ) dx
 
    */
    void poly_ln(int *poly,int *ln,int n){
        static int der[maxn],inv[maxn];
        memset(inv,0,sizeof(inv));
        poly_inverse(poly,inv,n);
        poly_derivative(poly,der,n);
        ntt(der,n<<1,1),ntt(inv,n<<1,1);
        for(int i=0;i<(n<<1);i++)
            ln[i]=1ll*der[i]*inv[i]%P;
        ntt(ln,n<<1,-1);
        memset(ln+n,0,sizeof(ln[0])*n);
        poly_integrate(ln,n);
        return;
    }
};
 
int n,m,g[maxn],f[maxn],fac[maxn];
 
inline ll C(ll i,ll j=2){
    return i*(i-1)>>1ll;
}
 
int main(){
	freopen("bzoj_3456.in","r",stdin);
	freopen("bzoj_3456.out","w",stdout);
    scanf("%d",&n);
    for(m=1;m<=(n<<1);m<<=1); m>>=1; 
    fac[0]=1;for(int i=1;i<m;i++) fac[i]=1ll*fac[i-1]*i%P;
    for(int i=0;i<m;i++) g[i]=1ll*polynomial::ksm(2,C(i,2))*polynomial::ksm(fac[i],P-2)%P;
    polynomial::poly_ln(g,f,m);
    f[n]=1ll*f[n]*fac[n]%P,(f[n]+=P)%=P,printf("%d",f[n]);
    return 0;
}