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 范围，要用类似于快速幂的方法来模拟。

参考代码：

#include <iostream>
#include <stdio.h>
#include <string.h>
using namespace std;
long long a[110001];
long long Mod;
long long c[110001];
long long w[1000];
long long Egcd(long long a,long long b,long long &x,long long &y){
    long long d;
    if(b==0){
        x=1;y=0;
        return a;
    }
    d=Egcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
long long cal(long long a, long long b, long long mod) {//模拟a*b%mod
    long long ret = 0;
    while(b) {
        if(b & 1) ret = (ret + a) % mod;
        a = (a + a) % mod;
        b >>= 1;
    }
    return ret;
}

long long c_r(int len){//中国剩余定理
    int i;
    long long d,x,y,n,m,ret;
    ret=0;
    n=1;
    for(i=0;i<len;i++)
        n*=c[i];
    for(i=0;i<len;i++){
        m=n/c[i];
        d=Egcd(m,c[i],y,x);
        y=(c[i]+y%c[i])%c[i];
        if( i&1) y -= c[i];
        //if(y>c[i]-y) y=y-c[i];
        ret=(ret+cal(y*m, w[i], n))%n;
    }
    return (n+ret%n)%n;
}

void init(){
    int i;
    memset(a, 0, sizeof(a));
    a[0]=1;
    for(i=1;i<Mod;i++)
        a[i]=a[i-1]*i%Mod;
}

long long gcd(long long a,long long b){
    if(b==0) return a;
    return gcd(b,a%b);
}

long long choose(long long n,long long m){
    if(m>n)return 0;
    else if(n==m) return 1;
    long long nn=a[n],mm=a[m]*a[n-m]%Mod;
    long long d=gcd(nn,mm);
    nn/=d;
    mm/=d;
    long long x,y;
    Egcd(mm,Mod,x,y);
    x=(x+Mod)%Mod;
    return (x*nn)%Mod;
}

long long work(long long n,long long m){//lucas
    long long ret=1;
    while(n&&m){
        ret*=choose(n%Mod,m%Mod);
        ret%=Mod;
        n/=Mod,m/=Mod;
    }
    return ret;
}

int main(){
    int T;
    long long n,m;
    int k;
    scanf("%d",&T);
    while(T--){
        cin >> n >> m >> k;
        if(m > n-m) m = n-m;
        for(int i=0;i<k;i++){
            cin >> c[i];
            Mod=c[i];
            init();
            w[i]=work(n,m);
        }
        cout << c_r(k) << endl;
    }
}