組合數求模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...0836060000151144 2002...355873344010000000000 11411...0000000000115123131 1210...212639300054546161515130 121231321...6496000000*/