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

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

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 1000000007

LL 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 1000000007

LL 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;
}