HDU 5446——Unknown Treasure——————【CRT+lucas+exgcd+快速乘+递推求逆元】

Each test case starts with three integers n,m,k(1≤m≤n≤1018,1≤k≤10) on a line where k is the number of primes. Following on the next line are kdifferent primes p1,...,pk. It is guaranteed that M=p1⋅p2⋅⋅⋅pk≤1018 and pi≤105 for every i∈{1,...,k}.

 

Output

For each test case output the correct combination on a line.

 

Sample Input

1

9 5 2

3 5

 

Sample Output

6

 

Source

2015 ACM/ICPC Asia Regional Changchun Online

 

题目大意：让求C(n,m)%(∏pi) 这个式子的值。

中国剩余定理：

解题思路：首先用lucas定理将求C(a,b)%p转化成求解∏C(bi,ai)，这样，我们可以得到c[i]数组。然后用中国剩余定理来求x0的值，即为答案。在求解的过程中需要用到扩展欧几里得来求解M_i的逆元，由于Mi比较大，所以在乘积的时候会爆数据范围，所以改成快速乘取模的方式代替直接乘积。

#include<bits/stdc++.h>
using namespace std;
typedef long long INT;
const int maxp=1e5+200;
INT p[15],c[15];
INT fac[maxp],inv[maxp];
INT powmod(INT a,INT n,INT mod){//快速幂取模
    INT ret=1;
    while(n){
        if(n&1){
            ret=ret*a%mod;
        }
        n>>=1;
        a = a*a%mod;
    }
    return ret;
}
INT mulmod(INT a,INT b,INT mod){//快速乘取模
    a = (a%mod + mod) % mod;    //用扩展欧几里得求出的值可能为负值
    b = (b%mod + mod) % mod;    //用扩展欧几里得求出的值可能为负值
    INT ret=0;
    while(b){
        if(b&1){
            ret = (ret+a)%mod;
        }
        b >>= 1;
        a = (a<<1) % mod;
    }
    return ret;
}
void init(INT n){   //递推出来阶乘和逆元数组
    fac[0]=1;
    for(int i=1;i<n;i++){
        fac[i]=fac[i-1]*i % n;
    }
    inv[n-1]=powmod(fac[n-1],n-2,n);
    for(int i=n-2;i>=0;i--){
        inv[i] = inv[i+1] * (i+1) % n;   
        //fac[n]*inv[fac[n]]≡1%p  ==>  fac[n-1]*(n*inv[fac[n]])≡1%p
    }
}
INT cm(INT n,INT m,INT mod){    //用逆元求组合数取模
    if(n<0||m<0||m>n){
        return 0;
    }
    return fac[n]*inv[n-m]%mod*inv[m]%mod;
}
INT lucas(INT n,INT m,INT mod){//lucas递归求P进制时的c
    if(m==0){
        return 1;
    }
    return lucas(n/mod,m/mod,mod) * cm(n%mod,m%mod,mod) % mod;
}
INT exgcd(INT a,INT b,INT &x,INT &y){   //求b关于模a的逆元。放在y中
    if(b==0) { x = 1; y = 0; return a; }
    INT d = exgcd(b, a%b , y, x);
    y -= x * (a / b);
    return d;
}
void CRT(INT k){//中国剩余定理求解一元线性同余方程组
    INT M=1,x,y;
    INT ans=0;
    for(int i=1;i<=k;i++){
        M *= p[i];
    }
    for(int i=1;i<=k;i++){
        INT Mi=M/p[i];
        exgcd(p[i],Mi,x,y);
        ans = (ans+mulmod(mulmod(y,Mi,M),c[i],M))%M                                             ;
    }
    printf("%I64d
",ans);
}
int main(){
    INT n,m,k;
    int t;
    scanf("%d",&t);
    while(t--){
        scanf("%I64d%I64d%I64d",&n,&m,&k);
        for(int i=1;i<=k;i++){
            scanf("%I64d",&p[i]);
            init(p[i]);
            c[i] = lucas(n,m,p[i]);
        }
        CRT(k);
    }
    return 0;
}