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

最新推荐文章于 2021-03-20 19:27:54 发布

原创最新推荐文章于 2021-03-20 19:27:54 发布 · 459 阅读

0 ·

CC 4.0 BY-SA版权

----数学题---- 专栏收录该内容

34 篇文章

订阅专栏

本文详细介绍了使用Lucas定理和快速乘法解决组合问题的方法，包括如何利用Lucas定理求解大数值组合数取模问题，以及如何通过快速乘法提高计算效率。同时，阐述了中国剩余定理在求解同余方程中的应用，为解决复杂组合数学问题提供了一套有效策略。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

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