之前看了乘法逆元(详见除法取模与逆元),发现不能处理不互质的情况,于是去找方法,最后找到了Lucas定理。。。
虽然与期待中的不一样,但是还是非常有用的。
(1)Lucas定理:
若p为素数,则有:
ni,mi 即为把n,m转换成p进制后对应的第i+1位数字。
因为a的p进制的最后一位为a%p,所以原公式可以转化为:
这个就是我们常使用的公式。
证明:
我们可以用归纳法证明整个定理。我们有下面的式子成立:
上式左右两边的x的某项
xm(m≤n)
的系数对模p同余。
其中左边的
xm
的系数是
Cmn
。 而由于
nmodp
和
mmodp
都小于
p
,因此右边的
(2)扩展Lucas:
若p不是素数,我们将p分解质因数,将
Cmn
分别按照(1)中的方法求对p的质因数的模,然后用中国剩余定理合并。
例如:
当我们需要计算
Cmnmodp
,其中
p=pq11×pq22×...×pqkk
,我们可以求出:
然后对于方程组:
我们可以求出满足条件的最小的 x ,记为
那么我们有:
但是,我们发现, 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/p)!
,于是递归即可。前半部分是以
pt
为周期的
(1×2×4×5×7×8)≡(10×11×13×14×16×17)(mod9)
。下面是孤立的19,可以知道孤立出来的长度不超过
pt
,于是直接计算即可。对于最后剩下的
36
这些数我们只要计算出
n!,m!,(n−m)!
里含有多少个
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;
}