FFT与NTT

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-15

uoj#34

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 #define maxn 266666
 4 #define pi acos(-1.0)
 5 typedef complex<double> C ;
 6 C a[maxn],b[maxn];
 7 void fft(C *x,int n,int flag){
 8     if(n==1)return;
 9     C l[n>>1],r[n>>1];
10     for(int i=0;i<n;i+=2)
11         l[i>>1]=x[i],r[i>>1]=x[i+1];
12     fft(l,n>>1,flag),fft(r,n>>1,flag);
13     C wn(cos(2*pi/n),sin(flag*2*pi/n)),w(1,0);
14     for(int i=0;i<(n>>1);i++,w*=wn)
15         x[i]=l[i]+w*r[i],x[i+(n>>1)]=l[i]-w*r[i];
16 }
17 int main(){
18     int n,m;
19     scanf("%d%d",&n,&m);
20     for(int i=0;i<=n;i++)
21         scanf("%lf",&a[i].real());
22     for(int i=0;i<=m;i++)
23         scanf("%lf",&b[i].real());
24     m+=n;
25     for(n=1;n<=m;n<<=1);
26     fft(a,n,1),fft(b,n,1);
27     for(int i=0;i<=n;i++)
28         a[i]*=b[i];
29     fft(a,n,-1);
30     for(int i=0;i<=m;i++)
31         printf("%d ",(int)(a[i].real()/n+0.5));
32     return 0;
33 }

View Code

sdoi2015序列统计

 1 #include<bits/stdc++.h>
 2 using namespace std;
 3 typedef long long ll;
 4 #define MOD 1004535809
 5 #define G 3
 6 #define maxm 8005
 7 #define maxn 16400
 8 ll inv;
 9 int N,M,X,S,_n,sta[maxm],top,idx[maxm];
10 ll qp(ll bs,int x,int mod){
11     ll sum=1;
12     while(x){
13         if(x&1)sum=sum*bs%mod;
14         bs=bs*bs%mod;
15         x>>=1;
16     }
17     return sum;
18 }
19 void NTT(int *x,int n,int flag){
20     if(n==1)return;
21     int l[n>>1],r[n>>1];
22     for(int i=0;i<n;i+=2)
23         l[i>>1]=x[i],r[i>>1]=x[i+1];
24     NTT(l,n>>1,flag),NTT(r,n>>1,flag);
25     ll wn=qp(G,(ll)flag*(MOD-1)/n%(MOD-1),MOD),w=1;
26     for(int i=0;i<(n>>1);i++,w=w*wn%MOD)
27         x[i]=(l[i]+w*r[i]%MOD)%MOD,x[i+(n>>1)]=(l[i]-w*r[i]%MOD+MOD)%MOD;
28 }
29 struct node{
30     int a[maxn];
31     int& operator[](int x){
32         return a[x];
33     }
34     node operator*(node b){
35         NTT(a,_n,1),NTT(b.a,_n,1);
36         for(int i=0;i<_n;i++)
37             b[i]=(ll)b[i]*a[i]%MOD;
38         NTT(b.a,_n,MOD-2);
39         for(int i=M-1;i<=(M-2)<<1;i++)
40             b[i-(M-1)]=(b[i-(M-1)]+b[i])%MOD,b[i]=0;
41         for(int i=0;i<M-1;i++)
42             b[i]=b[i]*inv%MOD;
43         return b;
44     }
45 }sum,c;
46 int gtrt(){
47     int MM=M-1;
48     for(int i=2;i<=M-1;i++){
49         if(MM%i==0){
50             sta[++top]=i;
51             while(MM%i==0)MM/=i;
52         }
53     }
54     for(int i=2;;i++){
55         bool flag=true;
56         for(int j=1;j<=top;j++)
57             if(qp(i,(M-1)/sta[j],M)==1){
58                 flag=false;
59                 break; 
60             }
61         if(flag)return i;
62     }
63 }
64 node qp(node bs,int x){
65     sum[0]=1;
66     while(x){
67         if(x&1)sum=sum*bs;
68         bs=bs*bs;
69         x>>=1;
70     }
71     return sum;
72 }
73 int main(){
74     scanf("%d%d%d%d",&N,&M,&X,&S);
75     for(_n=1;_n<=M<<1;_n<<=1);
76     inv=qp(_n,MOD-2,MOD);
77     int g=gtrt();
78     for(int i=0,x=1;i<M-1;i++,x=x*g%M)idx[x]=i;
79     for(int i=1;i<=S;i++){
80         int x;
81         scanf("%d",&x);
82         if(x)c[idx[x]]=1;
83     }
84     node ans=qp(c,N);
85     printf("%d
",ans[idx[X]]);
86     return 0;
87 }

View Code