组合数取模就是求的值。


这个问题比较简单,组合数的计算可以靠杨辉三角,那么由于n和m的范围小,直接两层循环即可。

c[0][0]=1;
for(long long i=1;i<=5003;i++)//j是上面,i是下面
{
    c[i][0]=1;
    for(long long j=1;j<=i;j++)
    {
        c[i][j]=(c[i-1][j]+c[i-1][j-1])%MOD;
        //printf("%lld %lld %lld\n",i,j,c[i][j]);
    }
}

,并且p是素数。

这个问题有个叫做Lucas的定理。

题意:,其中,并且p是素数。

long long PowerMod(long long a,long long b,long long c)
{
    long long ans=1;
    long long k;
    k=a;
    k=k%c;
    while(b>0)
    {
        if(b&1)
        ans=(ans*k)%c;
        b>>=1;
        k=(k*k)%c;
    }
    return ans;
}
long long C(long long n,long long m,long long p)
{
    if(m>n)return 0;
    long long ans = 1;
    for(int i=1;i<=m;i++)
    {
        long long a=(n+i-m)%p;
        long long b=i%p;
        ans=ans*(a*PowerMod(b,p-2,p)%p)%p;
    }
    return ans;
}
long long Lucas(long long n,long long m,long long p)
{
    if(m==0)return 1;
    return C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p;
}
int main()
{
    long long n,m;
    while(scanf("%lld%lld",&n,&m)!=EOF)
    {
        printf("%lld\n",Lucas(n,m,MOD));   //n为下面的数,m为上面的数
    }
    return 0;
}

,并且p可能为合数。

这样的话先采取暴力分解,然后快速幂即可。

#include <iostream>
#include <string.h>
#include <stdio.h>

using namespace std;
typedef long long LL;
const int N = 200005;

bool prime[N];
int p[N];
int cnt;

void isprime()
{
    cnt = 0;
    memset(prime,true,sizeof(prime));
    for(int i=2; i<N; i++)
    {
        if(prime[i])
        {
            p[cnt++] = i;
            for(int j=i+i; j<N; j+=i)
                prime[j] = false;
        }
    }
}

LL quick_mod(LL a,LL b,LL m)
{
    LL ans = 1;
    a %= m;
    while(b)
    {
        if(b & 1)
        {
            ans = ans * a % m;
            b--;
        }
        b >>= 1;
        a = a * a % m;
    }
    return ans;
}

LL Work(LL n,LL p)
{
    LL ans = 0;
    while(n)
    {
        ans += n / p;
        n /= p;
    }
    return ans;
}

LL Solve(LL n,LL m,LL P)
{
    LL ans = 1;
    for(int i=0; i<cnt && p[i]<=n; i++)
    {
        LL x = Work(n, p[i]);
        LL y = Work(n - m, p[i]);
        LL z = Work(m, p[i]);
        x -= (y + z);
        ans *= quick_mod(p[i],x,P);
        ans %= P;
    }
    return ans;
}

int main()
{
    int T;
    isprime();
    cin>>T;
    while(T--)
    {
        LL n,m,P;
        cin>>n>>m>>P;
        n += m - 2;
        m--;
        cout<<Solve(n,m,P)<<endl;
    }
    return 0;
}

逆元求取组合数

由定理可知:如果用C(n, r)表示n-元素集的r-组合的个数,有

C(n,r)=n!r!∗(n−r)!

而我们的目标就是计算 C(n,r)%mod 的值。

由数论的知识我们知道,模运算的加法,减法,乘法和四则运算类似,即:
模运算与基本四则运算有些相似,但是除法例外。其规则如下:

  • (a + b) % p = (a % p + b % p) % p
  • (a - b) % p = (a % p - b % p) % p
  • (a b) % p = (a % p b % p) % p

但对于除法却不成立,即(a / b) % p ≠ (a % p / b % p) % p 。

显然数学家们是不能忍受这种局面的,他们扔出了“逆元”来解决这个问题。那么什么是逆元? 逆元和模运算中的除法又有说明关系呢?

首先给出数论中的解释:

对于正整数 a 和 p,如果有 ax≡1(modp),那么把这个同余方程中 x 的最小正整数解叫做 a 模 p 的逆元。

什么意思呢? 就是指,如果 ax%p=1 , 那么x的最小正整数解就是 a 的逆元。

现在我们来解决模运算的除法问题。假设

ab

同时存在另一个数 x 满足

ax%p=m

由模运算对乘法成立,两边同时乘以 b ,得到:

a%p=(m(b%p))%p

如果 a 和 b 均小于模数 p 的话,上式可以改写为:

a=bm%p

等式两边再同时乘以 x, 得到:

ax%p=m%p=xbm%p

因此可以得到:

bx%p=1

哎,x是b的逆元呀(x 在模运算的乘法中等同于 1b, 这就是逆元的意义)

由以上过程我们看到,求取 (ab%p) 等同于 求取 (a∗(b的逆元)%p) 。 因此,求模运算的除法问题就转化为就一个数的逆元问题了。

而求取一个数的逆元,有两种方法

  1. 拓展欧几里得算法
  2. 费马小定理

对于利用拓展欧几里得算法求逆元,很显然,如果bx%p=1,那么 bx+py=1, 直接利用 exgcd(b, p, x, y)(代码实现在后面给出),则 (x%p+p)%p 即为 b 的逆元。

void init()
{
    fac[0]=fac[1]=1;
    for(int i=2;i<=MAXN;i++)
        fac[i]=fac[i-1]*i%MOD;
    inv[MAXN]=PowerMod(fac[MAXN],MOD-2,MOD);
    for(int i=MAXN-1;i>=0;i--)
        inv[i]=inv[i+1]*(i+1)%MOD;
}
ll calc(int n,int m)
{
    if(n<m)return 0;
    return fac[n]*inv[m]%MOD*inv[n-m]%MOD;
}

对于第二种方法,因为在算法竞赛中模数p总是质数,所以可以利用费马小定理

bp−1%p=1

可以直接得到 b 的逆元是 bp−2 , 使用 快速幂 求解即可。

明白了以上几个关键点,那么求取组合数 C(n,r) 的算法就呼之欲出了:

  1. 求取1到n的阶乘对 mod 取模的结果存入数组 JC[] 中;
  2. 求取 C(n,r) 时, 先利用“拓展欧几里得算法”或者“费马小定理+快速幂”求 JC[r]的逆元存入临时变量 x1 ;
  3. 然后计算JC[n]∗x1%mod 存入临时变量 x2;(x2 即为n!r!%mod 的值)
  4. 求取JC[n - r] 的逆元存入临时变量 x3;
  5. 则可以得到 C(n,r)=x2∗x3%mod

下面是方法二的代码片段:

typedef long long LL;
const LL maxn(1000005), mod(1e9 + 7);
LL Jc[maxn];
 
void calJc()    //求maxn以内的数的阶乘
{
   Jc[0] = Jc[1] = 1;
   for(LL i = 2; i < maxn; i++)
       Jc[i] = Jc[i - 1] * i % mod;
}
/*
//拓展欧几里得算法求逆元
void exgcd(LL a, LL b, LL &x, LL&y)    //拓展欧几里得算法
{
   if(!b) x = 1, y = 0;
   else
    {
       exgcd(b, a % b, y, x);
       y -= x * (a / b);
    }
}
 
LL niYuan(LL a, LL b)   //求a对b取模的逆元
{
   LL x, y;
   exgcd(a, b, x, y);
   return (x + b) % b;
}
*/
 
//费马小定理求逆元
LL pow(LL a, LL n, LL p)    //快速幂 a^n % p
{
   LL ans = 1;
   while(n)
    {
       if(n & 1) ans = ans * a % p;
       a = a * a % p;
       n >>= 1;
    }
   return ans;
}
 
LL niYuan(LL a, LL b)   //费马小定理求逆元
{
   return pow(a, b - 2, b);
}
 
LL C(LL a, LL b)    //计算C(a, b)
{
   return Jc[a] * niYuan(Jc[b], mod) % mod
       * niYuan(Jc[a - b], mod) % mod;
}

本博客所有文章除特别声明外,均采用 CC BY-SA 4.0 协议 ,转载请注明出处!

数位DP Previous
数论相关 Next