正題

題目鏈接: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)$，通過它的路徑有$\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$就有
$$\sum _{i=0}^{2n}\sum _{j=max\{0,i?n\}}^{min\{i,n\}}{a}^j{b}^{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; }

總結

以上是生活随笔為你收集整理的P4351-[CERC2015]Frightful Formula【组合数学,MTT】的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。