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


    组合数:
    在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();
    }
    

    卢卡斯定理:
    Cmn mod pC_m^n~mod~p
    m=a0p0+a1p1++akpkm={a_0}^{p_0}+{a_1}^{p_1}+cdots+{a_k}^{p_k}
    n=b0p0+b1p1++bkpkn={b_0}^{p_0}+{b_1}^{p_1}+cdots+{b_k}^{p_k}
    CmnCaibi(mod p)C_m^nequivprod{C_{a_i}^{b_i}}(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
    ", lucas(max - 1, min - 1));
        }
        return 0;
    }
    

    ExLucas扩展卢卡斯定理:
    求解Cnm%PC_{n}^{m}\% P,其中m,n较大,P较小且不一定为素数
    转化为CRT(中国剩余定理)
    好像这也不是什么定理,只是一个计算方法
    Cnm%Pp=p1q1×p2q2×pkqkCmn mod piqiCRTCmn mod piqiCmn=m!n!(mn)!m!,n!1,(mn)!1n!piqia×pib(gcd(a,p)=1)aan! mod pqn=19,p=3,q=2n=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 imes {p_i}^b(gcd(a,p)=1),\前面的aa用逆元,后面的次数加加减减一下就好了\ 问题转换成求n!~mod~p^q\ 例如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)

    上面这个式子分为四部分:
    (1×2×4×5×7×8)2(1×2×4×5×7×8)2pqpq1919pqpq3636m!,n!,(mn)!m!,n!,(mn)!            ppx,y,zx,y,zpxyzpxyz1×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)!\ 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⌋!,可以递归处理

    #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;
    }
    
  • 相关阅读:
    Angular之ngRoute与uiRoute
    Python实现网络测速--转载
    You are using pip version 10.0.1, however version 20.2.2 is available.
    Pycharm导出环境,实现环境迁移
    Python编辑器之pycharm2020.1破解汉化
    Rsyslog系统日志转发
    服务器DNS配置
    服务器系统配置初始化脚本
    系统巡检用到的参数
    inotify高效监控Linux文件系统事件
  • 原文地址:https://www.cnblogs.com/lunatic-talent/p/12798529.html
Copyright © 2020-2023  润新知