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
",(ll)a[m]*jie[m-1]%p);
87     return 0;
88 }
原文地址:https://www.cnblogs.com/ezyzy/p/6429785.html