[BZOJ3456] [多项式求逆] 城市规划

#include <cstdio>
#include <iostream>
#include <algorithm>
#define N 390010
#define P 1004535809
#define G 3

using namespace std;

int n,m,num;
int w[2]
;
int a
,b
,c
,fac
;

inline int Pow(int x,int y){
int ret=1;
for(;y;y>>=1,x=1LL*x*x%P) if(y&1) ret=1LL*ret*x%P;
return ret;
}

inline void Pre(int n){
num=n;
int g=Pow(G,(P-1)/num),ig=Pow(g,P-2);
w[0][0]=w[1][0]=1;
for(int i=1;i<num;i++)
w[0][i]=1LL*w[0][i-1]*ig%P,w[1][i]=1LL*w[1][i-1]*g%P;
}

int tmp
,rev
,invb
;

inline void FFT(int *a,int n,int r){
for(int i=0;i<n;i++) if(rev[i]>i) swap(a[i],a[rev[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++){
int x=a[j+k],y=1LL*w[r][num/(i<<1)*k]*a[i+j+k]%P;
a[j+k]=(x+y)%P; a[j+k+i]=(x+P-y)%P;
}
if(!r) for(int i=0;i<n;i++) a[i]=1LL*a[i]*Pow(n,P-2)%P;
}

void Inv(int *a,int *b,int n){
if(n==1){
b[0]=Pow(a[0],P-2); return ;
}
Inv(a,b,n>>1);
for(int i=0;i<n;i++) tmp[i]=a[i],tmp[i+n]=0;
int L=0; while(!(n>>L&1)) L++;
for(int i=1;i<(n<<1);i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<L);
FFT(tmp,n<<1,1); FFT(b,n<<1,1);
for(int i=0;i<(n<<1);i++)
tmp[i]=(1LL*b[i]*2+P-1LL*tmp[i]*b[i]%P*b[i]%P)%P;
FFT(tmp,n<<1,0);
for(int i=0;i<n;i++) b[i]=tmp[i],b[i+n]=0;
}

int main(){
scanf("%d",&n);
for(m=1;m<=n;m<<=1); Pre(m<<1);
fac[0]=1;
for(int i=1;i<=n;i++) fac[i]=1LL*fac[i-1]*i%P;
for(int i=1;i<=n;i++)
c[i]=1LL*Pow(2,(1LL*i*(i-1)>>1)%(P-1))*Pow(fac[i-1],P-2)%P;
for(int i=0;i<=n;i++)
b[i]=1LL*Pow(2,(1LL*i*(i-1)>>1)%(P-1))*Pow(fac[i],P-2)%P;
Inv(b,invb,m);
FFT(c,m<<1,1); FFT(invb,m<<1,1);
for(int i=0;i<m<<1;i++) a[i]=1LL*c[i]*invb[i]%P;
FFT(a,m<<1,0);
printf("%lld\n",1ll*a
*fac[n-1]%P);
return 0;
}