BZOJ 4451: [Cerc2015]Frightful Formula

Description

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

\(f_{i,j}=a\times f_{i,j-1}+b\times f_{i-1,j}+c\)

求\(f_{n,n} \mod (10^6+3)\)

Solution

递推/FFT.

写出来式子，一顿胡推..

当时化完没写博客...现在不想化了...

Code

/**************************************************************
    Problem: 4451
    User: BeiYu
    Language: C++
    Result: Accepted
    Time:752 ms
    Memory:29424 kb
****************************************************************/
 
#include <bits/stdc++.h>
using namespace std;
 
typedef long long LL;
const int N = 400050;
const LL p = 1e6+3;
 
inline LL in(LL x=0,char ch=getchar()) { while(ch>'9'||ch<'0') ch=getchar();
    while(ch>='0'&&ch<='9')x=x*10+ch-'0',ch=getchar();return x; }
LL Add(LL &x,LL y) { x=(x+y)%p; }
LL Pow(LL a,LL b,LL r=1) { for(;b;b>>=1,a=a*a%p) if(b&1) r=r*a%p;return r; }
 
LL n,m,a,b,c,ans;
LL t1[N],t2[N];
LL fac[N],ifac[N];
LL pw_a[N],pw_b[N];
LL f[N],g_a[N],g_b[N];
 
LL C(LL n,LL m) { return n<m?0:fac[n]*ifac[m]%p*ifac[n-m]%p; }
LL invv(LL x) { return x%p==0?1:Pow(x%p,p-2); }
void init() {
    m=n<<1;
    fac[0]=1;
    for(int i=1;i<=m;i++) fac[i]=fac[i-1]*i%p;
    ifac[m]=invv(fac[m]);
    for(int i=m-1;~i;--i) ifac[i]=ifac[i+1]*(i+1)%p;
    ifac[0]=1;
    pw_a[0]=1;for(int i=1;i<=n+1;i++) pw_a[i]=pw_a[i-1]*a%p;
    pw_b[0]=1;for(int i=1;i<=n+1;i++) pw_b[i]=pw_b[i-1]*b%p;
}
void work1() {
    for(int i=2;i<=n;i++) Add(ans,C(n-2+n-i,n-i)*pw_a[n-1]%p*pw_b[n-i]%p*t1[i]%p);
    for(int i=2;i<=n;i++) Add(ans,C(n-2+n-i,n-i)*pw_b[n-1]%p*pw_a[n-i]%p*t2[i]%p);
}
void work2() {
    LL inv_a=invv(1-a+p),inv_b=invv(1-b+p);
    if(a==0) { for(int i=0;i<=n-2;i++) Add(ans,pw_b[i]*c);return; }
    if(b==0) { for(int i=0;i<=n-2;i++) Add(ans,pw_a[i]*c);return;    }
    if(a==1) for(int i=0;i<=n-2;i++) g_a[i]=C(2*i+1,i);
    else { g_a[0]=1;for(int i=1;i<=n;i++) g_a[i]=(g_a[i-1]+C(2*i-1,i)*pw_a[i]%p-C(2*i,i)*pw_a[i+1]%p+p)%p*inv_a%p; }
    if(b==1) { for(int i=0;i<=n-2;i++) g_b[i]=C(2*i+1,i);    }
    else { g_b[0]=1;for(int i=1;i<=n;i++) g_b[i]=(g_b[i-1]+C(2*i-1,i)*pw_b[i]%p-C(2*i,i)*pw_b[i+1]%p+p)%p*inv_b%p; }
    f[0]=1;for(int i=1;i<=n;i++) f[i]=(f[i-1]+pw_a[i]*g_b[i]%p+pw_b[i]*g_a[i]%p-pw_a[i]*pw_b[i]*C(2*i,i)%p+p)%p;
    Add(ans,f[n-2]*c%p);
}
int main() {
    n=in(),a=in(),b=in(),c=in();
    for(int i=1;i<=n;i++) t1[i]=in();
    for(int i=1;i<=n;i++) t2[i]=in();
    init();
    work1();
    work2();
    printf("%lld\n",ans);
    return 0;
}