组合数求模C(n,m)%(10^10)
几个测试数据
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判断就过了。
还有一种方法是我参考了代码区里的，如下
贴上代码
[cpp]
#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 ex易做图(lld a,lld b,lld &x,lld &y)
{
if(b==0) {x=1;y=0;return a;}
else {
lld d=ex易做图(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;
ex易做图(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;
ex易做图(mod5,1024,x,y);
x=x*a%mod;
lld p = x, q;
ex易做图(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);
} <
补充：软件开发 , C++ ,
