标签：ll return ifac int 191NF 贡献 Frightful Formula mo

一个矩阵，给出第一行第一列，后面每个位置的数由此得到：\(F_{i,j}=F_{i,j-1}*a+F_{i-1,j}*b+c\)

求\(F_{n,n}\)。

答案对\(10^6+3\)取模。

\(n\le 2*10^5\)

如果\(c=0\)，考虑组合意义：对于从\((x,1)\)，贡献相当于\((x,1)\)到\((n,n)\)，每次可以向右或向下，向右有\(a\)的贡献，向下有\(b\)的贡献。第一步必须要向右。简单组合数即可计算。

如果\(c>0\)，类似地考虑除第一行第一列外的点，相当于它一开始有个\(c\)的贡献，向右或向下走到\((n,n)\)。这一部分显然可以列出式子：\(\sum_{x=0}^{n-2}\sum_{y=0}^{n-2}\binom{x+y}{x}a^xb^y\)。

这里其实可以直接卷积了。由于答案的模数比较小，所以可以直接FFT或大质数NTT（模\(29*2^{53}+1\)，原根为\(3\)）。

其实也不用。按照\(x+y\)的值分类计算，\(x+y\le n-2\)的贡献直接\(\sum_{s=0}^{n-2}(a+b)^s\)计算。

剩下的部分，考虑放到原来的题目中：先算出\(x+y=n-2\)的贡献，然后算下一个。乘\((a+b)\)相当于向右或向下走一格，那这个时候把走出矩阵的不合法贡献减掉即可。

NTT做法

using namespace std;
#include <cstdio>
#include <cstring>
#include <algorithm>
#define N 524288
#define ll long long
const int mo=1000003;
ll qpow(ll x,ll y=mo-2){
	ll r=1;
	for (;y;y>>=1,x=x*x%mo)
		if (y&1)
			r=r*x%mo;
	return r;
}
ll fac[N*2],ifac[N*2];
void initC(int n){
	fac[0]=1;
	for (int i=1;i<=n;++i)
		fac[i]=fac[i-1]*i%mo;
	ifac[n]=qpow(fac[n]);
	for (int i=n-1;i>=0;--i)
		ifac[i]=ifac[i+1]*(i+1)%mo;
}
ll C(int m,int n){return fac[m]*ifac[n]%mo*ifac[m-n]%mo;}
int n,a,b,c;
ll A[N],B[N],s[N];
int p[N],q[N];
namespace NTT{
	const ll mo=(29ll<<57)+1;
	ll multi(ll x,ll y){
//		return x*y%mo;
		//x%=mo,y%=mo
		x=(x+mo)%mo;
		y=(y+mo)%mo;
		ll d=(long double)x*y/mo;
		d=x*y-d*mo;
		if (d<0) d+=mo;
		if (d>=mo) d-=mo;
//		printf("%lld*%lld%%%lld=%lld\n",x,y,mo,d);
//		printf("%lf\n",((double)x*y-d)/mo);
		return d;
	}
	int re[N],nN;
	ll qpow(ll x,ll y=mo-2){
		ll r=1;
		for (;y;y>>=1,x=multi(x,x))
			if (y&1)
				r=multi(r,x);
		return r;
	}
	void setlen(int n){
		int bit=0;
		for (nN=1;nN<=n;++bit,nN<<=1);
		re[0]=0;
		for (int i=1;i<nN;++i)
			re[i]=re[i>>1]>>1|(i&1)<<bit-1;
	}
	void dft(ll A[],int flag){
		for (int i=0;i<nN;++i)
			if (i<re[i])
				swap(A[i],A[re[i]]);
		static ll wnk[N];
		for (int i=1;i<nN;i<<=1){
			ll wn=qpow(3,flag==1?(mo-1)/(2*i):mo-1-(mo-1)/(2*i));
			wnk[0]=1;
			for (int k=1;k<i;++k)
				wnk[k]=multi(wnk[k-1],wn);
			for (int j=0;j<nN;j+=i<<1)
				for (int k=0;k<i;++k){
					ll x=A[j+k],y=multi(A[j+k+i],wnk[k]);
					A[j+k]=(x+y)%mo;
					A[j+k+i]=(x-y)%mo;
				}
		}
		for (int i=0;i<nN;++i)
			A[i]=(A[i]+mo)%mo;
		if (flag==-1){
			ll invn=qpow(nN);
			for (int i=0;i<nN;++i)
				A[i]=multi(A[i],invn);
		}
	}
	void multi(ll c[],ll a[],ll b[]){
//		for (int i=0;i<=n;++i)
//			for (int j=0;j<=n;++j)
//				(c[i+j]+=multi(a[i],b[j]))%=mo;
		setlen(n*2);
		static ll A[N],B[N];
		for (int i=0;i<=n;++i)
			A[i]=a[i],B[i]=b[i];
		dft(A,1),dft(B,1);
		for (int i=0;i<nN;++i)
			A[i]=multi(A[i],B[i]);
		dft(A,-1);
		for (int i=0;i<=n*2;++i)
			c[i]=A[i];
	}
}
int main(){
//	freopen("in.txt","r",stdin);
//	freopen("out.txt","w",stdout);
	scanf("%d%d%d%d",&n,&a,&b,&c);
	A[0]=B[0]=1;
	for (int i=1;i<=n;++i){
		A[i]=A[i-1]*a%mo;
		B[i]=B[i-1]*b%mo;
	}
	for (int i=1;i<=n;++i) scanf("%d",&p[i]);
	for (int i=1;i<=n;++i) scanf("%d",&q[i]);
	initC(n*2);
	ll ans=0,sum=0;
	for (int i=2;i<=n;++i) (ans+=A[n-1]*B[n-i]%mo*p[i]%mo*C(n-2+n-i,n-2))%=mo;
	for (int i=2;i<=n;++i) (ans+=B[n-1]*A[n-i]%mo*q[i]%mo*C(n-2+n-i,n-2))%=mo;
	n-=2;
	for (int i=0;i<=n;++i){
		A[i]=A[i]*ifac[i]%mo;
		B[i]=B[i]*ifac[i]%mo;
	}
//	for (int i=0;i<=n;++i) printf("%lld ",A[i]);
//	printf("\n");
//	for (int i=0;i<=n;++i) printf("%lld ",B[i]);
//	printf("\n");
	NTT::multi(s,A,B);
//	for (int i=0;i<=n*2;++i) printf("%lld ",s[i]);
//	printf("\n");
	for (int i=0;i<=n*2;++i)
		(sum+=s[i]%mo*fac[i])%=mo;
	(ans+=sum*c)%=mo;
	printf("%lld\n",(long long)ans);
	return 0;
}

标签：ll,return,ifac,int,191NF,贡献,Frightful,Formula,mo
来源： https://www.cnblogs.com/jz-597/p/13916005.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。