专栏首页CSDN旧文『数学』--数论--组合数+卢卡斯定理+扩展卢卡斯定理

『数学』--数论--组合数+卢卡斯定理+扩展卢卡斯定理

组合数: 在N个数中选取M个数,问选的方式有几种?

直接递归暴力简单

#include<cstdio>
const int N = 2000 + 5;
const int MOD = (int)1e9 + 7;
int comb[N][N];//comb[n][m]就是C(n,m)
void init(){
    for(int i = 0; i < N; i ++){
        comb[i][0] = comb[i][i] = 1;
        for(int j = 1; j < i; j ++){
            comb[i][j] = comb[i-1][j] + comb[i-1][j-1];
            comb[i][j] %= MOD;
        }
    }
}
int main(){
    init();
}

卢卡斯定理: 求 C m n   m o d   p C_m^n~mod~p Cmn​ mod p 设 m = a 0 p 0 + a 1 p 1 + ⋯ + a k p k m={a_0}^{p_0}+{a_1}^{p_1}+\cdots+{a_k}^{p_k} m=a0​p0​+a1​p1​+⋯+ak​pk​ 设 n = b 0 p 0 + b 1 p 1 + ⋯ + b k p k n={b_0}^{p_0}+{b_1}^{p_1}+\cdots+{b_k}^{p_k} n=b0​p0​+b1​p1​+⋯+bk​pk​ 则 C m n ≡ ∏ C a i b i ( m o d   p ) C_m^n\equiv\prod{C_{a_i}^{b_i}}(mod~p) Cmn​≡∏Cai​bi​​(mod p)

证明我就不打了,百度百科上有,数学符号太多了~ _ ~ 前置知识: 1.欧几里得或者费马小定理求逆元 2.快速幂 实现代码:

#include<iostream>
using namespace std;
typedef long long LL;
const LL N=1e5+2;
LL a[N];
void init(LL p)
{
	a[1]=1;
	for(int i=2;i<=p;++i)a[i]=a[i-1]*i%p;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{
	if(!b){
		x=1;
		y=0;
		return;
	}
	exgcd(b,a%b,y,x);
	y-=a/b*x;
}
LL ksm(LL x,LL n,LL mod)
{
	LL ans=1;
	while(n){
		if(n&1)ans=ans*x%mod;
		n>>=1;
		x=x*x%mod;
	}
	return ans;
}
LL C(LL n,LL m,LL p)
{
	if(n==m||m==0)return 1;
	if(n<m)return 0;
	if(m*2>n)m=n-m;						  /*C(n,m)=c(n,n-m)*/
	return a[n]*ksm(a[m]*a[n-m],p-2,p)%p; /*求(a[m]*a[n-m])在(mod p)意义下的乘法逆元*/
										  /*拓展欧几里得与费马小定理均可*/ 
	/*LL x,y;
	exgcd(a[m]*a[n-m],p,x,y);
	return (a[n]*x%p+p)%p;*/ 
}
LL lucas(LL n,LL m,LL p)
{
	if(!m)return 1;
	return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
int main()
{
	ios::sync_with_stdio(false);
	LL T,n,m,p;
	cin>>T;
	while(T--){
		cin>>n>>m>>p;
		init(p);
		cout<<lucas(n+m,m,p)<<endl;
	}
	return 0;
}

为了快,放了一个数组,限制了P的大小,然后再发一个没有限制的卢卡斯。

#define MOD 1000000007
typedef long long LL;

LL quickPower(LL a, LL b)
{
    LL ans = 1;
    a %= MOD;
    while (b)
    {
        if (b & 1)
        {
            ans = ans * a % MOD;
        }
        b >>= 1;
        a = a * a % MOD;
    }
    return ans;
}

LL c(LL n, LL m)
{
    if (m > n)
    {
        return 0;
    }
    LL ans = 1;
    for (int i = 1; i <= m; i++)
    {
        LL a = (n + i - m) % MOD;
        LL b = i % MOD;
        ans = ans * (a * quickPower(b, MOD - 2) % MOD) % MOD;
    }
    return ans;
}

LL lucas(LL n, LL m)
{
    if (m == 0)
    {
        return 1;
    }
    return c(n % MOD, m % MOD) * lucas(n / MOD, m / MOD) % MOD;
}

int main(int argc, const char * argv[])
{
    LL n, m;
    while (~scanf("%lld %lld", &n, &m))
    {
        LL max, min;
        max = n + m - 3;
        min = m - 1;
        printf("%lld\n", lucas(max - 1, min - 1));
    }
    return 0;
}

ExLucas扩展卢卡斯定理: 求解 C n m % P C_{n}^{m}\% P Cnm​%P,其中m,n较大,P较小且不一定为素数 转化为CRT(中国剩余定理) 好像这也不是什么定理,只是一个计算方法 计 算 C n m % P , 其 中 p = p 1 q 1 × p 2 q 2 × ⋯ p k q k 时 , 我 们 可 以 先 求 出 C m n   m o d   p i q i , 然 后 用 C R T 合 并 。 那 么 怎 么 计 算 C m n   m o d   p i q i 呢 ? C m n = m ! n ! ( m − n ) ! , 我 们 只 需 要 算 出 m ! , n ! − 1 , ( m − n ) ! − 1 , 然 后 乘 在 一 起 。 注 : n ! 可 能 在 模 p i q i 的 意 义 下 没 有 逆 元 啊 , 那 这 就 是 错 的 了 其 实 这 里 求 得 不 是 逆 元 ( 可 能 没 有 逆 元 ) , 求 出 来 的 是 a × p i b ( g c d ( a , p ) = 1 ) , 前 面 的 a a 用 逆 元 , 后 面 的 次 数 加 加 减 减 一 下 就 好 了 问 题 转 换 成 求 n !   m o d   p q 例 如 n = 19 , p = 3 , q = 2 n = 19 , p = 3 , q = 2 : 计算C_{n}^{m}\% P,其中p=p1^{q1}×p2^{q2}×⋯pk^{qk}时,我们可以先求出C_m^n~mod~{p_i}^{q_i},然后用CRT合并。\\ 那么怎么计算C_m^n~mod~{p_i}^{q_i}呢?\\ C_m^n=\frac{m!}{n!(m-n)!},我们只需要算出m!,n!^{−1},(m−n)!^{−1} ,然后乘在一起。\\ 注:n! 可能在模{p_i}^{q_i}的意义下没有逆元啊,那这就是错的了\\ 其实这里求得不是逆元(可能没有逆元),求出来的是a\times {p_i}^b(gcd(a,p)=1),\\前面的aa用逆元,后面的次数加加减减一下就好了\\ 问题转换成求n!~mod~p^q\\ 例如n=19,p=3,q=2n=19,p=3,q=2: 计算Cnm​%P,其中p=p1q1×p2q2×⋯pkqk时,我们可以先求出Cmn​ mod pi​qi​,然后用CRT合并。那么怎么计算Cmn​ mod pi​qi​呢?Cmn​=n!(m−n)!m!​,我们只需要算出m!,n!−1,(m−n)!−1,然后乘在一起。注:n!可能在模pi​qi​的意义下没有逆元啊,那这就是错的了其实这里求得不是逆元(可能没有逆元),求出来的是a×pi​b(gcd(a,p)=1),前面的aa用逆元,后面的次数加加减减一下就好了问题转换成求n! mod pq例如n=19,p=3,q=2n=19,p=3,q=2: 参考 19 ! = 1 × 2 × 3 × ⋯ × 19 = ( 1 × 2 × 4 × 5 × 7 × 8 ⋯ × 16 × 17 × 19 ) × ( 3 × 6 × 9 × 12 × 15 × 18 ) = ( 1 × 2 × 4 × 5 × 7 × 8 ⋯ × 16 × 17 ) × 19 × 36 × ( 1 × 2 × 3 × 4 × 5 × 6 ) = ( 1 × 2 × 4 × 5 × 7 × 8 ) 2 × 19 × 36 × ( 1 × 2 × 3 × 4 × 5 × 6 ) 19!\\ =1×2×3×⋯×19\\ =(1×2×4×5×7×8⋯×16×17×19)×(3×6×9×12×15×18)\\ =(1×2×4×5×7×8⋯×16×17)×19×36×(1×2×3×4×5×6)\\ =(1×2×4×5×7×8)2×19×36×(1×2×3×4×5×6) 19!=1×2×3×⋯×19=(1×2×4×5×7×8⋯×16×17×19)×(3×6×9×12×15×18)=(1×2×4×5×7×8⋯×16×17)×19×36×(1×2×3×4×5×6)=(1×2×4×5×7×8)2×19×36×(1×2×3×4×5×6)

上面这个式子分为四部分: 第 一 部 分 : ( 1 × 2 × 4 × 5 × 7 × 8 ) 2 ( 1 × 2 × 4 × 5 × 7 × 8 ) 2 。 这 部 分 的 数 不 超 过 p q p q 个 , 可 以 暴 力 算 第 二 部 分 : 1919 。 这 部 分 的 数 不 超 过 p q p q 个 , 可 以 暴 力 算 第 三 部 分 : 3636 。 这 个 在 最 后 处 理 时 求 出 m ! , n ! , ( m − n ) ! m ! , n ! , ( m − n ) !              分 别 有 多 少 个 p p ( 设 为 x , y , z x , y , z ) , 则 答 案 要 乘 上 p x − y − z p x − y − z 第 四 部 分 : 1 × 2 × 3 × 4 × 5 × 61 × 2 × 3 × 4 × 5 × 6 。                      这 个 是 ⌊ n p ⌋ ! ⌊ n p ⌋ ! , 可 以 递 归 处 理 第一部分:(1×2×4×5×7×8)2(1×2×4×5×7×8)2。这部分的数不超过pqpq个,可以暴力算\\ 第二部分:1919。这部分的数不超过pqpq个,可以暴力算\\ 第三部分:3636。这个在最后处理时求出m!,n!,(m−n)!m!,n!,(m−n)!\\ \qquad\ \ \ \ \ \ \ \ \ \ \ \ 分别有多少个pp(设为x,y,zx,y,z),则答案要乘上px−y−zpx−y−z\\ 第四部分:1×2×3×4×5×61×2×3×4×5×6。\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 这个是⌊np⌋!⌊np⌋!,可以递归处理 第一部分:(1×2×4×5×7×8)2(1×2×4×5×7×8)2。这部分的数不超过pqpq个,可以暴力算第二部分:1919。这部分的数不超过pqpq个,可以暴力算第三部分:3636。这个在最后处理时求出m!,n!,(m−n)!m!,n!,(m−n)! 分别有多少个pp(设为x,y,zx,y,z),则答案要乘上px−y−zpx−y−z第四部分:1×2×3×4×5×61×2×3×4×5×6。 这个是⌊np⌋!⌊np⌋!,可以递归处理

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N=1e5+9;
LL A[N],M[N];
LL ksm(LL x,LL n,LL mod)
{
	LL ans=1;
	while(n){
		if(n&1)ans=ans*x%mod;
		n>>=1,x=x*x%mod;
	}
	return ans;
}
void exgcd(LL a,LL b,LL &x,LL &y)
{
	if(!b)x=1,y=0;
	else exgcd(b,a%b,y,x),y-=a/b*x;
}
LL inv(LL a,LL p)
{
	LL x,y;
	exgcd(a,p,x,y);
	return (x+p)%p?x:x+p;
}
LL get(LL n,LL pi,LL p)									/*求(与pi互素后的n!)%M[i]*/ 
{
	if(!n)return 1;
	LL ans=1;
	if(n/p){											/*判断有无循环节 */ 
		for(LL i=2;i<=p;++i)if(i%pi)ans=ans*i%p;
		ans=ksm(ans,n/p,p);
	}
	for(LL i=2;i<=n%p;++i)if(i%pi)ans=ans*i%p;			/*循环节剩余部分*/ 
	return ans*get(n/pi,pi,p)%p;
}
LL exlucas(LL n,LL m,LL pi,LL p)						/*求A[i]*/ 
{
	LL nn=get(n,pi,p);									/*求(与pi互素后的n)%M[i]*/ 
	LL mm=get(m,pi,p);									/*求(m!与pi互素后的m!)%M[i]*/ 
	LL nm=get(n-m,pi,p);								/*求(与pi互素后的(n-m)!)%M[i]*/ 
	LL k=0;												/*含质因数pi的数量*/ 
	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;
	return nn*inv(mm,p)*inv(nm,p)*ksm(pi,k,p)%p;
}
LL crt(LL len,LL Lcm)
{
	LL ans=0;
	for(LL i=1;i<=len;++i){
		LL Mi=Lcm/M[i];
		ans=((ans+A[i]*inv(Mi,M[i])*Mi)%Lcm+Lcm)%Lcm;
	}
	return ans;
}
int main()
{
	ios::sync_with_stdio(false);
	LL n,m,P,num;
	while(cin>>n>>m>>P){
		if(n<m){
			cout<<0<<endl;
			continue;
		}
		num=0;
		memset(A,0,sizeof(A));
		memset(M,0,sizeof(M));
		for(LL x=P,i=2;i<=P;++i)
			if(x%i==0){
				M[++num]=1;
				while(x%i==0){
					M[num]*=i;
					x/=i;
				}
				A[num]=exlucas(n,m,i,M[num])%P;
			} 
		cout<<crt(num,P)<<endl;
	}
	return 0;
}

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 数学--数论--组合数(卢卡斯+扩展卢卡斯)模板

    风骨散人Chiam
  • 数学--数论--中国剩余定理 拓展 HDU 1788

    问题描述 我知道部分同学最近在看中国剩余定理,就这个定理本身,还是比较简单的: 假设m1,m2,…,mk两两互素,则下面同余方程组: x≡a1(mod m...

    风骨散人Chiam
  • 数学--数论--中国剩余定理+扩展中国剩余定理(孙子定理)

    其中 m 1 ...

    风骨散人Chiam
  • 打表法——暴力破解方法之一

    特别声明:除特别标注,本站文章均为原创,本站文章原则上禁止转载,如确实要转载,请电联:wangyeuuu@qq.com,尊重他人劳动成果,谢过~

    lollipop72
  • POJ 2773 Happy 2006(容斥原理+二分)

    Happy 2006 Time Limit: 3000MS Memory Limit: 65536K Total Submissions: 1...

    ShenduCC
  • 11.6NOIP模拟赛解题报告

    很显然的一个贪心是从左往右扫,如果遇到一个不合法的点\(i\),那么升级\(i + R\)处的炮台。。

    attack
  • 洛谷P3799 妖梦拼木棒

    题目背景 上道题中,妖梦斩了一地的木棒,现在她想要将木棒拼起来。 题目描述 有n根木棒,现在从中选4根,想要组成一个正三角形,问有几种选法? 输入输出格式 输入...

    attack
  • 2017 Multi-University Training Contest - Team 9 1002&&HDU 6162 Ch’s gift【树链部分+线段树】

    Ch’s gift Time Limit: 6000/3000 MS (Java/Others)    Memory Limit: 65536/65536 K ...

    Angel_Kitty
  • 扩展中国剩余定理详解

    前言 阅读本文前,推荐先学一下中国剩余定理。其实不学也无所谓,毕竟两者没啥关系 扩展CRT 我们知道,中国剩余定理是用来解同余方程组 但是有一个非常令...

    attack
  • 洛谷P1516 青蛙的约会

    题目描述 两只青蛙在网上相识了,它们聊得很开心,于是觉得很有必要见一面。它们很高兴地发现它们住在同一条纬度线上,于是它们约定各自朝西跳,直到碰面为止。可是它们出...

    attack

扫码关注云+社区

领取腾讯云代金券