大组合数取模

最新推荐文章于 2019-03-01 16:38:00 发布
最新推荐文章于 2019-03-01 16:38:00 发布
阅读量994 收藏
点赞数
分类专栏： 数学
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/incincible/article/details/59117274
版权

考虑C(n,m)%P
情况一：n,m很大，P为素数
处理小范围的阶乘和阶乘的逆元 用卢卡斯定理即可。卢卡斯定理:
这里写图片描述

情况二： 当P=

p1∗p2∗p3∗...∗pn

求出 [Cmn] 分别在 [p1,p2,p3,…,pn] 模意义下的结果，记为 [m1,m2,m3,…,mn]

得到模方程组

xi≡mi(modpi)1≤i≤n
用中国剩余定理求解即可。
代码如下： 

#include<cstdio>
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<vector>
#include<algorithm>
#include<cmath>

#define LL long long

const LL M= 999911659;
const int maxn= 50000+5;

using namespace std;
template <typename T>
inline void _read(T &x){
    char ch=getchar(); bool mark=false;
    for(;!isdigit(ch);ch=getchar())if(ch=='-')mark=true;
    for(x=0;isdigit(ch);ch=getchar())x=x*10+ch-'0';
    if(mark)x=-x;
}

LL P[4]= {2,3,4679,35617};
LL N,G,fra[maxn][4],inv[maxn][4];
LL ans[4];

LL pow_mod(LL x,LL y,LL M){
    LL ans=1;
    for(x%=M;y;y>>=1,x=x*x%M)
        if(y&1) ans=ans*x%M;
    return ans;
}

// 解模方程组  x=A[i] (moc m[i]) 0<=i<n

void EXGCD(LL A,LL B,LL& d,LL& x,LL& y){
    if(B==0){d=A;x=1;y=0;}
    else {EXGCD(B,A%B,d,y,x); y-= x*(A/B);}
}

LL China(LL * A,LL * m,LL n){
    LL M= 1,ret=0,i;
    for(i=0;i<n;i++)M*= m[i];
    for(i=0;i<n;i++){
        LL d,x,y,t= M/m[i];
        EXGCD(m[i],t,d,x,y); 
        ret= (ret+ A[i]*y*t%M)%M;
    }
    return (ret+M)%M;
}

LL C(LL n,LL m,LL id){
    if(m>n) return 0; 
    if(!m||n==m) return 1;
    if(n<=P[id] && m<=P[id])
        return fra[n][id]*inv[m][id]%P[id]*inv[n-m][id]%P[id];
    return C(n/P[id],m/P[id],id)*C(n%P[id],m%P[id],id)%P[id];
}

void Divide(LL n){
    LL i,j;
    for(i=1;i*i<=n;i++){
        if(n%i!=0) continue;
        for(j=0;j<4;j++){
            ans[j]= ans[j]+C(n,i,j);
            if(ans[j]>=P[j]) ans[j]-=P[j];
        }
        if(i*i!=n){
            for(j=0;j<4;j++){
                ans[j]= ans[j]+ C(n,n/i,j);
                if(ans[j]>=P[j]) ans[j]-=P[j];
            }
        }
    }
}

void Init(){
    LL i,j;
    for(i=0;i<4;i++){   //分别预处理每一个质数  
        LL mod= P[i];
        inv[1][i]= fra[1][i]=1;
        inv[0][i]= fra[0][i]=1;
        for(j=2;j<=P[i];j++){
            fra[j][i]= fra[j-1][i]*j%mod;
            inv[j][i]= (mod-mod/j)*inv[mod%j][i]%mod;
        }
        for(j=2;j<=P[i];j++) inv[j][i]= inv[j][i]*inv[j-1][i]%mod;
    }
}

int main(){
    cin>>N>>G;
    Init(); 
    Divide(N);
    //for(int i=0;i<4;i++) cout<<ans[i]<<" ";cout<<endl;
    LL exp= China(ans,P,4);
    if(exp==0) exp= M-1;
    cout<<pow_mod(G,exp,M)<<endl;
    return 0;
}
/*
900951975 90523153
81793986
*/

终极情况：C(n,m)%p
n,m非常大 -> 卢卡斯定理
p=p1^c1 * p2^c2 * … * pn^cn
由之前的情况我们已经可以得到了一些结论，这种情况怎么办？
同理 卢卡斯定理 + CRT
但是前提条件是要求 C(n,m)% pi^ci
C(n,m)% pi^ci怎么求？
这里写图片描述

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cstring>
#define ll long long
using namespace std;
const ll label=1000000000000000LL;
template <typename T>
inline void _read(T& x){
    char t=getchar();bool sign=true;
    while(t<'0'||t>'9'){if(t=='-')sign=false;t=getchar();}
    for(x=0;t>='0'&&t<='9';t=getchar())x=x*10+t-'0';
    if(!sign)x=-x;
}
ll mont(ll x,ll y,ll z){
    ll temp=1;
    for(x%=z;y;y>>=1,x=x*x%z){
        if(y&1)temp=temp*x%z;
    }
    return temp;
}
ll exgcd(ll a,ll b,ll &x,ll &y){  
    if(b==0){  
        x=1;y=0;return a;  
    }
    else{  
        ll ans=exgcd(b,a%b,x,y);  
        ll t=x;  
        x=y;y=t-a/b*y;return ans;  
    }
}
ll rev(ll a,ll b){
    ll x,y;
    exgcd(a,b,x,y);
    return (x%b+b)%b;
}
ll C(ll n,ll m,ll mod){
    if(m>n)return 0;
    ll ans=1,i,a,b;
    for(i=1;i<=m;i++){
        a=(n+1-i)%mod;
        b=rev(i%mod,mod);
        ans=ans*a%mod*b%mod;
    }
    return ans;
}
ll C1(ll n,ll m,ll mod){
    if(m==0)return 1;
    return C(n%mod,m%mod,mod)*C1(n/mod,m/mod,mod)%mod;
}
ll cal(ll n,ll p,ll t){
    if(!n)return 1;
    ll x=mont(p,t,label),i,y=n/x,temp=1;
    for(i=1;i<=x;i++){
        if(i%p)temp=temp*i%x;
    }
    ll ans=mont(temp,y,x);
    for(i=y*x+1;i<=n;i++){
        if(i%p)ans=ans*i%x;
    }
    return ans*cal(n/p,p,t)%x;
}
ll C2(ll n,ll m,ll p,ll t){
    ll x=mont(p,t,label);
    ll a,b,c,ap=0,bp=0,cp=0,temp;
    for(temp=n;temp;temp/=p)ap+=temp/p;
    for(temp=m;temp;temp/=p)bp+=temp/p;
    for(temp=n-m;temp;temp/=p)cp+=temp/p;
    ap=ap-bp-cp;
    ll ans=mont(p,ap,x);
    a=cal(n,p,t);
    b=cal(m,p,t);
    c=cal(n-m,p,t);
    ans=ans*a%x*rev(b,x)%x*rev(c,x)%x;
    return ans;
}
ll CRT(ll A[],ll B[],ll r){
    ll M=1;
    ll i,j,k,ans=0,d,t,x0,y0;
    for(i=1;i<=r;i++){
        M*=A[i];
    }
    for(i=1;i<=r;i++){
        d=M/A[i];
        exgcd(d,A[i],x0,y0);
        ans=(ans+d*x0*B[i])%M;
    }
    while(ans<=0)ans+=M;
    return ans;
}
ll lucas(ll x,ll y,ll mod){
    ll i,j,k,d,cnt=0,t;
    ll A[205],M[205];
    for(i=2;i*i<=mod;i++){
        if(mod%i==0){
            t=0;
            while(mod%i==0){
                mod/=i;
                t++;
            }
            M[++cnt]=mont(i,t,label);
            if(t==1)A[cnt]=C1(x,y,i);
            else A[cnt]=C2(x,y,i,t);
        }
    }
    if(mod>1){
        M[++cnt]=mod;
        A[cnt]=C1(x,y,mod);
    }
    /*cout<<"divide:"<<endl;
    for(i=1;i<=cnt;i++){
        cout<<"M:"<<M[i]<<" residue:"<<A[i]<<endl;
    }*/
    return CRT(M,A,cnt);
}
int T;
int main(){
    ll x,y,temp;
    cin>>T;
    while(T--){
        scanf("%I64d%I64d%I64d",&x,&y,&temp);
        printf("%I64d\n",lucas(x,y,temp));
    }
}

#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<algorithm>
#include<cstring>
#define ll long long
using namespace std;
const ll label=1000000000000000LL;
template <typename T>
inline void _read(T& x){
    char t=getchar();bool sign=true;
    while(t<'0'||t>'9'){if(t=='-')sign=false;t=getchar();}
    for(x=0;t>='0'&&t<='9';t=getchar())x=x*10+t-'0';
    if(!sign)x=-x;
}
int T;
ll mod;
ll cnt;
ll n,m,n1,n2,c;
ll a[25];
ll A[205],M[205],prime[205],index[205];
ll F[50005];

ll mont(ll x,ll y,ll z){
    ll temp=1;
    for(x%=z;y;y>>=1,x=1ll*x*x%z){
        if(y&1)temp=1ll*temp*x%z;
    }
    return temp;
}
void pre(){
    ll MOD=mod,i,t;
    for(i=2;i*i<=MOD;i++){
        if(MOD%i==0){
            t=0;
            prime[++cnt]=i;
            while(MOD%i==0){
                MOD/=i;
                t++;
            }
            M[cnt]=mont(i,t,label);
            index[cnt]=t;
        }
    }
    if(MOD>1){
        M[++cnt]=MOD;
        prime[cnt]=MOD;
        index[cnt]=1;
    }
}
ll exgcd(ll a,ll b,ll &x,ll &y){  
    if(b==0){  
        x=1;y=0;return a;  
    }
    else{  
        ll ans=exgcd(b,a%b,x,y);  
        ll t=x;  
        x=y;y=t-a/b*y;return ans;  
    }
}

ll _fac(ll n,ll t){
    if(n<prime[t])return F[n];
    c+=n/prime[t];
    return 1ll*mont(F[M[t]-1],n/M[t],M[t])*F[n%M[t]]%M[t]*_fac(n/prime[t],t)%M[t];
}
ll C(ll n,ll m,ll t){
    if(n<m)return 0;
    if(m<0)return 0;
    if(n==m)return 1;
    ll i,p,temp,x,y,s;
    F[0]=1;
    for(i=1;i<M[t];i++){
        if(i%prime[t])F[i]=F[i-1]*i%M[t];
        else F[i]=F[i-1];
    }
    c=0;
    p=_fac(n,t);
    s=c;c=0;
    temp=1ll*_fac(m,t)*_fac(n-m,t)%M[t];
    exgcd(temp,M[t],x,y);
    x=(x+M[t])%M[t];
    p=1ll*p*x%M[t]*mont(prime[t],s-c,M[t])%M[t];
    return p;
}
ll lucas(ll n,ll m){
    if(n<m||m<0)return 0;
    if(n==m)return 1;

    //cout<<"lucas("<<n<<","<<m<<")"<<endl;
    ll temp=0,x,y,i,t;
    for(i=1;i<=cnt;i++){
        //cout<<"i:"<<i<<" prime:"<<prime[i]<<" mod:"<<M[i]<<" index:"<<index[i]<<endl;
        t=mod/M[i];
        exgcd(t,M[i],x,y);x=(x+M[i])%M[i];
        temp=(temp+1ll*x*C(n,m,i)%mod*t%mod)%mod;
    }
    return temp;
}
int main(){
    ll i,j,k;
    cin>>mod;
    pre();
}
INCINCIBLE
关注
  • 0
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
组合数学- 组合数取模- 逆元与递推打表.rar
09-16
组合数学- 组合数取模- 逆元与递推打表.rar
(扩展)Lucas定理——求大组合数取模
CaptainChen的博客
03-31 282
Lucas定理 若PPP为质数 Crn≡C⌊r/P⌋⌊n/P⌋×Cr&amp;amp;nbsp;mod&amp;amp;nbsp;Pn&amp;amp;nbsp;mod&amp;amp;nbsp;P&amp;amp;nbsp;(mod&amp;amp;nbsp;P)Cnr≡C⌊n/P⌋⌊r/P⌋×Cn&amp;amp;nbsp;mod&amp;amp;nbsp;Pr&amp;amp;nbsp;mod&amp;amp;nbsp
Lucas定理（大组合数的取模）
随心丶而遇的博客
09-15 1943
首先给出这个Lucas定理： A、B是非负整数，p是质数。AB写成p进制：A=a[n]a[n-1]...a[0]，B=b[n]b[n-1]...b[0]。 则组合数C(A,B)与C(a[n],b[n])*C(a[n-1],b[n-1])*...*C(a[0],b[0])  mod p同余 即：Lucas(n，m，p)=c(n%p，m%p)*Lucas(n/p，m/p，p)  这个定理的
Lucas定理——大组合数取模
AnnieL
02-22 488
大组合数取模，求C[n][m]%p 公式：C[n][m]%p == C[n%p][m%p]*C[n/p][m/p]%p 代码（内含粗略解释）： #define LL long long LL mont(LL a,LL b,LL c){//快速幂 LL ans=1; a%=c; while(b){ if(b&amp;amp;amp;amp;amp;amp;1)ans=ans*a%c; ...
教你一波Lucas(卢卡斯)定理在数论解题中的应用
LzyRapX
08-05 1万+
首先，Lucas（卢卡斯）定理是什么？有什么用？ Lucas定理是用来求 C(n,m) mod p，p为素数的值。（注意：p一定是素数） 有人会想，C(n,m)不能用C(n, m) = C(n - 1，m) + C(n - 1, m - 1)的公式来递推吗？ （ 提示：C(n, m) mod p = n!/(m!(n - m)!) mod p ） 可以是可以。但当n
组合数学- 组合数取模.rar
09-16
组合数学- 组合数取模.rar
利用扩展欧几里得算法快速计算组合数取模
08-17
利用扩展欧几里得方法，进行快速的计算C（n,m）%P;
组合数学- 组合数取模- 卢卡斯定理与扩展卢卡斯定理.rar
09-16
组合数学- 组合数取模- 卢卡斯定理与扩展卢卡斯定理.rar
PCtoLCD2002_取模软件
10-30
把文字输入框换成了文字输入组合框，这样就可以保存历史纪录。 5.输出数据会自动清除以前的数据 6.可隐藏自定义格式，简化操作 7.一些小的BUG修正。 版本为pctolcd2.03 5月3日凌晨发布第三测试版 更新如下: 1...
组合数取模总结
Hang_ccccc的博客
12-07 2819
题目描述： 在ACM竞赛现在叫JB竞赛？中，经常会遇到组合数取模的题目；就我现在的水平而言，大概分为以下三类，以后遇到新的方法会在做补充； 第一种： n和m都较小 （&lt;1000）,在这个数据范围内，我们可以直接用杨辉三角O(n^2)的复杂度内预处理出所有的组合数，然后直接输出即可；典型例题就是给你一个二项式，例如: 这样子的题目，组合数我们可以O(n2)预处理出来（利用递推式子来处理）...
大组合数取模——卢卡斯定理
dianshu1593的博客
03-01 213
我们学了O(n^2)的做法，加上逆元，我们又会了O(n)的做法， 在来了新问题，如果n和m很大呢， 比如求C(n, m) % p ，n<=1e18,m<=1e18,p<=1e5 看到没有，n和m这么大，但是p却很小，我们要利用这个p。 接下来进入正题：Lucas定理针对该取值范围较大又不特别大的情况 定理描述 ： 这样将组合数的求解分解为小问...
组合数取模
热门推荐
ACdreamer
10-03 3万+
Lucas定理： m，n是非负整数，p是质数。m，n分别写成p进制： 则组合数C(n,m)与C(n[k],m[k])*C(n[k-1],m[k-1])*...*C(n[0],m[0]) mod p同余  除了这个Lucas定理以外，还要用到快速幂和另外一个不知名定理： 若p为质数，则(a/b)(mod p) = a * b^(p-2) (mod p)。如果p不为质数就不能这样求。 a
组合(Lucas)
ivyChen的专栏
07-14 564
Description 给出组合数C(n,m), 表示从n个元素中选出m个元素的方案数。例如C(5,2) = 10, C(4,2) = 6.可是当n,m比较大的时候，C(n,m)很大！于是xiaobo希望你输出 C(n,m) mod p的值！ Input 输入数据第一行是一个正整数T,表示数据组数 (T Output 对于每组数据，输出一个正整数，表
Lucas定理推广（组合数取模）
guanhuaing
08-20 463
Lucas定理推广（组合数取模） 2011年4月21日Earthson3 条评论 比赛的时候，意外想到了个方法来处理组合数取模。一推，发现就是Lucas定理，并且，在模数是素数乘方的情况作了些推广。使得lucas能够处理模任意数的组合数取模。 Lucas定理 首先给出Lucas定理的递归描述 简要证明思路： 考虑连续的p...
大数组合数取模(逆元+打表)
Daniel
05-23 2767
将阶乘O(n)打表之后C(n,m)便可O(1)求出，除法取模用逆元解决 hdu5698瞬间移动#include<bits/stdc++.h> using namespace std; const int maxn=200000; const int p=1000000007;int f[maxn+10]; void init(){//阶乘 f[0]=1; for(int i=1;i
组合数模板（取模，大数）
wrwhahah的博客
08-18 1372
#include &lt;iostream&gt; #include &lt;bits/stdc++.h&gt; #define maxn 200005 typedef long long ll; using namespace std; const ll mod=998244353; ll fac[maxn],inv[maxn]; ll pow_mod(ll a,ll n) { ll ...
各种逆元求法 组合数取模 comb （组合数 Lucas）
w4149
06-02 4419
组合数取模（comb）【问题描述】 计算C（m,n）mod 9901的值 【输入格式】 从文件comb.in中输入数据。 输入的第一行包含两个整数，m和n 【输出格式】 输出到文件comb.out中。 输出一行，一个整数 【样例输入】 2 1 【样例输出】 2【数据规模与约定】 对于 20%的数据，n<=m<=20 对于 40%的数据，n<=m<=2000 对于 100%
大组合数取模模板
RCY_ZHU的博客
08-13 880
LL n,m,p = 1e9+7;LL quick_mod(LL a, LL b) { LL ans = 1; a %= p; while(b) { if(b & 1) { ans = ans * a % p; b--; } b >>= 1;
FZU - 2020 计算大组合数取模
sxy201658506207的博客
10-30 285
题意: 给出m,n. 计算C(n,m)%p,C(n,m)很大,p是素数  卢卡斯.不过不可以预处理,mod太大了 #include &lt;iostream&gt;//Lucas模板 #include &lt;cstdio&gt; #include &lt;algorithm&gt; #include &lt;cstring&gt; #include &lt;cmath&gt; #i...
快速幂求组合数python
最新发布
07-09
要使用快速幂求组合数，可以利用组合数的性质：C(n, k) = n! / (k! * (n-k)!)。具体的实现如下： ```python def binomial_coefficient(n, k): # 求阶乘 def factorial(num): result = 1 for i in range(1, num ...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值