把博客园图标替换成自己的图标
把博客园图标替换成自己的图标end

【洛谷4351】[CERC2015] Frightful Formula(坐标系走路)

点此看题面

  • 给定一个\(n\times n\)矩阵的第一行和第一列,其余格子满足\(f_{i,j}=a\times f_{i,j-1}+b\times f_{i-1,j}+c\),求\(f_{n,n}\)
  • \(n\le2\times10^5\)

坐标系走路+任意模数\(NTT\)

\(f_{i,j}\)看作是走到\((i,j)\)的方案数。

那么考虑这个式子的意义,可以认为是任何时候往右走有\(a\)种方案,往下走有\(b\)种方案,且有\(c\)种方式选择从某个点\((i,j)\)出发。

先特殊处理掉第一行第一列:\((1,1)\)是没用的,第一行的格子第一步必须往下,第一列的格子第一步必须往右。具体处理方式可参考一般情况,不过需要略加修改。

对于其他格子\((i,j)\),有\(c\)种方式从这里出发,然后共需要往右走\(n-j\)步,往下走\(n-i\)步,走的顺序可以任意安排,因此对答案的贡献是:

\[c\times a^{n-j}\times b^{n-i}\times C_{n-j+n-i}^{n-j} \]

拆开组合数:

\[c\times \frac{a^{n-j}}{(n-j)!}\times \frac{b^{n-i}}{(n-i)!}\times (n-j+n-i)! \]

发现这可以看成一个卷积的形式。设\(F_i=\sum_{j+k=i}\frac{a^j}{j!}\times\frac{b^k}{k!}\),那么答案就是:

\[c\times\sum_{i=0}^{2(n-2)}f_i\times i! \]

由于此题模数是\(10^6+3\),需要任意模数\(NTT\)

代码:\(O(nlogn)\)

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Rg register
#define RI Rg int
#define Cn const
#define CI Cn int&
#define I inline
#define W while
#define N 200000
#define X 1000003
#define DB long double
#define C(x,y) (1LL*Fac[x]*IFac[y]%X*IFac[(x)-(y)]%X)
using namespace std;
int n,a,b,c,p[N+5],q[N+5],A[N+5],B[N+5],f[2*N+5],Fac[2*N+5],IFac[2*N+5];
I int QP(RI x,RI y,CI p) {RI t=1;W(y) y&1&&(t=1LL*t*x%p),x=1LL*x*x%p,y>>=1;return t;}
namespace Poly//多项式模板
{
	#define PR 3
	#define X1 998244353
	#define X2 1004535809
	int P,L,R[N<<2],A[N<<2],B[N<<2],S1[N<<2],S2[N<<2];
	I void NTT(CI p,int* s,CI op)
	{
		RI i,j,k,x,y,U,S;for(i=0;i^P;++i) i<R[i]&&(swap(s[i],s[R[i]]),0);
		for(i=1;i^P;i<<=1) for(U=QP(QP(PR,op,p),(p-1)/(i<<1),p),j=0;j^P;j+=i<<1) for(S=1,
			k=0;k^i;++k,S=1LL*S*U%p) s[j+k]=((x=s[j+k])+(y=1LL*S*s[i+j+k]%p))%p,s[i+j+k]=(x-y+p)%p;
		if(op==p-2) for(x=QP(P,p-2,p),i=0;i^P;++i) s[i]=1LL*s[i]*x%p;
	}
	I int Calc(CI x,CI y) {RI k=1LL*(y-x+X2)*QP(X1,X2-2,X2)%X2;return (1LL*k*X1+x)%X;}
	I void Mul(int* f,int* a,int* b)//任意模数多项式乘法
	{
		RI i;P=1,L=0;W(P<=2*n) P<<=1,++L;for(i=0;i^P;++i) R[i]=(R[i>>1]>>1)|((i&1)<<L-1);
		for(i=0;i^P;++i) A[i]=B[i]=0;for(i=0;i<=n;++i) A[i]=a[i],B[i]=b[i];NTT(X1,A,1),NTT(X1,B,1);
		for(i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X1;for(NTT(X1,A,X1-2),i=0;i<=2*n;++i) S1[i]=A[i];
		for(i=0;i^P;++i) A[i]=B[i]=0;for(i=0;i<=n;++i) A[i]=a[i],B[i]=b[i];NTT(X2,A,1),NTT(X2,B,1);
		for(i=0;i^P;++i) A[i]=1LL*A[i]*B[i]%X2;for(NTT(X2,A,X2-2),i=0;i<=2*n;++i) S2[i]=A[i];
		for(i=0;i<=2*n;++i) f[i]=Calc(S1[i],S2[i]);
	}
}
int main()
{
	RI i,j,t=0;for(scanf("%d%d%d%d",&n,&a,&b,&c),i=1;i<=n;++i) scanf("%d",p+i);for(i=1;i<=n;++i) scanf("%d",q+i);
	for(Fac[0]=i=1;i<=2*n;++i) Fac[i]=1LL*Fac[i-1]*i%X;for(IFac[i=2*n]=QP(Fac[2*n],X-2,X);i;--i) IFac[i-1]=1LL*IFac[i]*i%X;
	for(i=2;i<=n;++i) t=(t+1LL*p[i]*QP(a,n-1,X)%X*QP(b,n-i,X)%X*C(n-i+n-2,n-i))%X;//特殊处理第一列
	for(j=2;j<=n;++j) t=(t+1LL*q[j]*QP(a,n-j,X)%X*QP(b,n-1,X)%X*C(n-2+n-j,n-j))%X;//特殊处理第一行
	for(i=0;i<=n-2;++i) A[i]=1LL*QP(a,i,X)*IFac[i]%X,B[i]=1LL*QP(b,i,X)*IFac[i]%X;//写出需要卷积的两个数组
	RI o=0;for(Poly::Mul(f,A,B),i=0;i<=2*(n-2);++i) o=(o+1LL*f[i]*Fac[i])%X;return printf("%d\n",(t+1LL*c*o)%X),0;//卷积后统计答案
}
posted @ 2021-07-22 19:18  TheLostWeak  阅读(69)  评论(0编辑  收藏  举报