组合数求模C(n,m)%(10^10)

提交地址：http://codeforces.com/gym/100110   B题

http://hi.baidu.com/aekdycoin/item/e051d6616ce60294c5d249d7

几个测试数据

1515151 1213
...0836060000
151144 2002
...3558733440
10000000000 11411
...0000000000
115123131 1210
...2126393000
54546161515130 121231321
...6496000000

由于本人数论实在太水，于是这题搞了我一天多的时间，最后发现是个浮点精度的问题，雷死。。。。。。

题目大意：求C(n,m)%(10^10), n <= 10^18  0<=m<=n

 看到这么大的范围，肯定是有一些数论知识的啦~

注意到10^10可以拆分为2^10  *  5^10

那么是不是可以分别对这两个东东取模然后再合并呢，可以的！！！！

于是现在的任务是求

C(n,m)%2^10  ; 

以及C(n,m)%5^10;

看到了么，其实就是求C(n,m)%(p^a)   p是素数

看了AC大神博客中的那个究极版之后肯定就会做这个了。

举个例子C(n,m)%1024

进而转换为n!%1024

分子分母分别求出来，然后取逆就好了

1*2*3*..*1024

1025*...*2048

...

一直乘到n，很明显，1024个数为一段，结果都是一样的（周期性！！！），最后还剩下几个，暴力求就好了。

注意，我们要想取分母的逆元，分母必须与2互质，所以我们先将2的倍数都剔除掉，显然，周期性仍然满足。

比如现在要求11!%(2^10)

假设calc(n,mod) 可以求n!中所有奇数乘积对mod取余

1*2*3*4*5*6*7*8*9*10*11 

1*3*5*7*9 *11 *    2*4*6*8*10 = 1*3*5*7*9*11 * 2^5*(1*2*3*4*5),   问题可以递归解决！！！！

所以可以写出式子来calc(11,1024) * 2^(11/2) *  calc(5) * (2^11/4)  * calc(2) * 2^(11/8)  * calc(1) 

复杂度是log级别的

这道题比较特殊，只有两个质因子，当P是一般的数的时候分解开来就好了，最后再用中国剩余定理合并

但是如果分解质因数后p^a很大的话就不能用上面的方法了，因为数组开不下，周期没法用。。。。

所以n m k都是10^9级别的数的时候貌似没法搞？？？？

我感觉  K如果是10^9级别的，m 必定要10^6级别才可搞

K如果是10^5级别的，n m 随便都可搞

这道题要判断C(n,m)是否大于10^10,中间用了double，开始的时候一直WA 9，没有检查出来，用eps判断就过了。

还有一种方法是我参考了代码区里的，如下

贴上代码

#include <cstdio>
#include <cstring>
#include <cmath>
typedef long long lld;
const lld mod5 = 9765625;
lld n , m ;
lld fac[10000000];
lld two[64];
lld five[64];
lld Pow(lld a,lld b,lld mod)
{
	lld ans = 1;
	while(b) {
		if(b&1) ans = ans * a % mod;
		a = a*a % mod; b >>= 1;
	}
	return ans;
}
lld exgcd(lld a,lld b,lld &x,lld &y)
{
	if(b==0) {x=1;y=0;return a;}
	else   {
		lld d=exgcd(b,a%b,x,y);
		lld t=x;
		x=y;
		y=t-a/b*y;
		return d;
	}
}
lld inverse(lld num,lld mod)
{
	lld x,y;
	exgcd(num,mod,x,y);
	while(x<0) x+=mod,y-=num;
	return x;
}
void init1()
{
	two[0] = 1;
	for(lld i = 1; i < 64; i++) two[i] = two[i-1]*2;
	fac[1] = 1;fac[0] = 1;
	for(lld i = 2; i <= 1024; i++) {
		lld num = i;
		fac[i] = fac[i-1];
		if(i&1)fac[i] = fac[i]*num % 1024;
	}
}
lld calc(lld n,lld mod)//n!去掉2或5的倍数后对mod取余
{
	lld t = n / mod;
	lld ans = 1 ;
	if(t>=1)ans = Pow(fac[mod],t,mod);
	int lim = n%mod;
	return ans*fac[lim]%mod;
}
lld gao(lld b[],lld mod,lld num) //计算c(n,m) % mod
{
	long long fenzi_ = 0 , fenmu_ = 0;
	lld tn = n , tm = m , tnm = n-m;
	for(lld i=1;b[i]<=tn;fenzi_+=(tn/b[i]),i++);
	for(lld i=1;b[i]<=tm;fenmu_+=(tm/b[i]),i++);
	for(lld i=1;b[i]<=tnm;fenmu_+=(tnm/b[i]),i++);
	lld cnt = fenzi_ - fenmu_;
	lld fenzi = calc(tn,mod);
	for(lld i = 1;b[i]<=tn;i++) fenzi = fenzi*calc(tn/b[i],mod)%mod;
	lld fenmu = calc(tnm,mod);
	for(lld i=1;b[i]<=tnm;i++) fenmu=fenmu*calc(tnm/b[i],mod)%mod;
	fenmu = fenmu * calc(tm,mod) % mod;
	for(lld i=1;b[i]<=tm;i++) fenmu=fenmu*calc(tm/b[i],mod)%mod;
	return fenzi*inverse(fenmu,mod)%mod*Pow(num,cnt,mod)%mod;
}
void init2()
{
	five[0] =1;
	for(lld i = 1; i <= 50; i++) five[i] = five[i-1] * 5;
	fac[1] = 1;  fac[0] = 1;
	for(lld i = 2; i <= mod5; i++) {
		lld num = i;
		fac[i] = fac[i-1];
		if(i%5!=0)  fac[i] = fac[i]*num % mod5;
	}
}
//x % 1024    = a;
//x % mod5    = b;
//上面两个式子用中国剩余定理合并
lld CRT(lld a,lld b)
{
	lld mod = 10000000000LL;
	lld x,y;
	exgcd(mod5,1024,x,y);
	x=x*a%mod;
	lld p = x, q;
	exgcd(1024,mod5,x,y);
	x=x*b%mod;
	q = x;
	lld ans = (mod5*p%mod+1024LL*q%mod)%mod;
	return (ans+mod)%mod;
}
void print(lld ans)
{
	if(n-m<m) m = n-m;
	double t=1;
	lld a=1;
	bool flag = false;
	for(lld i = 1; i <= m; i++)
	{
		t=t*(n-i+1)/i;
		a=a*(n-i+1)/i;
		if(a>=10000000000LL || t>=1e10)
		{
			flag = true;
			break;
		}
	}
	if(flag) printf("...%010I64d\n",ans);
	else printf("%I64d\n",ans);
}
int main()
{
	freopen("combi.in","r",stdin);
	freopen("combi.out","w",stdout);
	while(scanf("%I64d%I64d",&n,&m)!=EOF){
	init1();
	lld a = gao(two,1024,2);
	init2();
	lld b = gao(five,mod5,5);
	lld ans = CRT(a,b);
	print(ans);
	}
	return 0;
}
/*
1515151 1213
...0836060000
151144 2002
...3558733440
10000000000 11411
...0000000000
115123131 1210
...2126393000
54546161515130 121231321
...6496000000
*/