正题

题目链接:https://www.luogu.com.cn/problem/P4351

---

题目大意

$n\ast n$的矩形，给出第一行和第一列的数，剩下的满足$F_{i,j}=a\ast F_{i,j-1}+b\ast F_{i-1,j}+c$

求$F_{n,n}$

---

解题思路

第一眼看以为是水题，因为给出的数字的贡献通过组合数很好算，但是后来发现麻烦的是那个$c$。我们考虑每个格子的$c$产生的贡献。

下面为了方便我们先默认让所有格子横纵坐标减$1$
对于一个格子$(i,j)$，通过它的路径有$\left(\genfrac{}{}{0ex}{}{2n-i-j}{n-i}\right)$种，然后产生的贡献是$a^{n-i}{b}^{n-j}$。为了方便我们反过来表示，然后因为第一行第一列没有贡献所以$n$减一。
那么总共的贡献就是$$\sum _{i=0}^n\sum _{j=0}^n{a}^i{b}^i\left(\genfrac{}{}{0ex}{}{i+j}{i}\right)$$
然后因为有$i+j$很麻烦，考虑枚举$i+j$就有
∑i=02n∑j=max{0,i−n}min{i,n}ajbi−j(i+j)!j!(i−j)!\sum_{i=0}^{2n}\sum_{j=max\{0,i-n\}}^{min\{i,n\}}a^jb^{i-j}\frac{(i+j)!}{j!(i-j)!}i=0∑2n​j=max{0,i−n}∑min{i,n}​ajbi−jj!(i−j)!(i+j)!​
把$(i+j)!$拿出去就是一个卷积的形式了，模数比较丑所以要用$MTT$来做。（除了$MTT$部分全部自己推？）。

记得要预处理单位根，不然会被卡精度，时间复杂度$O(n\log n)$（常数巨大）

好像有更简单的做法就是用一个$x$满足$ax+bx+c=x$的来消掉$c$这个元就可以直接做了，但是我不会。

---

$code$

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define ll long long
using namespace std;
const double Pi=acos(-1);
const ll N=2e6+10,P=1e6+3;
ll n,a,b,c,ans,F[N],H[N],G[N];
ll inv[N],fac[N],r[N],apw[N],bpw[N];
ll power(ll x,ll b){ll ans=1;while(b){if(b&1)ans=ans*x%P;x=x*x%P;b>>=1;}return ans;
}
ll C(ll n,ll m)
{return fac[n]*inv[m]%P*inv[n-m]%P;}
namespace Poly{const ll seq=32768;struct complex{double x,y;complex(double xx=0,double yy=0){x=xx;y=yy;}}A[N],B[N],C[N],D[N],w[N];complex operator+(complex a,complex b){return complex(a.x+b.x,a.y+b.y);}complex operator-(complex a,complex b){return complex(a.x-b.x,a.y-b.y);}complex operator*(complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}void FFT(complex *f,ll n,ll op){for(ll i=0;i<n;i++)if(i<r[i])swap(f[i],f[r[i]]);for(ll p=2;p<=n;p<<=1){ll len=p>>1;for(ll k=0;k<n;k+=p){for(ll i=k;i<k+len;i++){complex tmp=w[n/len*(i-k)];if(op==-1)tmp.y=-tmp.y;complex tt=f[i+len]*tmp;f[i+len]=f[i]-tt;f[i]=f[i]+tt;}}}if(op==-1){for(ll i=0;i<n;i++)f[i].x=fabs(f[i].x/n+0.5);}return;}void MTT(ll *a,ll *b,ll *c,ll n,ll m){ll l=1;while(l<=n+m)l<<=1;for(ll i=0;i<l;i++)r[i]=(r[i>>1]>>1)|((i&1)?(l>>1):0);for (ll k=1;k<l;k<<=1)for (ll i=0;i<k;i++)w[l/k*i]=(complex){cos(i*Pi/k),sin(i*Pi/k)};for(ll i=0;i<n;i++)A[i].x=a[i]/seq,B[i].x=a[i]%seq;for(ll i=0;i<m;i++)C[i].x=b[i]/seq,D[i].x=b[i]%seq;FFT(A,l,1);FFT(B,l,1);FFT(C,l,1);FFT(D,l,1);complex t1,t2;for(ll i=0;i<l;i++){t1=A[i]*C[i];t2=B[i]*D[i];B[i]=A[i]*D[i]+B[i]*C[i];A[i]=t1;C[i]=t2;}FFT(A,l,-1);FFT(B,l,-1);FFT(C,l,-1);for(ll i=0;i<l;i++){c[i]=(c[i]+(ll)(A[i].x)*seq%P*seq%P)%P;c[i]=(c[i]+(ll)(B[i].x)*seq%P)%P;c[i]=(c[i]+(ll)(C[i].x))%P;}return;}
}
void init(){inv[1]=1;apw[0]=bpw[0]=1;for(ll i=1;i<=2*n;i++)apw[i]=apw[i-1]*a%P,bpw[i]=bpw[i-1]*b%P;for(ll i=2;i<=2*n;i++)inv[i]=(P-(P/i)*inv[P%i]%P)%P;inv[0]=fac[0]=1;for(ll i=1;i<=2*n;i++){fac[i]=fac[i-1]*i%P;inv[i]=inv[i-1]*inv[i]%P;}
}
signed main()
{scanf("%lld%lld%lld%lld",&n,&a,&b,&c);init();n--;for(ll i=0;i<=n;i++){ll x;scanf("%lld",&x);if(!i)continue;x=x*bpw[n-i]%P*apw[n]%P;(ans+=x*C(n+n-i-1,n-1)%P)%=P;}for(ll i=0;i<=n;i++){ll x;scanf("%lld",&x);if(!i)continue;x=x*apw[n-i]%P*bpw[n]%P;(ans+=x*C(n+n-i-1,n-1)%P)%=P;}//处理已知数列for(ll i=0;i<n;i++){F[i]=apw[i]*inv[i]%P;H[i]=bpw[i]*inv[i]%P;}Poly::MTT(F,H,G,n,n);for(ll i=0;i<2*n-1;i++)(ans+=G[i]*c%P*fac[i]%P)%=P;
//    for(ll i=0;i<=n;i++){
//        (ans+=P-apw[n-i]*bpw[n]%P*c%P*C(n+n-i,n)%P)%=P;
//		if(!i)continue;
//        (ans+=P-bpw[n-i]*apw[n]%P*c%P*C(n+n-i,n)%P)%=P;
//    }printf("%lld\n",ans);return 0;
}