3431. 【GDOI2014模拟】网格

Description

某城市的街道呈网格状，左下角坐标为A(0, 0)，右上角坐标为B(n, m)，其中n >= m。现在从A(0, 0)点出发，只能沿着街道向正右方或者正上方行走，且不能经过图示中直线左上方的点，即任何途径的点(x, y)都要满足x >= y，请问在这些前提下，到达B(n, m)有多少种走法。

Input

输入文件中仅有一行，包含两个整数n和m，表示城市街区的规模。

Output

输出文件中仅有一个整数和一个换行/回车符，表示不同的方案总数。

Sample Input

输入1：

6 6

输入2：

5 3

Sample Output

输出1：

132

输出2：

28

Data Constraint

50%的数据中，n = m，在另外的50%数据中，有30%的数据：1 <= m < n <= 100

100%的数据中，1 <= m <= n <= 5 000

Solution

对于n=m的情况：，即为卡特兰数的第n项。

因为必须有n步向上，m步向右。

向右的个数大于等于向上的个数，把向右看成左括号，向上看成右括号，即为括号序问题，答案为卡特兰数的第n项。

通项公式为f[ n ]=C(2n,n)-C(2n,n+1)。

对于n=m或者n!=m的情况：

考虑Ans=从(0,0)走到(n,m)的总方案数-不合法方案数。

从（0,0）走到（n,m）一共要向右走n次，向上走m次，一共n+m次，任意选出n次向右，其他向上，总方案数为C(n+m,n)。

考虑不合法方案数，即 穿过 直线y=x的路线的方案数。

考虑 穿过 直线y=x，就相当于碰到（即某一点走到或穿过直线上）直线y=x+1的路径的方案数。

那么将从起点（0,0）开始的路径中，最后一个与直线y=x+1接触的点之前的路径全部关于y=x+1对称过去,就得到了一个从（-1,1）开始的走到（n,m）的只能向上或向右走的路径，每一条路径都可以还原成一个从起点（0,0）开始的不合法的路径。

其方案数为C(n-(-1)+m-1,n-(-1))=C(n+m,n+1)。

答案为C(n+m,n)-C(n+m,n+1)。

考虑转化为阶乘，对乘积进行质因数分解，除法相当于直接将指数相减，最后做一次高精度乘法和高精度减法即可。

Code

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define I int
#define ll long long
#define F(i,a,b) for(ll i=a;i<=b;i++)
#define Fd(i,a,b) for(I i=a;i>=b;i--)
#define mem(a,b) memset(a,b,sizeof a)
#define M 1000000000
using namespace std;
I n,m,p[10010],bz[10010],f[10010],now;
ll a[10010],b[10010],c[10010],r;
I main(){
	freopen("grid.in","r",stdin);
	freopen("grid.out","w",stdout);
	scanf("%d%d",&n,&m);
	F(i,2,10000){
		if(!bz[i]) p[++p[0]]=i;
		F(j,1,p[0]){
			if(i*p[j]>10000) break;
			bz[i*p[j]]=1;
			if(i%p[j]==0) break;
		}
	}
	mem(f,0);
	F(i,1,p[0]){
		now=n+m;
		if(p[i]>now) break;
		while(now){
			f[i]+=now/p[i];
			now/=p[i];
		}
	}
	F(i,1,p[0]){
		now=n;
		if(p[i]>now) break;
		while(now){
			f[i]-=now/p[i];
			now/=p[i];
		}
	}
	F(i,1,p[0]){
		now=m;
		if(p[i]>now) break;
		while(now){
			f[i]-=now/p[i];
			now/=p[i];
		}
	}
	a[0]=a[1]=1;
	F(i,1,p[0]){
		while(f[i]){
			F(j,1,a[0]) a[j]*=p[i];
			F(j,1,a[0]){
				a[j+1]+=a[j]/M;
				a[j]%=M;
			}
			while(a[a[0]+1]){
				a[0]++;
				a[a[0]+1]+=a[a[0]]/M;
				a[a[0]]%=M;
			}
			f[i]--;
		}
	}
	mem(f,0);
	F(i,1,p[0]){
		now=n+m;
		if(p[i]>now) break;
		while(now){
			f[i]+=now/p[i];
			now/=p[i];
		}
	}
	F(i,1,p[0]){
		now=n+1;
		if(p[i]>now) break;
		while(now){
			f[i]-=now/p[i];
			now/=p[i];
		}
	}
	F(i,1,p[0]){
		now=m-1;
		if(p[i]>now) break;
		while(now){
			f[i]-=now/p[i];
			now/=p[i];
		}
	}
	b[0]=b[1]=1;
	F(i,1,p[0]){
		while(f[i]){
			F(j,1,b[0]) b[j]*=p[i];
			F(j,1,b[0]){
				b[j+1]+=b[j]/M;
				b[j]%=M;
			}
			while(b[b[0]+1]){
				b[0]++;
				b[b[0]+1]+=b[b[0]]/M;
				b[b[0]]%=M;
			}
			f[i]--;
		}
	}
	F(i,1,a[0]) c[i]=a[i]-b[i];
	F(i,1,a[0]) if(c[i]<0){
		c[i]+=M;
		c[i+1]--;
	}
	c[0]=a[0];
	while(!c[c[0]]&&c[0]) c[0]--;
	printf("%d",c[c[0]]);
	Fd(i,c[0]-1,1) printf("%09d",c[i]);
	return 0;
}