FZU 2020（LUCAS）@-CSDN博客

本文链接：https://blog.csdn.net/yjf3151731373/article/details/52134924

组合

Time Limit: 1000MS

Memory Limit: 32768KB

64bit IO Format: %I64d & %I64u

Submit Status

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 <= 100) 接下来是T组数据，每组数据有3个正整数 n, m, p (1 <= m <= n <= 10^9, m <= 10^4, m < p < 10^9, p是素数)

Output

对于每组数据，输出一个正整数，表示C(n,m) mod p的结果。

Sample Input

2
5 2 3
5 2 61

Sample Output

1
10

组合数取模就是求的值，根据，和的取值范围不同，采取的方法也不一样。
下面，我们来看常见的两种取值情况（m、n在64位整数型范围内）
（1）  ， 
     此时较简单，在O(n²)可承受的情况下组合数的计算可以直接用杨辉三角递推，边做加法边取模。
（2） ，   ，并且是素数
 　本文针对该取值范围较大又不太大的情况（2）进行讨论。
     这个问题可以使用Lucas定理，定理描述：
　　　　　　
　其中
     
     ^{这样将组合数的求解分解为小问题的乘积，下面考虑计算}C(ni, mi) %p.
^{已知C(n, m)} ^{mod p = n!/(m!(n - m)!) mod p。当我们要求（a/b）mod p的值，且a很大，无法直接求得a/b的值时，我们可以转而使用乘法逆元k，将a乘上k再模p，即(a*k) mod p。 其结果与(a/b) mod p等价。}
^{那么逆元是什么？}
<span style="font-size: 18px;">定义：满足a*k≡1 (mod p)的k值就是a关于p的乘法<strong>逆元</strong>（当p是1时，对于任意a，k都为1）</span>
^{除法取模，这里要用到m!(n-m)!的逆元。}
根据费马小定理：
已知gcd(a, p) = 1，则 a^p-1 ≡ 1 (mod p),  所以 a*a^p-2 ≡ 1 (mod p)。
也就是 (m!(n-m)!)的逆元为 (m!(n-m)!)^{p-2 ；}
上面的代码中用到了求幂取模操作来计算(m!(n-m)!)^p-2 % p.下面解释幂取模算法：
反复平方法 求a^b%m
通过研究指数b的二进制表示发现，对任意的整数b都可表示为：


   
   n表示b的实际二进制位数
b_i表示该位是0或1
因此，a^b可表示为：
即用b的每一位表示a的每一项，而对任意相邻的两项存在平方关系，即：
因此我们构造下面的算法：
   
   把b转换为二进制表示，并从右至左扫描其每一位（从低到高）
当扫描到第i位时，根据同余性质（2）计算a的第i项的模：

base变量表示第i-1位时计算出的模，通过递归能很容易地确定所有位的模。
如果第i位为1，即b_i=1,则表示该位需要参与模运算，计算结果 result = (result*base) mod m;其中result为前i-1次的计算结果；若b_i=0，则表示a的第i项为1，不必参与模运算
   
   
int mod(int a,int b,int m){
    int result = 1;
    int base = a;
    while(b>0){
         if(b & 1==1){
            result = (result*base) % m;
         }
         base = (base*base) %m;
         b >>=1;
    }
    return result;
}
其中运用了两个同余性质：
同余性质1：ab≡bc (mod m)
同余性质2：  a≡c (mod m) => a²≡c² (mod m)
理解要点：
   
   base记录了a的每项的模，无论b在该位是0还是1，该结果都记录，目的是给后续位为1的项使用，计算方式是前一结果的平方取模，这也是反复平方法的由来
result只记录了位为1的项的模结果，该计算方式使用了同余性质1
通过地把a使用二进制表示，并结合同余性质1，2，巧妙地化解了大数取模的运算。对1024位这样的大数，也最多进行1024次循环便可计算模值，性能非常快。
该方法是许多西方数学家努力的结果，通常也称为Montgomery算法。


   
   
    
    /这里简单介绍一下Lucas定理
   
   
   
   
    
    对于C(n, m) mod p。这里的n，m，p（p为素数）都很大的情况。就不能再用C(n, m) = C(n - 1，m) + C(n - 1, m - 1)的公式递推了。
    
    这里用到Lusac定理
    
    For non-negative integers m and n and a prime p, the following congruence relation holds:
    
    
     
     
      
      
     
     
    
    
    
    where
    
    
     
     
      
      
     
     
    
    
    
    and
    
    
     
     
      
      
     
     
    
    
    
    are the base p expansions of m and n respectively.
    
     
    
     对于单独的C(ni, mi) mod p，已知C(n, m) mod p = n!/(m!(n - m)!) mod p。显然除法取模，这里要用到m!(n-m)!的逆元。
    
    根据费马小定理：
    
    已知(a, p) = 1，则 a^p-1 ≡ 1 (mod p),  所以 a*a^p-2 ≡ 1 (mod p)。
    
    也就是 (m!(n-m)!)的逆元为 (m!(n-m)!)^{p-2 ；}

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
const int N = 2000100;
typedef long long LL;
LL quick(LL n,LL m,LL p);
LL lucas(LL n,LL m,LL p);
LL power(LL n,LL m,LL p);


int main()
{
    int t;
    scanf("%d", &t);
    while(t--)
    {
        LL n, m, p;
        scanf("%lld %lld %lld",&n, &m, &p);
        printf("%lld\n",lucas(n, m, p));
    }
    return 0;
}


LL lucas(LL n,LL m,LL p)
{
    LL ans=1;
    while(m&&n&&ans)
    {
        ans=((ans%p)*(power(n,m,p)%p))%p;
        m/=p;
        n/=p;
    }
    return ans;
}


LL power(LL n,LL m,LL p)
{
    LL sum1=1, sum2=1;
    for(int i=1;i<=m;i++)
    {
        sum1=((sum1%p)*(n-i+1)%p)%p;
        sum2=((sum2%p)*(i%p))%p;
    }
    LL ans=((sum1%p)*(quick(sum2,p-2,p)))%p;
    return ans;
}


LL quick(LL n,LL m,LL p)
{
    LL r=1;
    while(m!=0)
    {
        if(m&1)
        {
            r=((r%p)*(n%p))%p;
        }
        n=n*n%p;
        m>>=1;
    }
    return r%p;
}