基礎數論複習

時間 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;

}