Frightful Formula（杨辉三角求和 组合数学）

original link - https://www.luogu.org/problem/P4351

题意：

一个 $n*n$矩阵，给出第一行的数 $x_i$和第一列的数 $y_i$，其他的数 $F(i,j)=aF(i,j-1)+bF(j,i-1)+c$。求 $F(n,n)$。

解析：

简答推几项可以发现 $F$中的 $x_i,y_i,c$是独立的，所以考虑怎么快速求出每一项的系数矩阵。

假设左上角为 $1$，那么推出的每个位置的贡献应该是这样的：

$\begin{bmatrix} 1&a&a^2&a^3\\b&2ab&3a^2b&4a^3b\\b^2&3ab^2&6a^2b^2&10a^3b^2\\b^3&4ab^3&10a^2b^3&20a^3b^3 \end{bmatrix}$

可以发现是杨辉三角乘上 $a^ib^j$，坑的在于 $F(1,i)$对 $F(1,i+1)$并没有贡献，所以要将 $F(1,i)*b$作为 $1$，从 $(2,i)$开始。

这是比较简单的部分，还有 $c$的系数呢， $c$从 $F(2,2)$开始，贡献应该是：

$\begin{bmatrix} 1&a+1\\b+1&a(b+1)+b(a+1)+1 \end{bmatrix}$

这样不太好考虑，我们发现和原来的区别在于每一个位置多了一个 $1$，所以自己推一下可以发现就是原来的那个矩阵求和：

$\begin{bmatrix} 1&a&a^2&a^3\\b&2ab&3a^2b&4a^3b\\b^2&3ab^2&6a^2b^2&10a^3b^2\\b^3&4ab^3&10a^2b^3&20a^3b^3 \end{bmatrix}$

考虑斜着的每一列，发现副对角线以上部分，每一列的和都可以从上一列的和乘上 $(a+b)$得到，而过了副对角线后要减掉两边那两个的贡献。以第 $4$列（ $b^3,3ab^2,3a^2b,a^3$）为例，往下的 $b^3*b$已经超出矩阵了，所以要减掉。

代码：

/*
 *  Author : Jk_Chen
 *    Date : 2019-09-22-14.30.03
 */
#include<bits/stdc++.h>
using namespace std;
#define LL long long
#define rep(i,a,b) for(int i=(int)(a);i<=(int)(b);i++)
#define per(i,a,b) for(int i=(int)(a);i>=(int)(b);i--)
#define mmm(a,b) memset(a,b,sizeof(a))
#define pb push_back
#define pill pair<int, int>
#define fi first
#define se second
#define debug(x) cerr<<#x<<" = "<<x<<'\n';
const LL mod=1e6+3;
const int maxn=1e6+2;
LL rd(){ LL ans=0; char last=' ',ch=getchar();
    while(!(ch>='0' && ch<='9'))last=ch,ch=getchar();
    while(ch>='0' && ch<='9')ans=ans*10+ch-'0',ch=getchar();
    if(last=='-')ans=-ans; return ans;
}
/*_________________________________________________________begin*/

LL Pow(LL a,LL b,LL mod){
    LL res=1;
    while(b>0){

        if(b&1)res=res*a%mod;

        a=a*a%mod;
        b>>=1;

    }
    return res;
}
/*_________________________________________________________Pow*/

LL fac[maxn];
LL ifac[maxn];
LL C(int i,int j){return fac[i]*ifac[j]%mod*ifac[i-j]%mod; }

LL n,a,b,c;

LL Pa[maxn],Pb[maxn];
LL Get(int i,int j){

    return C(i+j-2,i-1)*Pa[j-1]%mod*Pb[i-1]%mod;

}
LL Row[maxn],Column[maxn];
int main(){
    fac[0]=1;
    rep(i,1,maxn-1)fac[i]=fac[i-1]*i%mod;
    ifac[maxn-1]=Pow(fac[maxn-1],mod-2,mod);
    per(i,maxn-2,0){
        ifac[i]=ifac[i+1]*(i+1)%mod;
    }

    n=rd(),a=rd(),b=rd(),c=rd();

    Pa[0]=Pb[0]=1;
    rep(i,1,maxn-1)Pa[i]=Pa[i-1]*a%mod,Pb[i]=Pb[i-1]*b%mod;



    rep(i,1,n){
        Column[i]=rd();
    }
    rep(i,1,n){
        Row[i]=rd();
    }

    LL ans=0;
    rep(i,2,n){
        ans=(ans+Row[i]*b%mod*Get(n-1,n-i+1)%mod)%mod;

    }
    rep(i,2,n){
        ans=(ans+Column[i]*a%mod*Get(n-i+1,n-1)%mod)%mod;
    }

    LL sum=1;
    n--;
    LL tmp=1;
    rep(i,2,n){
        tmp=tmp*(a+b)%mod;

        sum=(sum+tmp)%mod;
    }
    rep(i,2,n){
        tmp=(tmp*(a+b)-Get(n,i-1)*b-Get(i-1,n)*a)% mod;
        sum=(sum+tmp)%mod;
    }
    sum=sum*c%mod;


    ans%=mod;
    ans=((ans+sum)%mod+mod)%mod;
    printf("%lld\n",ans);
    return 0;
}

/*_________________________________________________________end*/