【模板】lucas定理和扩展lucas定理(组合数取模)

本文介绍了Lucas定理及其扩展版本,详细讲解了当模数为质数时的Lucas定理,并提供了对应的代码实现;同时探讨了模数不是质数时的扩展Lucas定理,给出了具体的算法思路与实现代码。

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

参考博客:

https://ksmeow.moe/lucas_theorem/

https://blog.csdn.net/clove_unique/article/details/54571216/

1、lucas定理(模数是质数)

ll mul(ll a, ll b) {return (a%mod)*(b%mod)%mod;}
ll qpow(ll a, ll b)
{
    ll res = 1; a%=mod;
    while(b)
    {
        if(b&1) res = mul(res, a);
        a = mul(a, a);
        b>>=1;
    }
    return res;
}
ll com(ll a, ll b)//Ca^b
{
    if(a<b) return 0;
    if(b==0 || a==b) return 1;
    b = min(b, a-b);
    ll ca = 1, cb = 1;
    for(int i = 0; i < b; i++)//此处可预处理阶乘和逆元
    {
        ca = mul(ca, a-i);
        cb = mul(cb, i+1);
    }
    return  mul(ca, qpow(cb, mod-2));
}
ll lucas(ll a, ll b)
{
    if(b==0) return 1;
    return com(a%mod, b%mod)*lucas(a/mod, b/mod)%mod;
}

2、扩展lucas定理(模数不是质数)

codeforces2015ICL,Finals,Div.1#J

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
#define LL long long

LL n,m,MOD,ans;

LL fast_pow(LL a,LL p,LL Mod)
{
    LL ans=1LL;
    for (;p;p>>=1,a=a*a%Mod)
        if (p&1)
            ans=ans*a%Mod;
    return ans;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{
    if (!b) x=1LL,y=0LL;
    else exgcd(b,a%b,y,x),y-=a/b*x;
}
LL inv(LL A,LL Mod)
{
    if (!A) return 0LL;
    LL a=A,b=Mod,x=0LL,y=0LL;
    exgcd(a,b,x,y);
    x=((x%b)+b)%b;
    if (!x) x+=b;
    return x;
}
LL Mul(LL n,LL pi,LL pk)
{
    if (!n) return 1LL;
    LL ans=1LL;
    if (n/pk)
    {
        for (LL i=2;i<=pk;++i)
            if (i%pi) ans=ans*i%pk;
        ans=fast_pow(ans,n/pk,pk);
    }
    for (LL i=2;i<=n%pk;++i)
        if (i%pi) ans=ans*i%pk;
    return ans*Mul(n/pi,pi,pk)%pk;
}
LL C(LL n,LL m,LL Mod,LL pi,LL pk)
{
    if (m>n) return 0LL;
    LL a=Mul(n,pi,pk),b=Mul(m,pi,pk),c=Mul(n-m,pi,pk);
    LL k=0LL,ans;
    for (LL i=n;i;i/=pi) k+=i/pi;
    for (LL i=m;i;i/=pi) k-=i/pi;
    for (LL i=n-m;i;i/=pi) k-=i/pi;
    ans=a*inv(b,pk)%pk*inv(c,pk)%pk*fast_pow(pi,k,pk)%pk;
    return ans*(Mod/pk)%Mod*inv(Mod/pk,pk)%Mod;
}
int main()
{
    scanf("%I64d%I64d%I64d",&n,&m,&MOD);
    for (LL x=MOD,i=2;i<=MOD;++i)
        if (x%i==0)
        {
            LL pk=1LL;
            while (x%i==0) pk*=i,x/=i;
            ans=(ans+C(n,m,MOD,i,pk))%MOD;
        }
    printf("%I64d\n",ans);
}

洛谷P2183

#include <cstdio>

typedef long long LL;

struct Num {
    LL num, pcnt; // pcnt: p因子数量,num: 提取p因子后乘起来的积
    Num(LL num = 0, LL pcnt = 0): num(num), pcnt(pcnt) {}
};

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF 
        : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    char c = fgc();
    while(c < '0' || c > '9') {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 100005, MAXM = 10;
const LL INF = 1e15;

inline LL fpow(LL n, LL k, LL p) {
    LL res = 1; n %= p;
    while(k) {
        if(k & 1) res = res * n % p;
        n = n * n % p;
        k >>= 1;
    }
    return res;
}

LL P, n, m, w[MAXM];

// calculate prime divisors in x
LL divi[MAXN], dcnt[MAXN], dpow[MAXN], dtot;

inline void calDivisor(LL x) {
    dtot = 0;
    for(LL i = 2; i * i <= P; i++) {
        if(x % i == 0) {
            divi[++dtot] = i;
            while(x % i == 0) {
                x /= i; dcnt[dtot]++;
            }
            dpow[dtot] = fpow(i, dcnt[dtot], INF);
        }
    }
    if(x != 1) {
        dtot++; 
        divi[dtot] = dpow[dtot] = x; dcnt[dtot] = 1;
    }
}

// calculate x! mod p^k
inline Num calFactorial(LL x, LL id) {
    if(x == 1 || x == 0) return Num(1, 0); // 1! = 0! = 1
    LL pk = dpow[id], res = 1;
    LL pnum = 0;
    for(LL i = x; i; i /= divi[id]) pnum += i / divi[id]; // 计算p因子的数量
    Num nxt = calFactorial(x / divi[id], id); // 递归计算提取了一个p因子的p倍数组成的阶乘
    res = res * nxt.num % pk; // 合并答案
    if(x >= pk) { // 如果有多段超过p^k的,预处理应用快速幂
        LL fac = 1;
        for(LL i = 1; i < pk; i++) {
            if(i % divi[id] == 0) continue;
            fac = fac * i % pk;
        }
        res = res * fpow(fac, x / pk, pk) % pk;
    }
    for(LL i = x; i >= 1; i--) { // 非整段p^k,暴力求解
        if(i % pk == 0) break;
        if(i % divi[id] == 0) continue;
        res = res * i % pk;
    }
    return Num(res, pnum);
}

inline LL exgcd(LL a, LL b, LL &x, LL &y) {
    if(!b) {
        x = 1; y = 0; return a;
    }
    LL res = exgcd(b, a % b, x, y);
    LL t = x; x = y; y = t - a / b * y;
    return res;
}

inline LL calInverse(LL n, LL p) {
    LL x, y;
    exgcd(n, p, x, y);
    return (x % p + p) % p;
}

LL crtm[MAXN];

// CRT求解组合数
inline LL crt() {
    LL res = 0;
    for(int i = 1; i <= dtot; i++) {
        res = (res + crtm[i] * P / dpow[i] % P * calInverse(P / dpow[i], dpow[i]) % P) % P;
    }
    return res;
}

int main() {
    P = readint(); n = readint(); m = readint();
    LL wsum = 0;
    for(int i = 1; i <= m; i++) {
        w[i] = readint();
        wsum += w[i];
    }
    if(wsum > n) {
        puts("Impossible"); return 0;
    }
    calDivisor(P);
    LL ans = 1;
    for(int i = 1; i <= m; i++) {
        for(int j = 1; j <= dtot; j++) {
            Num N = calFactorial(n, j), 
                M = calFactorial(w[i], j), 
                NM = calFactorial(n - w[i], j);
            // 分别代表n!、m!、(n-m)!
            N.pcnt -= M.pcnt + NM.pcnt;
            if(N.pcnt >= dcnt[j]) { // 如果p因子更多,说明能整除,模意义下为0
                crtm[j] = 0;
            } else {
                crtm[j] = N.num * calInverse(M.num, dpow[j]) % dpow[j] 
                    * calInverse(NM.num, dpow[j]) % dpow[j] 
                    * fpow(divi[j], N.pcnt, dpow[j]) % dpow[j]; // 把p因子乘进去
            }
        }
        ans = ans * crt() % P;
        n -= w[i];
    }
    printf("%lld", ans);
    return 0;
}

 

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值