先来一篇大佬的博客 Fading 的博客
然后如果你能认真研读上面大佬的博客的话,应该就比较轻松了,我比较蒟蒻,还是不要多说的好=w=
Code:
1 #include <bits/stdc++.h> 2 #define ll long long 3 using namespace std; 4 const int N = 1007; 5 ll n, m, p; 6 ll a[N], b[N]; 7 void exgcd(ll a, ll b, ll &x, ll &y) {//用于求逆元 8 if (b == 0) { 9 x = 1; 10 y = 0; 11 return; 12 } 13 exgcd(b, a % b, x, y); 14 ll t = x; 15 x = y; 16 y = t - a / b * y; 17 } 18 ll inv(ll a, ll p) {//求a在p意义下的逆元 19 ll x, y; 20 exgcd(a, p, x, y); 21 return (x + p) % p;//保证其为非负 22 } 23 ll mul(ll a, ll b, ll p) {//龟速乘防爆 24 ll re = 0; 25 a %= p; 26 b %= p; 27 for (; b; b >>= 1) { 28 if (b & 1) { 29 re = (re + a) % p; 30 } 31 a = (a + a) % p; 32 } 33 return re; 34 } 35 ll pw(ll x, ll y, ll p) {//快速幂 36 ll re = 1; 37 x %= p; 38 for (; y; y >>= 1) { 39 if (y & 1) { 40 re = re * x % p; 41 } 42 x = x * x % p; 43 } 44 return re; 45 } 46 ll fac(ll n, ll p, ll pk) {//求n! % p**k 47 if (n == 0) { 48 return 1; 49 } 50 ll rou = 1;//求循环节 51 ll rem = 1;//求余项 52 for (ll i = 2; i <= pk; i++) { 53 if (i % p) { 54 rou = rou * i % pk; 55 } 56 } 57 rou = pw(rou, n / pk, pk); 58 for (ll i = pk * (n / pk); i <= n; i++) { 59 if (i % p) { 60 rem = rem * (i % pk) % pk; 61 } 62 } 63 return fac(n / p, p, pk) * rou % pk * rem % pk; 64 } 65 ll G(ll n, ll p) {//算一个次幂 66 if (n < p) { 67 return 0; 68 } 69 return G(n / p, p) + (n / p); 70 } 71 ll C_pk(ll n, ll m, ll p, ll pk) {//求组合数C(n,m) % p**k 72 ll f = fac(n, p, pk); 73 ll inv1 = inv(fac(m, p, pk), pk); 74 ll inv2 = inv(fac(n - m, p, pk), pk); 75 ll mi = pw(p, G(n, p) - G(m, p) - G(n - m, p), pk); 76 return f * inv1 % pk * inv2 % pk * mi % pk; 77 } 78 ll exlucas(ll n, ll m, ll p) {//拓展卢卡斯定理 79 ll re = p, tot = 0; 80 for (ll i = 2; i * i <= p; i++) { 81 if (re % i == 0) { 82 ll pk = 1; 83 while (re % i == 0) { 84 pk *= i; 85 re /= i; 86 } 87 a[++tot] = C_pk(n, m, i, pk);//先将p分解,然后分别记下方程组 88 b[tot] = pk; 89 } 90 } 91 if (re != 1) { 92 a[++tot] = C_pk(n, m, re, re); 93 b[tot] = re; 94 } 95 ll ans = 0; 96 for (ll i = 1; i <= tot; i++) {//再用中国剩余定理合并求出答案 97 ll t = p / b[i]; 98 ll inv1 = inv(t, b[i]); 99 ans = (ans + a[i] * t % p * inv1 % p) % p; 100 } 101 return ans; 102 } 103 int main () { 104 scanf("%lld%lld%lld", &n, &m, &p); 105 printf("%lld\n", exlucas(n, m, p)); 106 return 0; 107 }