HDU 5446 Lucas + 中国剩余定理 + 快速乘法

HDU 5446
题目链接：
http://acm.hdu.edu.cn/showproblem.php?pid=5446
题意：
求，n，m的范围很大
思路：
Lucas + 中国剩余定理 + 快速乘法。
Lucas定理用来求当n和m很大时的C（n,m）%p（p是素数）。证明网上有，具体公式就是
Lucas(n,m,p) = C(n%p,m%p）*Lucas(n/p,m/p,p)
中国剩余定理用来求同余方程。对方程n的最小解。
设，，则
因为本题数据会爆long long，所以需要用快速乘法，具体见代码，类似快速幂。
源码：

#include <cstdio>#include <cstring>#include <cmath>#include <cstdlib>#include <algorithm>#include <iostream>#include <queue>#include <stack>using namespace std;#define LL long longconst int MAXN = 15;LL prime[MAXN];LL fac[100005], inv[100005];LL lv[MAXN];void exceed_gcd(LL a, LL b, LL &ans, LL &x, LL &y){    if(a < b){        exceed_gcd(b, a, ans, y, x);        return;    }    if(b == 0){        x = 1, y = 0;        ans = a;        return;    }    else{        exceed_gcd(b, a % b, ans, y, x);        y -= x * (a / b);    }}LL ppow(LL a, LL x, LL mod){    LL ans = 1;    while(x){//        printf("mod = %I64d, a = %I64d, x = %I64d, ans = %I64d\n", mod, a, x, ans);        if(x & 1)            ans = (ans * a) % mod;        a = (a * a) % mod;        x >>= 1;    }    return ans;}LL C(LL n, LL m, LL mod){    if(n < m || n < 0 || m < 0) return 0;//    printf("n = %I64d, m = %I64d, fac[n] = %I64d, in[m] = %I64d, inv[n - m] = %I64d\n", n, m, fac[n], inv[m], inv[n - m]);    return fac[n] * inv[m] % mod * inv[n - m] % mod;}LL lucas(LL n, LL m, LL mod){//    printf("n = %I64d, m = %I64d\n", n, m);    if(n < m)   {/*printf("n=%I64d,m=%I64d\n",n,m);*/return 0;}    if(n == 0 || m == 0)    return 1;//    printf("n = %I64d, m = %I64d, C = %I64d, lucas = %I64d\n", n, m, C(n % mod, m % mod, mod), lucas(n / mod, m / mod, mod));    return C(n % mod, m % mod, mod) * lucas(n / mod, m / mod, mod) % mod;}LL mul(LL a, LL b, LL mod){    a = (a % mod + mod) % mod;    b = (b % mod + mod) % mod;    LL ans = 0;    while(b){        if(b & 1)            ans = (ans + a) % mod;        b >>= 1;        a <<= 1;        a = a % mod;    }    return ans;}int main(){//    freopen("1006.in", "r", stdin);    LL n, m, cnt;    int t;    scanf("%d", &t);    while(t--){        scanf("%I64d%I64d%I64d", &n, &m, &cnt);        LL M = 1;        LL d, y;        LL ans = 0;        for(int i = 0 ; i < cnt ; i++){            scanf("%I64d", &prime[i]);            M *= prime[i];            fac[0] = 1;            for(int j = 1 ; j < prime[i] ; j++)                fac[j] = fac[j - 1] * j % prime[i];            inv[prime[i] - 1] = ppow(fac[prime[i] - 1], prime[i] - 2, prime[i]);            for(int j = prime[i] - 2 ; j >= 0 ; j--)                inv[j] = (inv[j + 1] * (j + 1)) % prime[i];//            for(int j = 0 ; j < prime[i] ; j++)//                printf("fac[%d] = %I64d, inv[%d] = %I64d\n", j, fac[j], j, inv[j]);            lv[i] =  lucas(n, m, prime[i]);        }        for(int i = 0 ; i < cnt ; i++){            LL M1 = M / prime[i];            LL M2;            exceed_gcd(prime[i], M1, d, y, M2);//            printf("i = %d, M1 = %I64d, M2 = %I64d, M = %I64d lv[%d] = %I64d\n", i, M1, M2, M, i, lv[i]);//            M2 = (M2 % M + M) % M;            ans = (ans + mul(mul(lv[i], M1, M), M2, M));        }        ans = (ans % M + M) % M;        printf("%I64d\n", ans);    }    return 0;}