【数论】组合数求模

Sprjdfn在努力的补习数学。
组合数求模是OI中常考的一个知识点，但是，我不会。。。
所以只好学习一下了，首先分享关于组合数求模的几个资源：
http://www.cnblogs.com/ykzou/p/4494902.html
http://yunpan.cn/cjUuySy5E9ZKB （提取码：f098）

然而如何求组合数对一个合数求模，我也没弄明白，但好像这篇文章里面个有所介绍，有兴趣的同学可以看一看。
http://www.mpqweb.com/blog/?p=123

首先我们先想一想暴力的方法：
1.C(n, m)=n!/(m!*(n-m)!)
暴力For一遍，但我们很容易发现，在n,m都比较大的时候，long long 都存不下了。所以这种方法只适用于较小的数据
2.将分子的每个数存在一个数组里，然后拿分母的数去约分，最后在将数组里的数乘起来就可以了。
这里附上一道题目：BJOI 2015 Day 1 T3
题目可以看dysnpv的博客：
http://blog.csdn.net/dysnpv/article/details/45464605
然后附上我的AC代码：

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

using namespace std;

typedef long long ll;

bool npr[100050];
ll n, m, k, p;
ll note[100050];
ll t[100050];

int main()
{
    scanf("%lld %lld %lld %lld", &n, &m, &k, &p);
    for (int i = 2; i <= m; i++)
        if (!npr[i]) {
            t[i]++;
            for (int j = i+i; j <= m; j += i) {
                npr[j] = true;
                int x = j;
                while (x % i == 0) {
                    x /= i;
                    t[i]++;
                }
            }
        }
    for (ll i = k; i < m+k; i++)
        note[i-k+1] = i;
    for (int i = 2; i <= m; i++)
        if (t[i]) {
            for (int j = ((k-1)/i*i+i)-k+1; j <= m; j += i) {
                while (note[j] % i == 0 && t[i]) {
                    note[j] /= i;
                    t[i]--;
                }
                if (!t[i])
                    break;
            }
        }
    ll c = 1;
    for (int i = 1; i <= m; i++) {
        c *= note[i];
        c %= p;
    }
    ll ans = 1;
    for (int i = c; i > c-n; i--) {
        ans *= i;
        ans %= p;
    }
    printf("%lld", ans);
    return 0;
}

说真的，我觉得出题人估计都没想到可以这么做，这种思路是我的同学提供的。

接下来说正经的做法：
首先先说一下Lucas定理：
对于任意质数p有
C(n, m) % p = Lucas(n, m, p)
Lucas(n, m, p) = C(n%p, m%p) * Lucas(n/p, m/p, p)
也就是:

然后借用jasaiq博客里的一句话：

我们已经知道了，Lucas定理的公式，但是关键是怎么编码来得出结果呢？
也许你会说，那还不简单，直接递归就搞定了，其实没那么简单。

毕竟C(n, m)在比较大的数据规模中也不是好求的东西。
接下来我们来看如何用乘法逆元求组合数对一个质数p取模的结果：

C(n, m) = n! / ( m! * (n-m)! )
x = m! * (n-m)! 的逆元
∴C(n, m) = n! / (m! * (n-m)! ) ≡ n! * x (mod p)
又 x * ( m! * (n-m)! ) ≡ 1 (mod p)
令 ( m! * (n-m)! ) = y
∴x * y + p * q = 1 = gcd(x, p)
又 y已知, p已知
∴可以用拓展欧几里得算法求出x

由此C(n, m) % p (p为质数)就可以求出了。
这种方法的适用范围大概在 n, m≤10^5 之内
代码如下：

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <algorithm>

using namespace std;

typedef long long ll;

ll fac[100050];

void getfac(ll p)
{
    fac[0] = 1;
    for (ll i = 1; i <= p; i++)
        fac[i] = fac[i-1]*i % p;
}

void exp_gcd(ll a, ll b, ll &x, ll &y)
{
    if (b == 0) {
        x = 1, y = 0;
        return;
    }
    exp_gcd(b, a%b, y, x);
    y -= x * (a/b);
}

ll inv(ll a, ll mod)
{
    a %= mod;
    ll x, y;
    exp_gcd(a, mod, x, y);
    while (x < 0)
        x += mod;
    return x % mod;
}

ll c(ll n, ll m, ll p)
{
    if (m > n)
        return 0;
    return fac[n] * inv(fac[m]*fac[n-m], p) % p;
}

ll Lucas(ll n, ll m, ll p)
{
    if (!m)
        return 1;
    return c(n%p, m%p, p) * Lucas(n/p, m/p, p) % p;
}

int main()
{
    ll n, m, p;
    scanf("%lld %lld %lld", &n, &m, &p);
    getfac(p);
    printf("%lld\n", Lucas(n, m, p));
    return 0;
}