拓展卢卡斯定理学习笔记（附拓展中国剩余定理）

前前言

很久之前学的了，但一直没有机会用到，就写个 $blog$ 防止忘记吧。

题意

求 $C_n^{m}\mod q$ 。
其中， $n,m\leq 10^{18}, q\leq 10^6$ ， $q$ 不一定为质数。
模板题在这里！

前言

$exLucas$ 定理和 $Lucas$ 定理一点关系都没有。。。。
所以根本不需要你会 $Lucas$ 定理，不过你得会中国剩余定理QAQ
好了下面进入正题！！

主要思路

先将 $q$ 分解质因数，变为 $q=p_1^{k_1}p_2^{k_2}...p_m^{k_m}$
分别求出 $C_{n}^{m}\mod p_i^{k_i}$
用中国剩余定理合并

分析

现在的主要问题就是求 $\dfrac{n!}{m!(n-m)!}\mod p_i^{k_i}$ 。
但是 $n,m$ 实在是太太太大了，让这一切非常难顶。
为了让下面存在逆元，我们要把所有的 $p$ 提出来。
令 $n!=f(n)*p^a$ ，其中 $\gcd(f(n),p)=1$ ，也就是把所有 $p$ 的次数都提出来。
于是原式可以写成这种形式： $\dfrac{n!}{m!(n-m)!}\mod p_i^{k_i}\Rightarrow\dfrac{f(n)}{f(m)f(n-m)}*p^{a-b-c} \mod p_i^{k_i}$

现在的问题就变成：

求出 $a,b,c$
求出 $f(n),f(m),f(n-m)\mod p_i^{k_i}$ ，然后求波逆元就完事儿了

第一个问题很好求嘛，令 $g(n)$ 表示 $n!$ 中 $p$ 的次数，有以下递推式： $g(n)=\lfloor\dfrac{n}{p}\rfloor+g(\lfloor\dfrac{n}{p}\rfloor)$

我们重点看第二个问题！！

第二个问题

我们要快速求出 $f(n) \mod p^k$ 。 $n\leq 10^{18},p^k\leq 10^6$ 。 $n!=1\times 2\times 3..\times p...\times n$
我们把 $n!$ 分成两部分，非 $p$ 的倍数和 $p$ 的倍数。

非 $p$ 的倍数

现在我们把所有 $p$ 的倍数拿掉，大概是这个样子： $f(n)=[1\times 2\times 3\times...\times (p-1)]\times[(p+1)\times(p+2)..\times(2p-1)]\times...$
我们把每 $p-1$ 个数放到一个中括号里，称这个为一组。其中，第 $x$ 组为 $[(x-1)p+1,xp-1]$ 。
考虑 $p^k-1$ 这个数，它是 $p^{k-1}$ 组的最后一个数。
那么它的下一组，也就是第 $p^{k-1}+1$ 组，就是 $[p^k+1,p^k+p-1]$ 。
而这一组，在模 $p^k$ 意义下，和第一组一模一样！！
也就是 $i\equiv p^k+i~(mod~p^k),(i\leq p-1)$ 。
这说明，出现了循环节！！
那我们把每 $p^{k-1}$ 组放在一节，乘积记作 $S$ 。
那么总共有 $\lfloor\dfrac{n}{p^k}\rfloor$ 节，这部分乘积为 $S^{\tiny{\lfloor\dfrac{n}{p^k}\rfloor}}$ 。
最后剩的部分的数不超过 $p^k$ 个，我们暴力求即可。
于是，我们就在 $O(p^k)$ 时间内求出了非 $p$ 的倍数的乘积，记作 $A$ 。

$p$ 的倍数

考虑 $p\times 2p\times 3p...\times \lfloor\dfrac{n}{p}\rfloor p$ 。
我们去掉所有 $p$ ，就变成 $(\lfloor\dfrac{n}{p}\rfloor)!$ 。
接下来我们就是要求 $f(\lfloor\dfrac{n}{p}\rfloor)$ 。这部分递归求就可以了。

综上

$f(n)=\begin{cases}1 &(n=1)\cr A \times f(n/p)&(otherwise)\end{cases}$
复杂度是 $O(p^klogn)$ 。

最后用中国剩余定理合并即可。

拓展中国剩余定理

考虑到可能有人不会 $CRT$ ，而且我怕自己忘了，就在这里把 $exCRT$ 介绍一下吧。
下面要求以下线性同余方程组的一个解： $\begin{cases}x\equiv a_1~(mod~b_1)\cr x\equiv a_2~(mod~b_2)\cr{...}\cr x\equiv a_m(mod~b_m)\end{cases}$

其中， $\{b_i\}$ 不一定两两互质。

考虑我们已经求出了前 $k-1$ 个方程的一个解 $X$ ，现在求前 $k$ 个方程的一个解。
令 $B=lcm(b_1,b_2,...,b_{k-1})$
那么前 $k$ 个方程的所有解为 $X+tB(t\in Z)$
考虑某个解满足第 $k$ 个方程，即： $X+tB\equiv a_{k}~(mod~b_k)$
我们的目标是求出某个满足的 $t$ ，移一下项，有： $tB\equiv a_k-X~(mod~b_k)$ 。
令 $g=\gcd(B,b_k)$ ，如果 $g\nmid a$ ，那么方程组无解。
否则，方程等价为： $tB'\equiv a'~(mod~b')$
其中 $B'=B/g,a'=(a_k-X)/g,b'=b_k/g$ 。
于是有 $\gcd(B',b')=1$ ， $B'$ 关于 $b'$ 的逆元存在，于是可以求出某个解 $t$ ，即： $t\equiv a'\times B'^{-1}~(mod~b')$ 。
于是我们可以在 $O(nlogV)$ 时间内求出线性同余方程组的解。

代码超级无敌清晰！

总代码如下

#include <bits/stdc++.h>
#define m_p make_pair
#define N 100005
using namespace std;
typedef long long LL;
typedef unsigned long long uLL;

LL z = 1;
LL read(){
    LL x, f = 1;
    char ch;
    while(ch = getchar(), ch < '0' || ch > '9') if(ch == '-') f = -1;
    x = ch - '0';
    while(ch = getchar(), ch >= '0' && ch <= '9') x = x * 10 + ch - 48;
    return x * f;
}

LL ksm(LL a, LL b, LL p){
    LL s = 1;
    while(b){
        if(b & 1) s = z * s * a % p;
        a = z * a * a % p;
        b >>= 1;
    }
    return s;
}

vector<pair<LL, LL> > d;

void get(LL x){//分解质因数，得到 p 和 p^k
    for(LL i = 2; i <= x / i; i++){
        if(x % i == 0){
            LL tot = 1;
            while(x % i == 0) tot *= i, x /= i;
            d.push_back(m_p(i, tot));
        }
    }
    if(x > 1) d.push_back(m_p(x, x));
}

LL f(LL n, LL p, LL pk){//求 f(n) 
    if(n == 0 || n == 1) return 1;
    LL i, s = 1;
    for(i = 1; i <= pk; i++) if(i % p) s = s * i % pk;
    s = ksm(s, n / pk, pk);
    for(i = n / pk * pk; i <= n; i++) if(i % p) s = s * (i % pk) % pk;
    return s * f(n / p, p, pk) % pk;
}

LL g(LL n, LL p){//求 g(n) 
    if(n < p) return 0;
    return n / p + g(n / p, p);
}

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

LL inv(LL a, LL b){
    LL x, y;
    exgcd(a, b, x, y);
    return (x + b) % b;//求逆元 
}

LL excrt(LL n, LL *a, LL *b){//excrt 合并解 
    LL x, B, M, t, k, g;
    __int128 z = 1;
    x = a[1], B = b[1];
    for(LL i = 2; i <= n; i++){
        g = __gcd(B, b[i]);
        M = B / g * b[i];
        if((a[i] - x) % g != 0) return -1;
        exgcd(B / g, b[i] / g, t, k);
        t = (z * t * (a[i] - x) / g % M + M) % M;
        x = (x + z * t * B % M) % M;
        B = M;
    }
    return x;
}

LL a[N], b[N];

LL exLucas(LL n, LL m, LL p){
    LL ans = 1, i, j, k, pk, w, cnt = 0;
    get(p);
    for(auto t: d){
        p = t.first, pk = t.second;
        i = f(n, p, pk);
        j = inv(f(m, p, pk), pk); k = inv(f(n - m, p, pk), pk);
        w = g(n, p) - g(m, p) - g(n - m, p);
        i = i * j % pk * k % pk;
        i = i * ksm(p, w, pk) % pk;
        a[++cnt] = i, b[cnt] = pk;//得到 C(n, m) % p^k 
    }
    return excrt(cnt, a, b);
}

int main(){
    int i, j;
    LL n, m, p;
    n = read(); m = read(); p = read();
    printf("%lld", exLucas(n, m, p));
    return 0;
}