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

HDU 5446

题目链接：

http://acm.hdu.edu.cn/showproblem.php?pid=5446

题意：

求C（n.m)%(p1 * p2 * .. * pn)，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的最小解。

设，(WPS竟然粘不上来)，则

因为本题数据会爆long long，所以需要用快速乘法，具体见代码，类似快速幂。

源码：

#include <cstdio>

#include <cstring>

#include <cmath>

#include <cstdlib>

#include <algorithm>

#include <iostream>

#include <queue>

#include <stack>

using namespace std;

#define LL long long

const 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;

}