NTT+多项式求逆
1 #include<iostream> 2 #include<cstdio> 3 #include<algorithm> 4 #include<cstring> 5 #include<cmath> 6 #define ll long long 7 #define N 300005 8 using namespace std; 9 const int ni = 3; 10 const int p = 1004535809; 11 ll pw(ll x,int y) 12 { 13 ll lst=1; 14 while(y) 15 { 16 if(y&1)lst=(lst*x)%p; 17 y>>=1; 18 x=(x*x)%p; 19 } 20 return lst; 21 } 22 int n,m; 23 int c[N],a[N],b[N],R[N]; 24 int invb[N]; 25 int ji[N],jie[N],nini[N]; 26 int tmp[N]; 27 void fft(int *a,int n,int f) 28 { 29 for(int i=0;i<n;i++)if(R[i]>i)swap(a[i],a[R[i]]); 30 for(int i=1;i<n;i<<=1) 31 { 32 int wn=pw(ni,((p-1)/(i<<1)*f+p-1)%(p-1)); 33 for(int j=0;j<n;j+=(i<<1)) 34 { 35 int w=1; 36 for(int k=0;k<i;k++,w=((ll)w*wn)%p) 37 { 38 int x=a[j+k];int y=((ll)a[j+i+k]*w)%p; 39 a[j+k]=(x+y)%p;a[j+k+i]=(x-y+p)%p; 40 } 41 } 42 } 43 if(f==-1) 44 { 45 ll nii=pw(n,p-2); 46 for(int i=0;i<n;i++)a[i]=(ll)a[i]*nii%p; 47 } 48 return ; 49 } 50 void get_inv(int *a,int *b,int n) 51 { 52 if(n==1) 53 { 54 b[0]=pw(a[0],p-2); 55 return ; 56 } 57 get_inv(a,b,n>>1); 58 for(int i=0;i<n;i++)tmp[i]=a[i]; 59 int l=0;int nn=1; 60 while(nn<n<<1)nn<<=1,l++; 61 for(int i=0;i<nn;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(l-1)); 62 fft(tmp,nn,1);fft(b,nn,1); 63 for(int i=0;i<nn;i++) 64 b[i]=((ll)b[i]*(2-(ll)tmp[i]*b[i]%p+p))%p; 65 fft(b,nn,-1); 66 memset(b+n,0,sizeof(int)*n); 67 } 68 int main() 69 { 70 scanf("%d",&n); 71 int l=0; 72 for(m=n,n=1;n<=m;n<<=1)l++; 73 nini[0]=jie[0]=ji[0]=1; 74 for(int i=1;i<n;i++)jie[i]=((ll)jie[i-1]*i)%p; 75 for(int i=0;i<n;i++) 76 b[i]=pw(2,((ll)i*(i-1)>>1)%(p-1))*pw(jie[i],p-2)%p; 77 for(int i=1;i<n;i++) 78 c[i]=pw(2,((ll)i*(i-1)>>1)%(p-1))*pw(jie[i-1],p-2)%p; 79 get_inv(b,invb,n); 80 l++;int tp=n<<1; 81 for(int i=0;i<tp;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(l-1)); 82 fft(invb,tp,1); 83 fft(c,tp,1); 84 for(int i=0;i<tp;i++)a[i]=(ll)invb[i]*c[i]%p; 85 fft(a,tp,-1); 86 printf("%d\n",(ll)a[m]*jie[m-1]%p); 87 return 0; 88 }