P4351-[CERC2015]Frightful Formula【組合數學,MTT】
阿新 • • 發佈:2021-01-10
正題
題目連結: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}\)
那麼總共的貢獻就是 \[\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; }