Lucas定理与扩展Lucas

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

（1）Lucas定理：

若p为素数，则有:

Cnm≡∏i=0kCnimi(modp)

C_m^n \equiv \prod_{i=0}^kC_{m_i}^{n_i} \pmod p

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

n = n_kp^k+n_{k-1}p^{k-1}+...+n_0

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

m = m_kp^k+m_{k-1}p^{k-1}+...+m_0
ni,mi n_i,m_i即为把n,m转换成p进制后对应的第i+1位数字。
因为a的p进制的最后一位为a%p，所以原公式可以转化为：

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

C_n^m \equiv C_{[\frac np]}^{[\frac mp]} \times C_{n\mod p}^{m\mod p} \pmod p
这个就是我们常使用的公式。
证明：
我们可以用归纳法证明整个定理。我们有下面的式子成立：

(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)

(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的某项 xm(m≤n) x^m(m\le n)的系数对模p同余。
其中左边的 xm x^m的系数是 Cmn C_n^m。而由于 nmodp n\mod p和 mmodp m\mod p都小于 p p，因此右边的xmx^m一定是由 x[mp]p x^{[\frac mp]p} 和 xmmodp x^{m\mod p} (即 i=[mp],j=mmodp i=[\frac mp] , j=m\mod p ) 相乘而得,因此有： Cmn=C[mp][np]×Cmmodpnmodp(modp) C_n^m=C_{[\frac np]}^{[\frac mp]} \times C_{n\mod p}^{m\mod p} \pmod p。

（2）扩展Lucas：

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

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

C_n^m\equiv a_i\pmod {p_i^{q_i}} (1\lt i \lt k)
然后对于方程组：

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

x\equiv a_i \pmod {p_i^{q_i}}(1\lt i\lt k)
我们可以求出满足条件的最小的 x x,记为x0x_0。
那么我们有：

Cmn≡x0(modp)

C_n^m\equiv x_0\pmod p
但是，我们发现， pqii p_i^{q_i}并不是一个素数，它是某个素数的某次方。
下面我们介绍如何计算 Cmnmodpt(t≥2,p为素数) C_n^m \mod p^t(t\ge2,p为素数)。
计算 Cmnmodpt C_n^m\mod p^t：
我们知道， Cmn=n!m!(n−m)! C_n^m=\frac {n!}{m!(n-m)!}，若我们可以计算出 n!modpt n!\mod p^t，我们就能计算出 m!modpt m!\mod p^t以及 (n−m)!modpt (n-m)!\mod p^t。我们不妨设 x=n!modpt,y=m!modpt,z=(n−m)!modpt, x=n!\mod p^t,y=m!\mod p^t,z=(n-m)!\mod p^t,那么答案就是 x×reverse(y,pt)×reverse(z,pt) x\times reverse(y,p^t)\times reverse(z,p^t)( reverse(a,b) reverse(a,b)表示计算a对b的乘法逆元)。那么下面问题就转化成如何计算 n!modpt n!\mod p^t。例如 p=3,t=2,n=19 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!=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)! (n/p)!，于是递归即可。前半部分是以 pt p^t为周期的 (1×2×4×5×7×8)≡(10×11×13×14×16×17)(mod9) (1\times2\times4\times5\times7\times8)\equiv(10\times11\times13\times14\times16\times17)\pmod9。下面是孤立的19，可以知道孤立出来的长度不超过 pt p^t,于是直接计算即可。对于最后剩下的 36 3^6这些数我们只要计算出 n!,m!,(n−m)! n!,m!,(n-m)!里含有多少个 p p(不妨设x,y,zx,y,z)，那么 x−y−z x-y-z就是 Cmn C_n^m中 p p的个数，直接计算就行。

下面是计算CmnmodpC_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;
}