這個復習沒有什么順序的... 想到什么復習什么而已qwq
Miller-Rabin 質數測試
問題描述
測試一個正整數 \(p\) 是否為素數，需要較快且准確的測試方法。\((p \le 10^{18})\)
算法解決
費馬小定理
首先需要了解 費馬小定理 。
費馬小定理：對於質數 \(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\) 測試相比，增加了一個 二次探測定理 。
二次探測定理：如果 \(p\) 是奇素數，則 \(x^2 \equiv 1 \pmod p\) 的解為 \(x \equiv 1\) 或 \(x \equiv p - 1 \pmod p\) 。
這是很容易證明的：
\[\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}\) 。
那么概率就是
假設我們帶入 \(n\) 等於 \(60\) ，那么這個概率就是 \(0.0058\) ，也就是說基本上不可能發生。
好像是因為和大眾的常識有些違背，所以稱作生日悖論。
假設一年有 \(N\) 天，那么只要 \(n \ge \sqrt N\) 存在相同生日的概率就已經大於 \(50 \%\) 了。
利用偽隨機數判斷
本段參考了 Doggu 大佬的博客 。
我們考慮利用之前那個性質來判斷，生成了許多隨機數，然后兩兩枚舉判斷。
假設枚舉的兩個數為 \(i, j\) ，假設 \(i > j\) ，那么判斷一下 \(\gcd(i - j, n)\) 是多少，如果非 \(1\) ，那么我們就已經找出了一個 \(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\) 輪）
我們考慮用數學歸納法證明：
顯然當只有 \(1\) 個候選人時，該候選人就是當選者，並且他的編號為 \(0\) 。
假設我們已經求解出了 \(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)\) 的解。
接下來我們利用前面的朴素歐幾里得定律推一波式子。
不難發現此時 \(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)\) 擴展成一個解系。
如果學過不定方程的話，就可以輕易得到這個解系的，在此不過多贅述了。
要記住，\((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(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; }
}
}
}
擴展歐拉定理
證明看 這里 ，有點難證，記結論算啦qwq
這個常常可以用於降冪這種操作。它的一個擴展就是最前面的那個 費馬小定理 。
求解模線性方程組
問題描述
給定了 \(n\) 組除數 \(m_i\) 和余數 \(r_i\) ，通過這 \(n\) 組 \((m_i,r_i)\) 求解一個 \(x\) ，使得 \(x \bmod m_i = r_i\) 這就是 模線性方程組 。
數學形式表達就是 ：
求解一個 \(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\) 開始遞推。
也就等價於
此處 \(k_1, k_2 \in \mathbb{N}\) 。聯立后就得到了如下一個方程：
我們令 \(a = m_1, b = m_2, c = r_2 - r_1\) 就變成了 \(ax+by=c\) 的形式，用之前講過的 擴展歐幾里得 ，可以求解。
首先先判斷有無解。假設存在解，我們先解出一組 \((k_1,k_2)\) ，然后帶入解出 \(x=x_0\) 的一個特解。
我們將這個可以擴展成一個解系：
由於前面不定方程的結論， \(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];
}
擴展盧卡斯
問題描述
求
\(1 \le m \le n \le 10^{18}, 2 \le p \le 10^6\) 不保證 \(p\) 為質數。
問題求解
雖說是擴展盧卡斯，但是和盧卡斯定理沒有半點關系。
用了幾個性質，首先可以考慮 \(p\) 分解質因數。
假設分解后
我們可以對於每個 \(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\) 時:
不難發現每次把 \(p_i\) 倍數提出來的東西，就是 \(\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 博客上抄的）。（不可靠）
代碼解決
#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;
}