51nod1172 Partial Sums V2

推一下式子发现是裸的FFT,$ans[k]=\sum_{i}\sum_{j}[i+j=k]a[i]*C_{m-1+j}^{j}$

比较坑爹的就是这个模数,于是我们上任意模数的FFT

任意模数的FFT目的就是降低卷积中的元素上界,我们设$P=\lfloor \sqrt{mod} \rfloor$,

我们将原函数中的系数变成两个$a1[i]=a[i]/P,a2[i]=a[i]%P,b1[i]=b[i]/P,b2[i]=b[i]%P$

这样我们新卷积出来的值的上限是$n*mod$,

$P^2$的系数是$a1*b1$,$P$的系数是$a1*b2+a2*b1$,$1$的系数是$a2*b2$,然后我们再分别卷积,最后统计答案即可

一开始炸精了。。改成预处理单位复数根就好了。

 1 #include <cstdio>
 2 #include <cstring>
 3 #include <iostream>
 4 #include <algorithm>
 5 #include <cmath>
 6 #define mod 1000000007
 7 #define MAXN 131072
 8 using namespace std;
 9 const double pi=acos(-1.0);
10 const int P=31622;
11 struct cmx{
12     double x,y;
13     cmx(){}
14     cmx(double a,double b){x=a,y=b;}
15     cmx operator + (cmx a){return cmx(x+a.x,y+a.y);}
16     cmx operator - (cmx a){return cmx(x-a.x,y-a.y);}
17     cmx operator * (cmx a){return cmx(x*a.x-y*a.y,x*a.y+y*a.x);}
18     cmx operator / (double a){return cmx(x/a,y/a);}
19 }A[MAXN],B[MAXN],C[MAXN],D[MAXN],E[MAXN],F[MAXN],G[MAXN],wn[MAXN],inv[MAXN];
20 int n,m,N,rev[MAXN];
21 long long ans[MAXN],a[MAXN],b[MAXN];
22 void fft(cmx *a,int o, cmx *wn){
23     register int i,j,k;
24     cmx t;
25     for(i=0;i<N;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
26     for(k=2;k<=N;k<<=1)
27         for(j=0;j<N;j+=k)
28             for(i=0;i<(k>>1);i++){
29                 t=wn[N/k*i]*a[i+j+(k>>1)];
30                 a[i+j+(k>>1)]=a[i+j]-t;
31                 a[i+j]=a[i+j]+t;
32             }
33     if(o==-1)for(i=0;i<N;i++)a[i]=a[i]/N;
34 }
35 void FFT(long long *a,long long *b,long long *ans){
36     register int i;
37     for(i=0;i<N;i++){
38         A[i]=cmx(a[i]/P,0);
39         B[i]=cmx(a[i]%P,0);
40         C[i]=cmx(b[i]/P,0);
41         D[i]=cmx(b[i]%P,0);
42     }
43     fft(A,1,wn);fft(B,1,wn);fft(C,1,wn);fft(D,1,wn);
44     for(i=0;i<N;i++){
45         E[i]=A[i]*C[i];
46         F[i]=A[i]*D[i]+B[i]*C[i];
47         G[i]=B[i]*D[i];
48     }
49     fft(E,-1,inv);fft(F,-1,inv);fft(G,-1,inv);
50     register long long tmp;
51     for(i=0;i<n;i++){
52         tmp=1ll*round(E[i].x);
53         ans[i]=tmp%mod*P%mod*P%mod;
54         tmp=1ll*round(F[i].x);
55         (ans[i]+=tmp%mod*P%mod)%=mod;
56         (ans[i]+=1ll*round(G[i].x))%=mod;
57     }
58 }
59 long long qp(long long a,int b){
60     long long c=1;
61     while(b){
62         if(b&1)c=c*a%mod;
63         a=a*a%mod; b>>=1;
64     }return c;
65 }
66 int main(){
67     scanf("%d%d",&n,&m);
68     register int i;
69     for(i=0;i<n;i++)scanf("%lld",&a[i]);
70     b[0]=1;
71     for(i=1;i<n;i++)
72         b[i]=1ll*b[i-1]*(m-1+i)%mod*qp(i,mod-2)%mod;
73     for(N=1;N<(n<<1);N<<=1);
74     for(i=0;i<N;i++){
75         if(i&1)rev[i]=(rev[i>>1]>>1)|(N>>1);
76         else rev[i]=rev[i>>1]>>1;
77     }
78     for(i=0;i<N;i++)
79         wn[i]=cmx(cos(2*pi/N*i),sin(2*pi/N*i)),
80         inv[i]=cmx(cos(2*pi/N*i),-sin(2*pi/N*i));
81     FFT(a,b,ans);
82     for(i=0;i<n;i++)printf("%lld\n",ans[i]);
83     return 0;
84 }
View Code

 

posted @ 2018-02-27 08:09  Ren_Ivan  阅读(110)  评论(0编辑  收藏