bzoj 4451 : [Cerc2015]Frightful Formula FFT

4451: [Cerc2015]Frightful Formula

Time Limit: 10 Sec  Memory Limit: 64 MB
Submit: 177  Solved: 57
[Submit][Status][Discuss]

Description

给你一个n*n矩阵的第一行和第一列，其余的数通过如下公式推出： 

F[i,j]=a*f[i,j-1]+b*f[i-1,j]+c 

求f[n][n]%(10^6+3) 

Input

第一行三个数n，a，b，c 

第二行n个数，第i个表示f[i][1] 

第三行n个数，第i个表示f[1][i] 

Output

仅一个数表示f[n][n]%(10^6+3) 

Sample Input

Sample Input1:
3 0 0 0
0 0 2
0 3 0

Sample Input2:
4 3 5 2
7 1 4 3
7 4 4 8

Sample Output

Sample Output1:
0

Sample Output2:
41817
数据范围：
2<=n<=200000
其余的数大于等于0小于等于10^6

 

首先这道题每个变量对答案的贡献需要分开考虑。
第一行第i个数贡献为$a[1][i]\times a^{n-i} \times b^{n-1}\times C^{n-i}_{2n-i-2}$
相当于每个点先转到第二列对应点上，然后类似杨辉三角形向右下拓展，乘上对应的a和b。
列同理。
因为每个点会多产生出来一个c，所以还有一部分是c的贡献。
$$c\times \sum^{n}_{i=2} \sum^{n}_{j=2} a^{n-i}b^{n-j}C^{n-i}_{2n-i-j}=c\times \sum^{2n}_{i=4}(2n-i)!\sum_{j=2}^{n} \frac{a^{n-i+j}}{(n-i+j)!}\times \frac{b^{n-j}}{(n-j)!}$$
上FFT.
为了避免掉精度，把每个数拆成$a*1024+b$的两部分别卷积再合起来。

 

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120

#include<iostream> #include<cstdio> #include<algorithm> #include<cstring> #include<cmath> #define pi acos(-1) #define ll long long #define N 524289 #define NX 200005 using namespace std; struct E {     double x,y;     E (){;}     E (double _x,double _y){x=_x,y=_y;}     E operator + (const E a){return E(a.x+x,a.y+y);};     E operator - (const E a){return E(x-a.x,y-a.y);};     E operator * (const E a){return E(x*a.x-y*a.y,x*a.y+y*a.x);};   }a[N],b[N]; int n,R[N]; void FFT(E *a,int f) {     for(int i=0;i<n;i++)if(i<R[i])swap(a[i],a[R[i]]);     for(int i=1;i<n;i<<=1)     {         E wn(cos(pi/i),f*sin(pi/i));         for(int j=0;j<n;j+=(i<<1))         {             E w(1,0);             for(int k=0;k<i;k++,w=w*wn)             {                 E x=a[j+k],y=w*a[j+k+i];                 a[j+k]=x+y;a[j+k+i]=x-y;             }         }     }     if(f==-1)for(int i=0;i<n;i++)a[i].x/=n;     return ; } const int p = 1000003; const int M = 1024; ll pw(ll x,ll y) {     ll lst=1;     while(y)     {         if(y&1)lst=lst*x%p;         y>>=1;         x=x*x%p;        }     return lst; } int l[NX],r[NX],pwa[NX],pwb[NX],jie[2*NX],ni[NX]; ll A0[N],B0[N],A1[N],B1[N],T1[N],T2[N],T3[N],T4[N]; void mul(ll *ans,ll *a1,ll *b1) {     for(int i=0;i<n;i++)a[i].x=a[i].y=b[i].x=b[i].y=0;     for(int i=0;i<n;i++)a[i].x=a1[i];     for(int i=0;i<n;i++)b[i].x=b1[i];     FFT(a,1);FFT(b,1);     for(int i=0;i<n;i++)a[i]=a[i]*b[i];     FFT(a,-1);     for(int i=0;i<n;i++)ans[i]=((ll)(a[i].x+0.5))%p;     return ; } int main() {     int tn,ta,tb,tc;     scanf("%d%d%d%d",&tn,&ta,&tb,&tc);           for(int i=1;i<=tn;i++)scanf("%d",&l[i]);     for(int i=1;i<=tn;i++)scanf("%d",&r[i]);           jie[0]=1;ni[0]=ni[1]=1;     for(int i=1;i<=2*tn;i++)jie[i]=1LL*jie[i-1]*i%p;     for(int i=2;i<=tn;i++)ni[i]=(1LL*(-p/i)*ni[p%i]%p+p)%p;     for(int i=1;i<=tn;i++)ni[i]=1LL*ni[i]*ni[i-1]%p;           pwa[0]=1;pwb[0]=1;     for(int i=1;i<=tn;i++)pwa[i]=1LL*pwa[i-1]*ta%p;     for(int i=1;i<=tn;i++)pwb[i]=1LL*pwb[i-1]*tb%p;           for(int i=1;i<=tn;i++)T1[i]=1LL*pwa[tn-i]*ni[tn-i]%p;     for(int i=1;i<=tn;i++)T2[i]=1LL*pwb[tn-i]*ni[tn-i]%p;           ll ans=0;           for(int i=2;i<=tn;i++)ans+=1LL*r[i]*pwb[tn-1]%p*pwa[tn-i]%p*jie[2*tn-i-2]%p*ni[tn-2]%p*ni[tn-i]%p;     for(int i=2;i<=tn;i++)ans+=1LL*l[i]*pwa[tn-1]%p*pwb[tn-i]%p*jie[2*tn-i-2]%p*ni[tn-2]%p*ni[tn-i]%p;       int l=0;n=1;     while(n<=tn<<1)n<<=1,l++;     for(int i=0;i<n;i++)R[i]=(R[i>>1]>>1)|((i&1)<<(l-1));                 for(int i=2;i<=tn;i++)     {         A0[i]=T1[i]%M;A1[i]=T1[i]/M;         B0[i]=T2[i]%M;B1[i]=T2[i]/M;     }     memset(T1,0,sizeof(T1));memset(T2,0,sizeof(T2));     mul(T1,A0,B0);mul(T2,A1,B0);mul(T3,A0,B1);mul(T4,A1,B1);           for(int i=0;i<n;i++)     {         T2[i]+=T3[i];         (T1[i]+=T2[i]*M%p+T4[i]*M%p*M%p)%=p;     }           ll tmp=0;     for(int i=4;i<=2*tn;i++)     {         tmp+=1LL*jie[2*tn-i]*T1[i]%p;         tmp%=p;     }tmp=tmp*tc%p;           ans=(ans+tmp)%p;     printf("%lld\n",ans);     return 0; }