正题
题目链接:https://www.luogu.com.cn/problem/P4351
题目大意
\(n*n\)的矩形，给出第一行和第一列的数，剩下的满足\(F_{i,j}=a*F_{i,j-1}+b*F_{i-1,j}+c\)
求\(F_{n,n}\)
解题思路
第一眼看以为是水题，因为给出的数字的贡献通过组合数很好算，但是后来发现麻烦的是那个\(c\)。我们考虑每个格子的\(c\)产生的贡献。
下面为了方便我们先默认让所有格子横纵坐标减\(1\)
对于一个格子\((i,j)\)，通过它的路径有\(\binom{2n-i-j}{n-i}\)种，然后产生的贡献是\(a^{n-i}b^{n-j}\)。为了方便我们反过来表示，然后因为第一行第一列没有贡献所以\(n\)减一。
那么总共的贡献就是
\[\sum_{i=0}^n\sum_{j=0}^na^ib^i\binom{i+j}{i}
\]
然后因为有\(i+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+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;
}