卢卡斯定理及拓展

卢卡斯定理

当$m,n$比较大时达到了$longlong$时，需要使用卢卡斯定理

时间复杂度$O(lo{g}_{p}n\ast lo{g}_{2}p)$

取模数p为质数

$A,B$是非负整数，$p$是质数。$A,B$写成$p$进制：$A=a[n]a[n-1]\dots a[0]$，$B=b[n]b[n-1]\dots b[0]$

则组合数$C_A^B$与$(C_{a[n]}^{b[n]}\ast C_{a[n-1]}^{b[n-1]}\dots \ast C_{a[0]}^{b[0]})modp$同余

即：$Lucas(n,m,p)=C(n\%p,m\%p)\ast Lucas(n/p,m/p,p)$

ll qkp(ll x, ll n, ll p) {
    ll ans = 1;
    x %= p;
    while (n) {
        if (n & 1) ans = ans * x % p;
        x = x * x % p;
        n >>= 1;
    }
    return ans;
}

ll inv(ll x, ll p) {  //求逆元
    return qkp(x, p - 2, p);
}

ll cal(ll n, ll m, ll p) {
    if (m > n) return 0;
    ll u = 1, d = 1;
    for (int i = n - m + 1; i <= n; i++) u = u * i % p;
    for (int i = 1; i <= m; i++) d = d * i % p;
    return u * inv(d, p) % p;
}

ll lucas(ll n, ll m, ll p) {
    if (!m) return 1;
    return cal(n % p, m % p, p) * lucas(n / p, m / p, p) % p;
}

如果题目给的$p$较小，可能使算法的时间复杂度退化到$(logn{)}^2$

观察可知卢卡斯定理每次计算的组合数的$n,m$都小于$p$，那么可以预处理保存$p$的范围内的阶乘对$p$取模的结果，此时的时间复杂度为$O(logn)$

int num;
vector<int> prime;
bool is_prime[maxn];
ll fac[maxn][1300];
unordered_map<int, int> mp;

void getPrime() {  //欧拉筛求素数
    int m = sqrt(maxn + 0.5);
    memset(is_prime, 1, sizeof is_prime);
    is_prime[0] = is_prime[1] = 0;
    for (int i = 2; i < maxn; i++) {
        if (is_prime[i]) prime.push_back(i);
        for (int j = 0; j < prime.size() && 1LL * i * prime[j] < maxn; j++) {
            is_prime[i * prime[j]] = 0;
            if (i % prime[j] == 0) break;
        }
    }
    mp.clear(), num = 0;
    for (auto i: prime) mp[i] = ++num;
}

void init() {
    for (auto p: prime) {
        int j = mp[p];
        fac[0][j] = 1;
        for (int i = 1; i < maxn; i++)
            fac[i][j] = fac[i - 1][j] * i % p;
    }
}

ll qkp(ll x, ll n, ll p) {
    ll ans = 1;
    x %= p;
    while (n) {
        if (n & 1) ans = ans * x % p;
        x = x * x % p;
        n >>= 1;
    }
    return ans;
}

ll inv(ll x, ll p) {  //求逆元
    return qkp(x, p - 2, p);
}

ll cal(ll n, ll m, ll p) {
    if (m > n) return 0;
    return fac[n][mp[p]] * inv(fac[m][mp[p]], p) % p * inv(fac[n - m][mp[p]], p) % p;
}

ll lucas(ll n, ll m, ll p) {
    if (!m) return 1;
    return cal(n % p, m % p, p) * lucas(n / p, m / p, p) % p;
}

取模数p为合数

求$C_n^{m}modp$，其中$p$为合数

原理

将$p$质因数分解为$p=p_1^{a_1}\ast p_2^{a_2}...p_k^{a_k}$，问题就变成了求同余方程组：

$\{\begin{array}{r}x\equiv C_n^m(mod{p}_1^{a_1})\\ x\equiv C_n^m(mod{p}_2^{a_2})\\ ......\\ x\equiv C_n^m(mod{p}_n^{a_n})\end{array}$

因为$p_1,p_2,...p_n$互质，最后使用中国剩余定理合并答案

问题就变成了求$C_n^{m}mod{p}^k$

思路

根据组合数的定义，上述问题转化为：$\frac{n!}{m!(n-m)!}mod{p}^k$

因为求逆元的充要条件为$gcd(a,p)=1$，因此不能对$m!,(n-m)!$求逆元，考虑不能求逆元时就是$x!$中含有质因数$p$，考虑除掉这个质因子，则问题变为：

$\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}{p}^{x-y-z}mod{p}^k$，其中$p^x$代表$n!$中质因子$p$的个数

设$f(n)=\frac{n!}{p^x}$，那么问题又转化为$\frac{f(n)}{f(m)f(n-m)}{p}^{x-y-z}mod{p}^k$

只要能求出$f(n)$就大功告成

对$n!$进行变形，将其中$p$的倍数和非倍数分开看，可以得到$n!=(p\ast 2p\ast 3p...tp)(1\ast ...)$

又因为$1-n$中有$\lfloor \frac{n!}{p}\rfloor$个$p$的倍数，那么又能得到$p^{\lfloor \frac{n!}{p}\rfloor }(1\ast 2...\ast t)(1\ast ...)$

即$p^{\lfloor \frac{n!}{p}\rfloor }(\lfloor \frac{n!}{p}\rfloor )!({\prod }_{i=1,i\%p!=0}^{n}i)①$

其中${\prod }_{i=1,i\%p!=0}^{n}i$这个式子是有循环节的，例如$15!mod7$，那么后面的式子为：

$$\begin{gathered} (1\ast 2\ast 3\ast 4\ast 5\ast 6)\ast (8\ast 9\ast 10\ast 11\ast 12\ast 13)\ast 15 \\ =(1\ast 2\ast 3\ast 4\ast 5\ast 6)\ast (8\%7\ast 9\%7\ast 10\%7\ast 10\%7\ast 11\%7\ast 12\%7\ast 13\%7)\ast 15 \\ =(1\ast 2\ast 3\ast 4\ast 5\ast 6{)}^2\ast 15 \end{gathered}$$

因此$①$式又可以化为$p^{\lfloor \frac{n!}{p}\rfloor }(\lfloor \frac{n!}{p}\rfloor )!({\prod }_{i=1,i\%p!=0}^{p^k}i{)}^{\lfloor \frac{n}{p^k}\rfloor }({\prod }_{i=\lfloor \frac{n}{p^k}\rfloor p^k,i\%p!=0}^{n}i)$

因为$p^{\lfloor \frac{n!}{p}\rfloor }$是要提取出来化为$p^{x-y-z}$，但是$(\lfloor \frac{n!}{p}\rfloor )!$还可能包含$p$，那么需要这样定义$f(n)$：

$f(n)=f(\lfloor \frac{n}{p}\rfloor )({\prod }_{i=1,i\%p!=0}^{p^k}i{)}^{\lfloor \frac{n}{p^k}\rfloor }({\prod }_{i=\lfloor \frac{n}{p^k}\rfloor p^k,i\%p!=0}^{n}i)$，递归边界为$f(0)=1$

最后还需要求$p^{x-y-z}$，设$g(n)=\lfloor \frac{n}{p}\rfloor +g(\lfloor \frac{n}{p}\rfloor )$，递归边界为$if(n<p)g(n)=0$

这样分别求出$g(n),g(m),g(n-m)$，则$p^{x-y-z}=p^{g(n)-g(m)-g(n-m)}$

综上，总的公式为：$\frac{f(n)}{f(m)f(n-m)}{p}^{g(n)-g(m)-g(n-m)}mod{p}^k$

inline void exgcd(ll a, ll b, ll &x, ll &y) {
    if (!b) {
        x = 1, y = 0;
        return;
    }
    exgcd(b, a % b, y, x);
    y -= a / b * x;
}

inline ll inv(ll a, ll p) {
    ll x, y;
    exgcd(a, p, x, y);
    return (x % p + p) % p;
}

inline ll mul(ll x, ll n, ll p) {
    ll ans = 0;
    x %= p, n %= p;
    while (n) {
        if (n & 1) ans = (ans + x) % p;
        x = (x + x) % p;
        n >>= 1;
    }
    return ans;
}

inline ll qkp(ll x, ll n, ll p) {
    ll ans = 1;
    x %= p;
    while (n) {
        if (n & 1) ans = ans * x % p;
        x = x * x % p;
        n >>= 1;
    }
    return ans;
}

inline ll f(ll n, ll p, ll mod) {
    if (!n) return 1;
    ll tmp = 1, res = 1;
    for (ll i = 1; i <= mod; i++)
        if (i % p) {
            tmp = tmp * i % mod;
        }
    tmp = qkp(tmp, n / mod, mod);
    for (ll i = mod * (n / mod); i <= n; i++)
        if (i % p) {
            res = res * (i % mod) % mod;
        }
    return f(n / p, p, mod) * tmp % mod * res % mod;
}

inline ll g(ll n, ll p) {
    if (n < p) return 0;
    return g(n / p, p) + n / p;
}

inline ll cal(ll n, ll m, ll p, ll mod) {
    ll fn = f(n, p, mod), fm = inv(f(m, p, mod), mod), fnm = inv(f(n - m, p, mod), mod);
    ll res = qkp(p, g(n, p) - g(m, p) - g(n - m, p), mod);
    return fn * fm % mod * fnm % mod * res % mod;
}

ll A[maxn], B[maxn];

inline ll exlucas(ll n, ll m, ll p) {
    ll res = p, cnt = 0;
    for (ll i = 2; i * i <= p && res > 1; i++)
        if (res % i == 0) {
            ll pk = 1;
            while (res % i == 0) pk *= i, res /= i;
            A[++cnt] = pk, B[cnt] = cal(n, m, i, pk);
        }
    if (res != 1) A[++cnt] = res, B[cnt] = cal(n, m, res, res);
    ll ans = 0;
    for (ll i = 1; i <= cnt; i++) {
        ll M = p / A[i], M1 = inv(M, A[i]);
        ans = (ans + B[i] * M % p * M1 % p) % p;
    }
    return ans;
}