扩展卢卡斯定理

  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 }