組合數求模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 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
*/ 

#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
*/