基础数论复习

时间 2019-11-09

标签基础数论复习繁體版

原文原文链接

目录php

这个复习没有什么顺序的... 想到什么复习什么而已qwqhtml

Miller-Rabin 质数测试

问题描述

测试一个正整数 $p$ 是否为素数，须要较快且准确的测试方法。$(p \le 10^{18})$c++

算法解决

费马小定理

首先须要了解 费马小定理 。算法

费马小定理：对于质数 $p$ 和任意整数 $a$ ，有 $a^p \equiv a \pmod p$ 。函数

反之，若知足 $a^p \equiv a \pmod p$，$p$ 也有很大几率为质数。学习

将两边同时约去一个 $a$ ，则有 $a^{p-1} \equiv 1 \pmod p$ 。测试

也便是说：假设咱们要测试 $n$ 是否为质数。咱们能够随机选取一个数 $a$ ，而后计算 $a^{n-1} \pmod n$ ，若是结果不为 $1$ ，咱们能够 $100\%$ 判定 $n$ 不是质数。不然咱们再随机选取一个新的数 $a$ 进行测试。如此反复屡次，若是每次结果都是 $1$ ，咱们就假定 $n$ 是质数。优化

二次探测定理

该测试被称为 $Fermat$ 测试。$Miller$ 和 $Rabin$ 在 $Fermat$ 测试上，创建了 $Miller-Rabin$ 质数测试算法。与 $Fermat$ 测试相比，增长了一个 二次探测定理 。ui

二次探测定理：若是 $p$ 是奇素数，则 $x^2 \equiv 1 \pmod p$ 的解为 $x \equiv 1$ 或 $x \equiv p - 1 \pmod p$ 。spa

这是很容易证实的：
\[ \begin{align} x^2 \equiv 1 \pmod p \\ x^2 -1 \equiv 0 \pmod p \\ (x-1) (x+1) \equiv 0 \pmod p \end{align} \]

又 $\because$ $p$ 为奇素数，有且仅有 $1,p$ 两个因子。

$\therefore$ 只有两解 $x \equiv 1$ 或 $x \equiv p - 1 \pmod p$ 。

若是 $a^{n-1} \equiv 1 \pmod n$ 成立，$Miller-Rabin$ 算法不是当即找另外一个 $a$ 进行测试，而是看 $n-1$ 是否是偶数。若是 $n-1$ 是偶数，另 $\displaystyle u=\frac{n-1}{2}$ ，并检查是否知足二次探测定理即 $a^u \equiv 1$ 或 $a^u \equiv n - 1 \pmod n$ 。

假设 $n=341$ ，咱们选取的 $a=2$ 。则第一次测试时，$2^{340} \bmod 341 = 1$ 。因为 $340$ 是偶数，所以咱们检查 $2^{170}$ ，获得 $2^{170} \bmod 341=1$ ，知足二次探测定理。同时因为 $170$ 仍是偶数，所以咱们进一步检查$2^{85} \bmod 341=32$ 。此时不知足二次探测定理，所以能够断定 $341$ 不为质数。

将这两条定理合起来，也就是最多见的 $Miller-Rabin$ 测试 。

单次测试失败的概率为 $\displaystyle \frac{1}{4}$ 。那么假设咱们作了 $times$ 次测试，那么咱们总错误的概率为 $\displaystyle (\frac{1}{4})^{times}$ 若是 $times$ 为 $20$ 次，那么错误几率约为 $0.00000000000090949470$ 也就是意味着不可能发生的事件。

单次时间复杂度是 $O(\log n)$ ，总复杂度就是 $O(times \log n)$ 了。注意 long long 会爆，用 __int128 就好了。

代码实现

inline bool Miller_Rabin(ll p) {
    if (p <= 2) return p == 2;
    if (!(p & 1)) return false;
    ll u = p - 1; int power = 0;
    for (; !(u & 1); u >>= 1) ++ power;
    For (i, 1, times) {
        ll a = randint(2, p - 1), x = fpm(a, u, p), y;
        for (int i = 1; i <= power; ++ i, x = y) {
            if ((y = mul(x, x, p)) == 1 && x != 1 && x != p - 1) return false;
        }
        if (x != 1) return false;
    }
    return true;
}

Pollard-Rho 大合数分解

问题描述

分解一个合数 $n$ ，时间效率要求较高。（比试除法 $O(\sqrt n)$ 要快许多）

算法解决

生日悖论

在一个班级里，假设有 $60$ 人，全部人生日不一样几率是多少？

依次按人考虑，第一我的有 $1$ 的几率，第二我的要与第一我的不一样就是 $\displaystyle 1 - \frac{1}{365}$ ，第三我的与前两人不一样，那么就是 $\displaystyle 1 - \frac{2}{365}$ 。那么第 $i$ 我的就是 $\displaystyle 1 - \frac{i-1}{365}$ 。

那么几率就是
\[ \displaystyle \prod_{i=1}^{n} (1 - \frac{i-1}{365}) = \prod _{i=0}^{n-1}\frac{365-i}{365}=\frac{365^{\underline{n}}}{365^n} \]

假设咱们带入 $n$ 等于 $60$ ，那么这个几率就是 $0.0058$ ，也就是说基本上不可能发生。

好像是由于和大众的常识有些违背，因此称做生日悖论。

假设一年有 $N$ 天，那么只要 $n \ge \sqrt N$ 存在相同生日的几率就已经大于 $50 \%$ 了。

利用伪随机数判断

本段参考了 Doggu 大佬的博客。

咱们考虑利用以前那个性质来判断，生成了许多随机数，而后两两枚举判断。

假设枚举的两个数为 $i, j$ ，假设 $i > j$ ，那么判断一下 $\gcd(i - j, n)$ 是多少，若是非 $1$ ，那么咱们就已经找出了一个 $n$ 的因子了，这样的几率随着枚举的组数变多会愈来愈快。

若是咱们一开始把这些用来判断的数都造出来这样，这样就很坑了。不只内存有点恶心，并且其实效率也不够高。

考虑优化，咱们用一个伪随机数去判断，不断生成两个随机数，而后像流水线同样逐个判断就好了。

有一个随机函数彷佛很优秀，大部分人都用的这个：
\[ f(x)=(x^2+a) \bmod n \]

$a$ 在单次判断中是一个常数。

而后通常这种简单的随机数都会出现循环节，咱们判断一下循环节就好了。

但为了节省内存，须要接下来这种算法。

Floyd 判圈法

两我的同时在一个环形跑道上起跑，一我的是另一我的的两倍速度，那么这我的绝对会追上另一我的。追上时，意味着其中一我的已经跑了一圈了。

而后就能够用这个算法把退出条件变松，只要有环，咱们就直接退出就好了。

若是此次失败了，那么继续找另一个 $a$ 就好了。

用 Miller-Rabin 优化

这个算法对于因子多，因子值小的数 $n$ 是较为优异的。

也就是对于它的一个小因子 $p$ ， $Pollard-Rho$ 指望在 $O(\sqrt p)$ 时间内找出来。

但对于 $n$ 是一个质数的状况，就没有什么较好的方法快速判断了，那么复杂度就是 $O(\sqrt n)$ 的。

此时就退化成暴力试除法了。

此时咱们考虑用前面学习的 $Miller-Rabin$ 算法来优化这个过程就好了。

此时把这两个算法结合在一块儿，就能解决这道题啦qwq

代码实现

我彷佛要特判因子为 $2$ 的数，那个彷佛一直在里面死循环qwq

namespace Pollard_Rho {
    ll c, Mod;

    ll f(ll x) { return (mul(x, x, Mod) + c) % Mod; }

    ll find(ll x) {
        if (!(x & 1)) return 2; Mod = x;
        ll a = randint(2, x - 1), b = a; c = randint(2, x - 1);
        do {
            a = f(a); b = f(f(b));
            ll p = __gcd(abs(a - b), x); 
            if (p > 1) return p;
        } while (b != a);
        return find(x);
    }

    void ReSolve(ll x, vector<ll> &factor) {
        if (x <= 1) return;
        if (Miller_Rabin(x)) { factor.pb(x); return; }
        ll fac = find(x); ReSolve(fac, factor); ReSolve(x / fac, factor);
    }
};

约瑟夫问题

问题描述

首先 $n$ 个候选人围成一个圈，依次编号为 $0\dots n-1$ 。而后随机抽选一个数 $k$ ，并 $0$ 号候选人开始按从 $1$ 到 $k$ 的顺序依次报数，$n-1$号候选人报数以后，又再次从 $0$ 开始。当有人报到 $k$ 时，这我的被淘汰，从圈里出去。下一我的从 $1$ 开始从新报数。

如今询问：最后一个被淘汰在哪一个位置。 $(n \le 10^9, k \le 1000)$ 。

算法解决

暴力模拟

最直观的解法是用循环链表模拟报数、淘汰的过程，复杂度是 $O(nk)$ 。

递推一

令 $f[n]$ 表示当有 $n$ 个候选人时，最后当选者的编号。 （注意是有 $n$ 我的，而不是 $n$ 轮）
\[ \begin{cases} f[1] = 0\\ f[n] = (f[n - 1] + k) \pmod n &n > 1 \end{cases} \]

咱们考虑用数学概括法证实：

\[f[1]=0\]

显然当只有 $1$ 个候选人时，该候选人就是当选者，而且他的编号为 $0$ 。
\[f[n] = (f[n - 1] + k) \pmod n\]

假设咱们已经求解出了 $f[n - 1]$ ，而且保证 $f[n - 1]$ 的值是正确的。

如今先将 $n$ 我的按照编号进行排序：
0 1 2 3 ... n-1
那么第一次被淘汰的人编号必定是 $k-1$ (假设 $k < n$ ，若 $k > n$ 则为 $(k-1) \bmod n$ )。将被选中的人标记为 $\underline{\#}$ ：
0 1 2 3 ... k-2 # k k+1 k+2 ... n-1
第二轮报数时，起点为 $k$ 这个候选人。而且只剩下 $n-1$ 个选手。假如此时把 $k$ 看做 $0'$ ，$k+1$ 看做 $1'$ ... 则对应有：
0  1 2 3 ... k-2  # k  k+1 k+2 ... n-1
n-k'   ...     n-2'   0'  1'  2' ... n-k-1'
此时在 $0',1',...,n-2'$ 上再进行一次 $k$ 报数的选择。而 $f[n-1]$ 的值已经求得，所以咱们能够直接求得当选者的编号 $s'$ 。

可是，该编号 $s'$ 是在 $n-1$ 个候选人报数时的编号，并不等于 $n$ 我的时的编号，因此咱们还须要将$s'$ 转换为对应的 $s$ 。经过观察，$s$ 和 $s'$ 编号相对偏移了 $k$ ，又由于是在环中，所以获得 $s = (s'+k) \bmod n$ ，即 $f[n] = (f[n - 1] + k) \pmod n$ 。

而后就证实完毕了。

这个时间复杂度就是 $O(n)$ 的了，算比较优秀了，但对于这题来讲仍是不行。

递推二

所以当 $n$ 小于 $k$ 时，就只能采用第一种递推的算法来计算了。

那么咱们就能够用第二种递推，解决的思路仍然和上面相同，而区别在于咱们每次减小的 $n$ 的规模再也不是 $1$ 。

一样用一个例子来讲明，初始 $n=10,k=4$ :

初始序列：
0 1 2 3 4 5 6 7 8 9
当 $7$ 号进行过报数以后：
0 1 2 - 4 5 6 - 8 9
而对于任意一个 $n,k$ 来讲，退出的候选人数量为 $\displaystyle \lfloor \frac{n}{k} \rfloor$ 。

因为此时起点为 $8$ ，则等价于：
2 3 4 - 5 6 7 - 0 1
所以咱们仍然能够从 $f[8]$ 的结果来推导出 $f[10]$ 的结果。

咱们须要根据 $f[8]$ 的值进行分类讨论。假设 $f[8]=s$ ，则根据 $s$ 和 $n \bmod k$ 的大小关系有两种状况：
\[ \begin{cases} s' = s - n \bmod k + n & (s< n \bmod k) \\ \displaystyle s' = s - n \bmod k + \lfloor \frac{(s - n \bmod k)}{k-1} \rfloor & (s \ge n \bmod k) \end{cases} \]

上面那种就是最后剩的的几个数，对于上面例子就是 $s=\{0, 1\}$ ，$s' = \{8, 9\}$ 。

下面那种就是其他的几个数，对于上面例子就是 $s = \{5,6,7\}$ ，$s' = \{4,5,6\}$ 。

这两种就是先定位好初始位置，而后再加上前面的空隙，获得新的位置。

因为咱们不断的在减少 $n$ 的规模，最后必定会将 $n$ 减小到小于 $k$ ，此时 $\displaystyle \lfloor \frac{n}{k} \rfloor = 0$ 。

所以当 $n$ 小于 $k$ 时，就只能采用第一种递推的算法来计算了。

每次咱们将 $\displaystyle n \to n - \lfloor \frac{n}{k} \rfloor$ ，咱们能够近似看作把 $n$ 乘上 $\displaystyle 1-\frac{1}{k}$ 。

按这样分析的话，复杂度就是 $\displaystyle O(-\log_{\frac{k-1}{k}}n+k)$ 的。这个对于 $k$ 很小的时候会特别快。

代码实现

int Josephus(int n, int k) {
    if (n == 1) return 0;
    int res = 0;
    if (n < k) {
        For (i, 2, n)
            res = (res + k) % i;
        return res;
    }
    res = Josephus(n - n / k, k);
    if (res < n % k) 
        res = res - n % k + n;
    else 
        res = res - n % k + (res - n % k) / (k - 1);
    return res;
}

扩展欧几里得

问题描述

对于三个天然数 $a,b,c$ ，求解 $ax+by = c$ 的 $(x, y)$ 的整数解。

算法解决

首先咱们要判断是否存在解，对于这个这个存在整数解的充分条件是 $\gcd(a, b)~ |~c$ 。

也就是说 $c$ 为 $a,b$ 最大公因数的一个倍数。

朴素欧几里得

对于求解 $\gcd(a, b)$ 咱们须要用朴素欧几里得定理。

$\gcd(a,b) = \gcd(b, a \bmod b)$ 。

这个是比较好证实的：

假设 $a = k*b+r$ ，有 $r = a \bmod b$ 。不妨设 $d$ 为 $a$ 和 $b$ 的一个任意一个公约数，则有 $a \equiv b \equiv 0 \pmod d$ 。
因为同余的性质 $a-kb \equiv r \equiv 0 \pmod d$ 所以 $d$ 是 $b$ 和 $a \bmod b$ 的公约数。
而后 $\forall~ d~|\gcd(a, b)$ 都知足这个性质，因此这个定理成立啦。

因此咱们就能够获得算 $\gcd$ 的一个简单函数啦。

int gcd(int a, int b) {
    return !b ? a : gcd(b, a % b);
}

这个复杂度是 $O(\log)$ 的，十分迅速。

而后断定是否有解后，咱们须要在这个基础上求一组解 $(x, y)$ ，因为 $a,b,c$ 都是 $\gcd(a,b)$ 的倍数。

对于 $a, b$ 有负数的状况，咱们须要将他们其中一个负数加上另一个数直到非负。（因为前面朴素欧几里得定理是不会影响的）两个负数，直接将整个式子反号，而后放到 $c$ 上就好了。

咱们将它们都除以 $\gcd(a, b)$ ，不影响后面的计算。

也就是咱们先求对于 $ax + by = 1$ 且 $a \bot b$ （互质）求 $(x, y)$ 的解。

接下来咱们利用前面的朴素欧几里得定律推一波式子。

\[ \begin{align} ax+by&=\gcd(a,b)\\ &=\gcd(b, a \bmod b) \\ &\Rightarrow bx + (a \bmod b)~ y \\ &= bx + (a - \lfloor \frac{a}{b} \rfloor ~b)~y \\ &= a y + b~(x - \lfloor \frac{a}{b} \rfloor ~y) \end{align} \]

不难发现此时 $x$ 变成了 $y$ ， $y$ 变成了 $\displaystyle x - \lfloor \frac{a}{b} \rfloor ~y$ ，利用这个性质，咱们能够递归的去求解 $(x,y)$ 。

边界条件其实和前面朴素欧几里得是同样的 $b=0$ 的时候，咱们有 $a=1,ax+by=1$ 那么此时 $x = 1, y = 0$ 。

这样作完的话咱们用 $O(\log)$ 的时间就会获得一组 $(x, y)$ 的特殊解。

解系扩展

但经常咱们要求对于 $x ~or~ y$ 的最小非负整数解，这个的话咱们须要将单个 $(x, y)$ 扩展成一个解系。

若是学过不定方程的话，就能够轻易获得这个解系的，在此不过多赘述了。
\[ d\in \mathbb{Z}\\ \begin{cases} x = x_0 + db \\ y = y_0 - da \\ \end{cases} \]

要记住，$(x, y)$ 都须要乘上 $c$ 。

而后咱们直接令 $x$ 为 $(x \bmod b + b) \bmod b$ 就好了。

代码实现

超简短版本：

递归调用的话， $y=x', x = y'$ ，只须要修改 $y$ 就好了。~~（是否是很好背啊）~~

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;
}

欧拉函数

定义

咱们定义 $\varphi (x)$ 为小于 $x$ 的正整数中与 $x$ 互质的数的个数，称做欧拉函数。数学方式表达就是
\[ \varphi(x) = \sum_{i < x} [i \bot x] \]

但须要注意，咱们定义 $\varphi(1) = 1$ 。

性质

若 $x$ 为质数，$\varphi(x) = x - 1$ 。

证实：这个很显然了，由于除了质数自己的数都与它互质。
若 $x = p^k$ （ $p$ 为质数， $x$ 为单个质数的整数幂），则 $\varphi(x) = (p - 1) \times p ^ {k - 1}$ 。

证实：不难发现全部 $p$ 的倍数都与 $x$ 不互质，其余全部数都与它互质。

$p$ 的倍数恰好有 $p^{k-1}$ 个（包括了 $x$ 自己）。

那么就有 $\varphi(x) = p^k - p^{k-1} = (p - 1) \times p^{k - 1}$ 。
若 $p, q$ 互质，则有 $\varphi(p \times q) = \varphi(p) \times \varphi(q)$ ，也就是欧拉函数是个积性函数。

证实：

若是 $a$ 与 $p$ 互质 $(a<p)$ ， $b$ 与 $q$ 互质 $(b<q)$ ， $c$ 与 $pq$ 互质 $(c<pq)$ ，则 $c$ 与数对 $(a,b)$ 是一一对应关系。因为 $a$ 的值有 $\varphi (p)$ 种可能， $b$ 的值有 $\varphi(q)$ 种可能，则数对 $(a,b)$ 有 $φ(p)φ(q)$ 种可能，而 $c$ 的值有 $φ(pq)$ 种可能，因此 $φ(pq)$ 就等于 $φ(p)φ(q)$ 。

具体来讲这一条须要 中国剩余定理 以及 彻底剩余系 才能证实，感性理解一下它的思路吧。~~（后面再填）~~
对于一个正整数 $x$ 的质数幂分解 $\displaystyle x = {p_1}^{k_1} \times {p_2}^{k_2} \times \cdots \times {p_{n} }^{k_n} = \prod _{i=1}^{n} {p_i}^{k_i}$ 。
\[\displaystyle \varphi(x) = x \times (1 - \frac{1}{p_1}) \times (1 - \frac{1}{p_2}) \times \cdots \times (1 - \frac{1}{p_n}) = x~\prod_{i=1}^{n} (1-\frac{1}{p_i})\]

证实：

咱们考虑用前几条定理一块儿来证实。
\[ \begin{align} \varphi(x) &= \prod_{i=1}^{n} \varphi(p_i^{k_i}) \\ &= \prod_{i=1}^{n} (p_i-1)\times {p_i}^{k_i-1}\\ &=\prod_{i=1}^{n} {p_i}^{k_i} \times(1 - \frac{1}{p_i})\\ &=x~ \prod_{i=1}^{n} (1- \frac{1}{p_i}) \end{align} \]
若 $p$ 为 $x$ 的约数（ $p$ 为质数, $x$ 为任意正整数），咱们有 $\varphi(x \times p) = \varphi(x) \times p $ 。

证实：

咱们利用以前的第 $4$ 个性质来证实就好了。

不难发现 $\displaystyle \prod _{i=1}^{n} (1 - \frac{1}{p_i})$ 是不会变的，前面的那个 $x$ 会变成 $x \times p$ 。

因此乘 $p$ 就好了。
若 $p$ 不是 $x$ 的约数（ $p$ 为质数, $x$ 为任意正整数），咱们有 $\varphi(x \times p) = \varphi(x) \times (p - 1)$ 。

证实：$p, x$ 互质，因为 $\varphi$ 积性直接获得。

求欧拉函数

求解单个欧拉函数

咱们考虑质因数分解，而后直接利用以前的性质 $4$ 来求解。

快速分解的话能够用前面讲的 $Pollard~Rho$ 算法。

求解一系列欧拉函数

首先须要学习线性筛，咱们将其改一些地方。

质数 $p$ 的 $\varphi(p) = p-1$ 。

对于枚举一个数 $i$ 和另一个质数 $p$ 的积 $x$ 。

咱们在线性筛，把合数筛掉会分两种状况。

$p$ 不是 $i$ 的一个约数，因为性质 $5$ 就有 $\varphi(x) = \varphi(i) \times (p - 1)$ 。
$p$ 是 $i$ 的一个约数，此时咱们会跳出循环，因为性质 $6$ 有 $\varphi(x) = \varphi(i) \times p$ 。

代码实现

bitset<N> is_prime; int prime[N], phi[N], cnt = 0;
void Init(int maxn) {
    is_prime.set(); is_prime[0] = is_prime[1] = false; phi[1] = 1;
    For (i, 2, maxn) {
        if (is_prime[i]) phi[i] = i - 1, prime[++ cnt] = i;
        For (j, 1, cnt) {
            int res = i * prime[j];
            if (res > maxn) break;
            is_prime[res] = false;
            if (i % prime[j]) phi[res] = phi[i] * (prime[j] - 1);
            else { phi[res] = phi[i] * prime[j]; break; }
        }
    }
}

扩展欧拉定理

\[ a^b\equiv \begin{cases} a^{b \bmod \varphi(p)} & a \bot p\\ a^b & a \not \bot p,b<\varphi(p)\\ a^{b \bmod \varphi(p) + \varphi(p)} & a \not \bot p,b\geq\varphi(p) \end{cases} \pmod p \]

证实看这里，有点难证，记结论算啦qwq

这个经常能够用于降幂这种操做。它的一个扩展就是最前面的那个 费马小定理 。

求解模线性方程组

问题描述

给定了 $n$ 组除数 $m_i$ 和余数 $r_i$ ，经过这 $n$ 组 $(m_i,r_i)$ 求解一个 $x$ ，使得 $x \bmod m_i = r_i$ 这就是 模线性方程组 。

数学形式表达就是：
\[ \begin{cases} x \equiv r_1 \pmod {m_1}\\ x \equiv r_2 \pmod {m_2}\\ ~~~~\vdots \\ x \equiv r_n \pmod {m_n} \end{cases} \]

求解一个 $x$ 知足上列全部方程。

算法解决

中国剩余定理

中国剩余定理提供了一个较为通用的解决方法。（彷佛是惟一一个以中国来命名的定理）

若是 $m_1, m_2, \dots , m_n$ 两两互质，则对于任意的整数 $r_1, r_2 , \dots , r_n$ 方程组 $(S)$ 有解，而且能够用以下方法来构造：

令 $\displaystyle M = \prod_{i=1} ^ {n} m_i$ ，并设 $\displaystyle M_i = \frac{M}{m_i}$ ，令 $t_i$ 为 $M_i$ 在模 $m_i$ 意义下的逆元（也就是 $t_i \times M_i \equiv 1 \pmod {m_i}$ ）。

不难发现这个 $t_i$ 是必定存在的，由于因为前面要知足两两互质，那么 $M_i \bot m_i$ ，因此必有逆元。

那么在模 $M$ 意义下，方程 $(S)$ 有且仅有一个解： $\displaystyle x = (\sum_{i=1}^{n} r_i t_i M_i) \bmod M$ 。

不难发现这个特解是必定知足每个方程的，只有一个解的证实能够考虑特解 $x$ 对于第 $i$ 个方程，下个继续知足条件的 $x$ 为 $x + m_i$ 。那么对于总体来讲，下个知足条件的数就是 $x + \mathrm{lcm} (m_1, m_2, \dots, m_n) = x + M$ 。

严谨数学证实请见百度百科。

利用扩欧合并方程组

咱们考虑合并方程组，好比从 $n=2$ 开始递推。
\[ \begin{cases} x \equiv r_1 \pmod {m_1}\\ x \equiv r_2 \pmod {m_2} \end{cases} \]

也就等价于
\[ \begin{cases} x = m_1 \times k_1 + r_1 \\ x = m_2 \times k_2 + r_2 \\ \end{cases} \]

此处 $k_1, k_2 \in \mathbb{N}$ 。联立后就获得了以下一个方程：
\[ \begin{align} m_1 \times k_1 + r_1 &= m_2 \times k_2 + r_2 \\ m_1 \times k_1 - m_2 \times k_2 &= r_2 - r_1 \end{align} \]

咱们令 $a = m_1, b = m_2, c = r_2 - r_1$ 就变成了 $ax+by=c$ 的形式，用以前讲过的扩展欧几里得，能够求解。

首先先判断有无解。假设存在解，咱们先解出一组 $(k_1,k_2)$ ，而后带入解出 $x=x_0$ 的一个特解。

咱们将这个能够扩展成一个解系：
\[ x = x_0 + k \times \mathrm{lcm} (m_1,m_2) ~,~ k \in \mathbb{N} \]

因为前面不定方程的结论， $k_1$ 与其相邻解的间距为 $\displaystyle \frac{m_2}{\gcd(m_1,m_2)}$ ，又有 $x=m_1 \times k_1 + r_1$ 。因此 $x$ 与其相邻解的距离为 $\displaystyle \frac{m_1 m_2}{\gcd(m_1,m_2)} = \mathrm{lcm}(m_1,m_2)$ 。

因此咱们令 $M=\mathrm{lcm}(m_1, m_2), R = x_0$ 则又有新的模方程 $x \equiv R \pmod M$ 。

而后咱们就将两个方程合并成一个了，只要不断重复这个过程就能作完了。

这个比 中国剩余定理 优秀的地方就在于它这个不须要要求 $m$ 两两互质，而且能够较容易地判断无解的状况。

代码实现

注意有些细节，好比求两个数的 $\mathrm{gcd}$ 的时候，必定要先除再乘，防止溢出！！

ll mod[N], rest[N];
ll Solve() {
    For (i, 1, n - 1) {
        ll a = mod[i], b = mod[i + 1], c = rest[i + 1] - rest[i], gcd = __gcd(a, b), k1, k2;
        if (c % gcd) return - 1;
        a /= gcd; b /= gcd; c /= gcd;
        Exgcd(a, b, k1, k2);

        k1 = ((k1 * c) % b + b) % b;
        mod[i + 1] = mod[i] / __gcd(mod[i], mod[i + 1]) * mod[i + 1] ;
        rest[i + 1] = (mod[i] * k1 % mod[i + 1] + rest[i]) % mod[i + 1];
    }
    return rest[n];
}

扩展卢卡斯

问题描述

求
\[ {n \choose m} \bmod p \]

$1 \le m \le n \le 10^{18}, 2 \le p \le 10^6$ 不保证 $p$ 为质数。

问题求解

虽然说是扩展卢卡斯，可是和卢卡斯定理没有半点关系。

用了几个性质，首先能够考虑 $p$ 分解质因数。

假设分解后
\[ p = \prod_i {p_i}^{k_i} \]

咱们能够对于每一个 $p_i$ 单独考虑，假设咱们求出了 $\displaystyle {n \choose m} \bmod {{p_i}^{k_i}}$ 的值。

而后咱们能够考虑用前文提到的 $CRT$ 来进行合并（须要用扩欧求逆元）。

问题就转化成如何求 $\displaystyle {n \choose m} \bmod {{p_i}^{k_i}}$ 了。

首先咱们考虑它有多少个关于 $p_i$ 的因子（也就是多少次方）。

咱们令 $f(n)$ 为 $n!$ 中包含 $p_i$ 的因子数量，那么 $\displaystyle {n \choose m} = \frac{n!}{m!(n - m)!}$ ，因此因子数量就为 $f(n) - f(m) - f(n - m)$ 。

那如何求 $f(n)$ 呢。

咱们举出一个很简单的例子来讨论。

$n=19,p_i=3,k_i=2$ 时:

\[ \begin{align} n!&=1*2*3*\cdots*19\\ &=(1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19)∗3^6∗6!\\ &\equiv (1*2*4*5*7*8)^2*19*3^6*6!\\ &\equiv (1∗2∗4∗5∗7∗8)^2∗19∗{3^6}∗6! \pmod {3^2} \\ \end{align} \]

不难发现每次把 $p_i$ 倍数提出来的东西，就是 $\displaystyle \lfloor \frac{n}{p_i} \rfloor !$ ，那么很容易获得以下递推式：

\[ f(n) = f(\displaystyle \lfloor \frac{n}{p_i} \rfloor) + \displaystyle \lfloor \frac{n}{p_i} \rfloor \]

不难发现这个求单个的复杂度是 $O(\log n)$ 的，十分优秀。

而后考虑剩下与 $p_i$ 互不相关的如何求，不难发如今 ${p_i}^{k_i}$ 意义下会存在循环节（好比前面的 $(1∗2∗4∗5∗7∗8)^2$ ，可能最后多出来一个或者多个不存在循环节中的数。

不难发现循环节长度是 $< {p_i}^{k_i}$ ，由于同余的在 $> {p_i}^{k_i}$ 以后立刻出现，而后把多余的 $\displaystyle \lfloor \frac{n}{p_i} \rfloor$ 的部分递归处理就好了。

而后就能够快速处理出与 $p_i$ 无关的了，最后合并一下就好了。

简介算法流程：

对于每一个 ${p_i}^{k_i}$ 单独考虑。
- 用 $f(n)$ 处理出有关于 $p_i$ 的因子个数。
- 而后递归处理出无关 $p_i$ 的组合数的答案。
- 把这两个乘起来合并便可。
最后用 $CRT$ 合并全部解便可。

瞎分析的复杂度（ Great_Influence 博客上抄的）。（不可靠）
\[ O(\sum {p_i}^{k_i} \log n(\log_{p_i} n - k)+p_i \log p) \sim O(p \log p) \]

代码解决

#include <bits/stdc++.h>

#define For(i, l, r) for(register int i = (l), i##end = (int)(r); i <= i##end; ++i)
#define Fordown(i, r, l) for(register int i = (r), i##end = (int)(l); i >= i##end; --i)
#define Set(a, v) memset(a, v, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define debug(x) cout << #x << ": " << (x) << endl
#define DEBUG(...) fprintf(stderr, __VA_ARGS__)

using namespace std;

typedef long long ll;

template<typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, T b) { return b > a ? a = b, 1 : 0; }

void File() {
#ifdef zjp_shadow
    freopen ("P4720.in", "r", stdin);
    freopen ("P4720.out", "w", stdout);
#endif
}

ll n, m, 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;
}

inline ll fpm(ll x, ll power, ll Mod) {
    ll res = 1;
    for (; power; power >>= 1, (x *= x) %= Mod)
        if (power & 1) (res *= x) %= Mod;
    return res;
}

inline ll fac(ll n, ll pi, ll pk) {
    if (!n) return 1;
    ll res = 1;
    For (i, 2, pk) if (i % pi) (res *= i) %= pk;
    res = fpm(res, n / pk, pk);
    For (i, 2, n % pk) if (i % pi) (res *= i) %= pk;
    return res * fac(n / pi, pi, pk) % pk;
}

inline ll Inv(ll n, ll Mod) {
    ll x, y; Exgcd(n, Mod, x, y);
    return (x % Mod + Mod) % Mod;
}

inline ll CRT(ll b, ll Mod) {
    return b * Inv(p / Mod, Mod) % p * (p / Mod) % p;
}

inline ll factor(ll x, ll Mod) {
    return x ? factor(x / Mod, Mod) + (x / Mod) : 0;
}

inline ll Comb(ll n, ll m, ll pi, ll pk) {
    ll k = factor(n, pi) - factor(m, pi) - factor(n - m, pi);
    if (!fpm(pi, k, pk)) return 0;
    return fac(n, pi, pk) * Inv(fac(m, pi, pk), pk) % pk * Inv(fac(n - m, pi, pk), pk) % pk * fpm(pi, k, pk) % pk;
}

inline ll ExLucas(ll n, ll m) {
    ll res = 0, tmp = p;
    For (i, 2, sqrt(p + .5)) if (!(tmp % i)) {
        ll pk = 1; while (!(tmp % i)) pk *= i, tmp /= i;
        (res += CRT(Comb(n, m, i, pk), pk)) %= p;
    }
    if (tmp > 1) (res += CRT(Comb(n, m, tmp, tmp), tmp)) %= p;
    return res;
}

int main () {

    File();

    cin >> n >> m >> p;

    printf ("%lld\n", ExLucas(n, m));

    return 0;

}