HDU 5446 Unknown Treasure(lucas + 中国剩余定理 + 模拟乘法)

题目链接： http://acm.hdu.edu.cn/showproblem.php?pid=5446

题目大意：求C(n, m) % M，其中M为不同素数的乘积,即M=p1*p2*...*pk, 1≤k≤10。1≤m≤n≤10^18。

分析：如果M是素数，则可以直接用lucas定理来做，但是M不是素数，而是素数的连乘积。令C(n, m)为 X ,则可以利用lucas定理分别计算出 X%p1，X%p2, ... , X % pk的值，然后用中国剩余定理来组合得到所求结果。

比较坑的地方是，中间结果会超long long 范围，要用类似于快速幂的方法来模拟。

参考代码：

 1 #include <iostream>
 2 #include <stdio.h>
 3 #include <string.h>
 4 using namespace std;
 5 long long a[110001];
 6 long long Mod;
 7 long long c[110001];
 8 long long w[1000];
 9 long long Egcd(long long a,long long b,long long &x,long long &y){
10     long long d;
11     if(b==0){
12         x=1;y=0;
13         return a;
14     }
15     d=Egcd(b,a%b,y,x);
16     y-=a/b*x;
17     return d;
18 }
19 long long cal(long long a, long long b, long long mod) {//模拟a*b%mod
20     long long ret = 0;
21     while(b) {
22         if(b & 1) ret = (ret + a) % mod;
23         a = (a + a) % mod;
24         b >>= 1;
25     }
26     return ret;
27 }
28 
29 long long c_r(int len){//中国剩余定理
30     int i;
31     long long d,x,y,n,m,ret;
32     ret=0;
33     n=1;
34     for(i=0;i<len;i++)
35         n*=c[i];
36     for(i=0;i<len;i++){
37         m=n/c[i];
38         d=Egcd(m,c[i],y,x);
39         y=(c[i]+y%c[i])%c[i];
40         if( i&1) y -= c[i];
41         //if(y>c[i]-y) y=y-c[i];
42         ret=(ret+cal(y*m, w[i], n))%n;
43     }
44     return (n+ret%n)%n;
45 }
46 
47 void init(){
48     int i;
49     memset(a, 0, sizeof(a));
50     a[0]=1;
51     for(i=1;i<Mod;i++)
52         a[i]=a[i-1]*i%Mod;
53 }
54 
55 long long gcd(long long a,long long b){
56     if(b==0) return a;
57     return gcd(b,a%b);
58 }
59 
60 long long choose(long long n,long long m){
61     if(m>n)return 0;
62     else if(n==m) return 1;
63     long long nn=a[n],mm=a[m]*a[n-m]%Mod;
64     long long d=gcd(nn,mm);
65     nn/=d;
66     mm/=d;
67     long long x,y;
68     Egcd(mm,Mod,x,y);
69     x=(x+Mod)%Mod;
70     return (x*nn)%Mod;
71 }
72 
73 long long work(long long n,long long m){//lucas
74     long long ret=1;
75     while(n&&m){
76         ret*=choose(n%Mod,m%Mod);
77         ret%=Mod;
78         n/=Mod,m/=Mod;
79     }
80     return ret;
81 }
82 
83 int main(){
84     int T;
85     long long n,m;
86     int k;
87     scanf("%d",&T);
88     while(T--){
89         cin >> n >> m >> k;
90         if(m > n-m) m = n-m;
91         for(int i=0;i<k;i++){
92             cin >> c[i];
93             Mod=c[i];
94             init();
95             w[i]=work(n,m);
96         }
97         cout << c_r(k) << endl;
98     }
99 }