从欧几里得算法到组合数取模

来源：互联网发布：怎么开好淘宝店? 编辑：程序博客网时间：2024/06/05 19:36

欧几里得算法→扩展欧几里得算法→求模运算下的乘法逆元→大组合数取模

1、欧几里得算法-----gcd(a,b)

辗转相除法无需多言

#include<stdio.h>typedef long long LL;LL gcd(LL a,LL b)         //欧几里得辗转相除法{    return (b==0)? a : gcd(b,a%b);}

2、拓展欧几里得算法-----extend_gcd(a,b,d,x,y)

在求出gcd的基础上，求整数x,y，使得ax+by=gcd(a,b)

只需在gcd函数上做些改动，考虑gcd递归上下层之间的关系，从更深一层返回回来的x'和y'满足那一层的贝祖定理，即

本层： ax+by = gcd(a,b)=d

更深一层：bx'+(a%b)y '= gcd(b,a%b)=d

联立之，由于ab任意，得到上下层x，y的关系

x=y'

y=x'-(a/b)y'

于是代码如下：

#include<stdio.h>typedef long long LL;void extend_gcd(LL a,LL b,LL &d,LL &x,LL &y)   //扩展欧几里得算法，求整数x,y，使得ax+by=d{    if(b==0){d=a;x=1;y=0;}    else    {        extend_gcd(b,a%b,d,x,y);        LL t=x;        x=y;                  //可以根据下一层推出这一层的x，y值（联立ax+by=d 和 bx'+(a%b)y'=d 即可）        y=t-(a/b)*(y);    }}

3、求模运算下的乘法逆元-----inv(a,b)

主要是利用上述的拓展欧几里得算法

注意只有当gcd(a,b)==1时，a在模b下才存在乘法逆元！

ax+by=1 对b取模 → ax=1（mod b）

故调用extend——gcd 求得的x，即为a在b下的乘法逆元

#include<stdio.h>typedef long long LL;void extend_gcd(LL a,LL b,LL &d,LL &x,LL &y)   //扩展欧几里得算法，求整数x,y，使得ax+by=d{    if(b==0){d=a;x=1;y=0;}    else    {        extend_gcd(b,a%b,d,x,y);        LL t=x;        x=y;                  //可以根据下一层推出这一层的x，y值（联立ax+by=d 和 bx'+(a%b)y'=d 即可）        y=t-(a/b)*(y);    }}LL inv(LL a,LL b)        //求a在模b下的乘法逆元（注意只有在a,b互素时才有乘法逆元）,若不存在逆元则返回-1{    LL x,y,d;    extend_gcd(a,b,d,x,y);    return (d==1)? (x+b)%b : -1;  //因为x可能是负数，所以+b保证为正}

4、大组合数取模------C(n,m)

如果n,m都在1e6以内，可以先打一个阶乘表，组合数的计算直接按阶乘公式即可，当然这里同样要用到求逆元

#include<stdio.h>typedef long long LL;#define mod 1000000007LL table[1000005];void cal_table()      //计算阶乘表{    table[0]=1;    for(int i=1;i<=1000000;i++)        table[i]=table[i-1]*i%mod;}void extend_gcd(LL a,LL b,LL &d,LL &x,LL &y)   //扩展欧几里得算法，求整数x,y，使得ax+by=d{    if(b==0){d=a;x=1;y=0;}    else    {        extend_gcd(b,a%b,d,x,y);        LL t=x;        x=y;                  //可以根据下一层推出这一层的x，y值（联立ax+by=d 和 bx'+(a%b)y'=d 即可）        y=t-(a/b)*(y);    }}LL inv(LL a,LL b)        //求a在模b下的乘法逆元（注意只有在a,b互素时才有乘法逆元）,若不存在逆元则返回-1{    LL x,y,d;    extend_gcd(a,b,d,x,y);    return (d==1)? (x+b)%b : -1;  //因为x可能是负数，所以+b保证为正}LL C(LL n,LL m)       //组合数{    if(n<0 || m<0 || m>n) return 0;    return table[n]*inv(table[m],mod)%mod*inv(table[n-m],mod)%mod;}int main(){    cal_table();    LL n=10,m=3;    printf("%lld\n",C(n,m));    return 0;}

练习题：HDU - 5894

题意：n个座位，m个人，每个人之间至少隔k个位置。

解法：列出答案的组合数学表达式，用大组合数取模模板求解。（注意单独讨论m==1的情况）

#include<stdio.h>typedef long long LL;#define mod 1000000007LL table[1000005];void cal_table()      //计算阶乘表{    table[0]=1;    for(int i=1;i<=1000000;i++)        table[i]=table[i-1]*i%mod;}void extend_gcd(LL a,LL b,LL &d,LL &x,LL &y)   //扩展欧几里得算法，求整数x,y，使得ax+by=d{    if(b==0){d=a;x=1;y=0;}    else    {        extend_gcd(b,a%b,d,x,y);        LL t=x;        x=y;                  //可以根据下一层推出这一层的x，y值（联立ax+by=d 和 bx'+(a%b)y'=d 即可）        y=t-(a/b)*(y);    }}LL inv(LL a,LL b)        //求a在模b下的乘法逆元（注意只有在a,b互素时才有乘法逆元）,若不存在逆元则返回-1{    LL x,y,d;    extend_gcd(a,b,d,x,y);    return (d==1)? (x+b)%b : -1;  //因为x可能是负数，所以+b保证为正}LL C(LL n,LL m)       //组合数{    if(n<0 || m<0 || m>n) return 0;    return table[n]*inv(table[m],mod)%mod*inv(table[n-m],mod)%mod;}int main(){    cal_table();    int T,n,m,k;    LL ans;    scanf("%d",&T);    while(T--)    {        scanf("%d%d%d",&n,&m,&k);        ans = (m==1)? n : n*C(n-k*m-1,m-1)%mod*inv(m,mod)%mod;        printf("%lld\n",ans);    }    return 0;}

1 0