Lucas定理与扩展Lucas

之前看了乘法逆元(详见除法取模与逆元)，发现不能处理不互质的情况，于是去找方法，最后找到了Lucas定理。。。
虽然与期待中的不一样，但是还是非常有用的。

（1）Lucas定理：

若p为素数，则有:

$$C_m^n \equiv \prod_{i=0}^kC_{m_i}^{n_i} \pmod p$$
$$n = n_kp^k+n_{k-1}p^{k-1}+...+n_0$$
$$m = m_kp^k+m_{k-1}p^{k-1}+...+m_0$$
$n_i,m_i$即为把n,m转换成p进制后对应的第i+1位数字。
因为a的p进制的最后一位为a%p，所以原公式可以转化为：
$$C_n^m \equiv C_{[\frac np]}^{[\frac mp]} \times C_{n\mod p}^{m\mod p} \pmod p$$
这个就是我们常使用的公式。
证明：
我们可以用归纳法证明整个定理。我们有下面的式子成立：
$$(1+x)^n \equiv (1+x)^{p[\frac n p]}(1+x)^{n\mod p}\equiv (1+x^p)^{[\frac n p]}(1+x)^{n\mod p}\equiv \{\sum_{i = 0}^{[\frac np]}C_{[\frac np]}^ix^{pi}\}\{\sum_{j=0}^{n\mod p}C_{n\mod p}^jx^j\} \pmod p$$

上式左右两边的x的某项$x^m(m\le n)$的系数对模p同余。
其中左边的$x^m$的系数是 $C_n^m$。 而由于$n\mod p$和$m\mod p$都小于$p$，因此右边的$x^m$一定是由 $x^{[\frac mp]p}$ 和 $x^{m\mod p}$ (即$i=[\frac mp] , j=m\mod p$ ) 相乘而得,因此有：$C_n^m=C_{[\frac np]}^{[\frac mp]} \times C_{n\mod p}^{m\mod p} \pmod p$。

（2）扩展Lucas：

若p不是素数，我们将p分解质因数，将$C_n^m$分别按照（1）中的方法求对p的质因数的模，然后用中国剩余定理合并。
例如：
当我们需要计算$C_n^m\mod p$，其中$p = p_1^{q_1}\times p_2^{q_2}\times ...\times p_k^{q_k}$，我们可以求出：

$$C_n^m\equiv a_i\pmod {p_i^{q_i}} (1\lt i \lt k)$$
然后对于方程组：
$$x\equiv a_i \pmod {p_i^{q_i}}(1\lt i\lt k)$$
我们可以求出满足条件的最小的 $x$,记为$x_0$。
那么我们有：
$$C_n^m\equiv x_0\pmod p$$
但是，我们发现， $p_i^{q_i}$并不是一个素数，它是某个素数的某次方。
下面我们介绍如何计算 $C_n^m \mod p^t(t\ge2,p为素数)$。
计算 $C_n^m\mod p^t$：
我们知道， $C_n^m=\frac {n!}{m!(n-m)!}$，若我们可以计算出 $n!\mod p^t$，我们就能计算出 $m!\mod p^t$以及 $(n-m)!\mod p^t$。我们不妨设 $x=n!\mod p^t,y=m!\mod p^t,z=(n-m)!\mod p^t,$那么答案就是 $x\times reverse(y,p^t)\times reverse(z,p^t)$( $reverse(a,b)$表示计算a对b的乘法逆元)。那么下面问题就转化成如何计算 $n!\mod p^t$。例如 $p=3,t=2,n=19$

$$n!=1\times2\times3\times4\times5\times6\times7\times8\times ...\times19=(1\times2\times4\times5\times7\times8\times...\times16\times17\times19)\times(3\times6\times9\times12\times15\times18)=(1\times2\times4\times5\times7\times8\times...\times16\times17\times19)\times3^6\times(1\times2\times3\times4\times5\times6)$$

后面的部分恰好是$(n/p)!$，于是递归即可。前半部分是以$p^t$为周期的$(1\times2\times4\times5\times7\times8)\equiv(10\times11\times13\times14\times16\times17)\pmod9$。下面是孤立的19，可以知道孤立出来的长度不超过 $p^t$,于是直接计算即可。对于最后剩下的$3^6$这些数我们只要计算出$n!,m!,(n-m)!$里含有多少个$p$(不妨设$x,y,z$)，那么$x-y-z$就是$C_n^m$中$p$的个数，直接计算就行。

下面是计算$C_n^m\mod p$的代码：

#include<cstdio>
#include<algorithm>
using namespace std;

long long pow(long long a, long long b, long long p) {
    long long res = 1;
    while(b) {
        if(b&1) res = res*a%p;
        a = a*a%p;
        b >>= 1;
    }
    return res;
}

long long exgcd(long long a, long long b, long long& x, long long& y) {
    if(!b) {
        x = 1;
        y = 0;
        return a;
    }
    long long res = exgcd(b, a%b, y, x);
    y -= (a/b)*x;
    return res;
}

long long reverse(long long a, long long p) {
    long long x, y;
    exgcd(a, p, x, y);
    return (x % p + p)%p;
}

long long C(long long n, long long m, long long p) {
    if(m>n) return 0;
    long long res = 1, i, a, b;
    for(i = 1; i <= m; i++) {
        a = (n+1-i) % p;
        b = reverse(i%p, p);
        res = res*a%p*b%p;
    }
    return res;
}

long long Lucas(long long n, long long m, long long p) {
    if(m == 0) return 1;
    return Lucas(n/p, m/p, p)*C(n%p, m%p, p) % p;
}

long long cal(long long n, long long a, long long b, long long p) {
    if(!n) return 1;
    long long i, y = n/p, tmp = 1;
    for(i = 1; i <= p; i++) if(i%a) tmp = tmp*i%p;
    long long ans = pow(tmp, y, p);
    for(i = y*p+1; i <= n; i++) if(i%a) ans = ans*i%p;
    return ans * cal(n/a, a, b, p)%p;
}

long long multiLucas(long long n, long long m, long long a, long long b, long long p) {
    long long i, t1, t2, t3, s = 0, tmp;
    for(i = n; i; i/=a) s += i/a;
    for(i = m; i; i/=a) s -= i/a;
    for(i = n-m; i; i/=a) s -= i/a;
    tmp = pow(a, s, p);
    t1 = cal(n, a, b, p);
    t2 = cal(m, a, b, p);
    t3 = cal(n-m, a, b, p);
    return tmp*t1%p*reverse(t2, p)%p*reverse(t3, p)%p;
}


long long exLucas(long long n, long long m, long long p) {
    long long i, d, c, t, x, y, q[100], a[100], e = 0;
    for(i = 2; i*i <= p; i++) {
        if(p % i == 0) {
            q[++e] = 1;
            t = 0;
            while(p%i==0) {
                p /= i;
                q[e] *= i;
                t++;
            }
            if(q[e] == i) a[e] = Lucas(n, m, q[e]);
            else a[e] = multiLucas(n, m, i, t, q[e]);
        }
    }
    if(p > 1) {
        q[++e] = p;
        a[e] = Lucas(n, m, p);
    }
    for(i = 2; i <= e; i++) {
        d = exgcd(q[1], q[i], x, y);
        c = a[i]-a[1];
        if(c%d) exit(-1);
        t = q[i]/d;
        x = (c/d*x%t+t)%t;
        a[1] = q[1]*x+a[1];
        q[1] = q[1]*q[i]/d;
    }
    return a[1];
}

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