FFT的常数优化

卡得一手好常数。。学习了。。(似乎只对FFT有效)

JZOJ 4349

 

 1 #include <bits/stdc++.h>
 2 #define LL long long
 3 #define DB long double
 4 using namespace std;
 5 const int mo=100003;
 6 const DB pi=acos(-1);
 7 int K,T,p[360005],q[360005],f[70005],n,m,k,a[70005],rev[70005],ans;
 8 DB COS[70005],SIN[70005],Cos[70005],Sin[70005];
 9 struct cp{
10     DB x,y;
11     cp (DB p=0,DB q=0){x=p,y=q;}
12 }w[70005],N[35005];
13 cp operator+(cp a,cp b){
14     return cp(a.x+b.x,a.y+b.y);
15 }
16 cp operator-(cp a,cp b){
17     return cp(a.x-b.x,a.y-b.y);
18 }
19 cp operator*(cp a,cp b){
20     return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
21 }
22 LL po(LL x,int y){
23     LL z=1;
24     for (;y;y>>=1,x=x*x%mo)
25     if (y&1) z=z*x%mo;
26     return z;
27 }
28 LL C(int n,int m){
29     return n<m?0:(n<mo?1ll*p[n]*q[m]%mo*q[n-m]%mo:C(n%mo,m%mo)*C(n/mo,m/mo)%mo);
30 }
31 void fft(cp a[],int n,int d=1){
32     for (int i=0;i<n;++i) if (i<rev[i]) swap(a[i],a[rev[i]]);
33     for (int m=2,k=1;m<=n;k=m,m<<=1){
34         cp wn;
35         wn=d>0?cp(COS[m],SIN[m]):cp(Cos[m],Sin[m]);
36         N[0]=cp(1,0);
37         for (int i=1;i<k;++i) N[i]=N[i-1]*wn;
38         for (int i=0;i<n;i+=m)
39         for (int j=i;j<i+k;++j){
40             cp u=a[j],v=a[j+k]*N[j-i];
41             a[j]=u+v; a[j+k]=u-v;
42         }
43     }
44     if (d==-1) for (int i=0;i<=n;++i) a[i].x/=n;
45 }
46 void preF(){
47     f[1]=9; f[2]=243; f[3]=16224;
48     int t=12,X=299,Y=1<<t;
49     for (int i=4;i<=60000;++i){
50         int x=1,y=0;
51         for (int j=0;j<=4;++j,x*=3)
52             (f[i]+=(C(4,j)*x*(Y-X+y)%mo))%=mo,(y+=C(t,i-j-1))%=mo;
53         f[i]=(f[i]%mo+mo)%mo;
54         for (int j=1;j<=6;++j)
55             X=(X*2-C(t,i-1))%mo,++t;
56         X=(X+C(t,i))%mo; Y=Y*64%mo;
57     }
58 }
59 int main(){
60     for (int i=2;i<=65536;i<<=1)
61         COS[i]=cos(pi*2/i),Cos[i]=cos(pi*-2/i),
62         SIN[i]=sin(pi*2/i),Sin[i]=sin(pi*-2/i);
63     scanf("%d",&T);
64     p[0]=q[0]=1;
65     for (int i=1,x;i<mo;++i){
66         for (x=i;x%mo==0;x/=mo);
67         p[i]=1ll*p[i-1]*x%mo;
68     }
69     q[mo-1]=po(p[mo-1],mo-2);
70     for (int i=mo-2,x;i;--i){
71         for (x=i+1;x%mo==0;x/=mo);
72         q[i]=1ll*q[i+1]*x%mo;
73     }
74     preF();
75     while (T--){
76         scanf("%d%d",&n,&K);
77         for (int i=1;i<=n;++i) a[i]=f[i];
78         for (m=1;m<n;m<<=1);
79         for (int i=n+1;i<=m;++i) a[i]=0; n=m;
80         for (m=2,k=1;m<=n;k=m,m<<=1){
81             for (int i=0;i<m;++i) rev[i]=rev[i>>1]>>1|(i&1)*(m>>1);
82             for (int i=1;i<=n;i+=m){
83                 w[0]=cp(2,0);
84                 for (int j=i;j<i+k;++j) w[j-i+1]=cp(a[j]+a[j+k],a[j]-a[j+k]);
85                 for (int j=k+1;j<m;++j) w[j]=cp(0,0);
86                 fft(w,m);
87                 for (int j=0;j<m;++j) w[j]=w[j]*w[j];
88                 fft(w,m,-1);
89                 for (int j=i;j<i+m-1;++j) a[j]=(LL)round(w[j-i+1].x/4)%mo;
90                 a[i+m-1]=(LL)round(w[0].x/4-1)%mo;
91             }
92         }
93         ans=0;
94         for (int i=K;i<=n;++i) (ans+=a[i])%=mo;
95         printf("%d\n",(ans%mo+mo)%mo);
96     }
97     return 0;
98 }
DAISHI

 

posted @ 2018-01-07 11:58  cyz666  阅读(494)  评论(0编辑  收藏  举报