Lucas定理与扩展Lucas

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

（1）Lucas定理：

若p为素数，则有:

Cnm≡∏i=0kCnimi(modp)

n=nkpk+nk−1pk−1+...+n0

m=mkpk+mk−1pk−1+...+m0

ni,mi即为把n,m转换成p进制后对应的第i+1位数字。
因为a的p进制的最后一位为a%p，所以原公式可以转化为：

Cmn≡C[mp][np]×Cmmodpnmodp(modp)

这个就是我们常使用的公式。
证明：
我们可以用归纳法证明整个定理。我们有下面的式子成立：

(1+x)n≡(1+x)p[np](1+x)nmodp≡(1+xp)[np](1+x)nmodp≡{∑i=0[np]Ci[np]xpi}{∑j=0nmodpCjnmodpxj}(modp)

上式左右两边的x的某项xm(m≤n)的系数对模p同余。
其中左边的xm的系数是 Cmn。 而由于nmodp和mmodp都小于p，因此右边的xm一定是由 x[mp]p 和 xmmodp (即i=[mp],j=mmodp ) 相乘而得,因此有：Cmn=C[mp][np]×Cmmodpnmodp(modp)。

（2）扩展Lucas：

若p不是素数，我们将p分解质因数，将Cmn分别按照（1）中的方法求对p的质因数的模，然后用中国剩余定理合并。
例如：
当我们需要计算Cmnmodp，其中p=pq11×pq22×...×pqkk，我们可以求出：

Cmn≡ai(modpqii)(1<i<k)

然后对于方程组：

x≡ai(modpqii)(1<i<k)

我们可以求出满足条件的最小的x,记为x0。
那么我们有：

Cmn≡x0(modp)

但是，我们发现，pqii并不是一个素数，它是某个素数的某次方。
下面我们介绍如何计算Cmnmodpt(t≥2,p为素数)。
计算Cmnmodpt：
我们知道，Cmn=n!m!(n−m)!，若我们可以计算出n!modpt，我们就能计算出m!modpt以及(n−m)!modpt。我们不妨设x=n!modpt,y=m!modpt,z=(n−m)!modpt,那么答案就是x×reverse(y,pt)×reverse(z,pt)(reverse(a,b)表示计算a对b的乘法逆元)。那么下面问题就转化成如何计算n!modpt。例如p=3,t=2,n=19

n!=1×2×3×4×5×6×7×8×...×19=(1×2×4×5×7×8×...×16×17×19)×(3×6×9×12×15×18)=(1×2×4×5×7×8×...×16×17×19)×36×(1×2×3×4×5×6)

后面的部分恰好是(n/p)!，于是递归即可。前半部分是以pt为周期的(1×2×4×5×7×8)≡(10×11×13×14×16×17)(mod9)。下面是孤立的19，可以知道孤立出来的长度不超过 pt,于是直接计算即可。对于最后剩下的36这些数我们只要计算出n!,m!,(n−m)!里含有多少个p(不妨设x,y,z)，那么x−y−z就是Cmn中p的个数，直接计算就行。

下面是计算Cmnmodp的代码：

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