【快速因數分解】Pollard's Rho 演算法

**Pollard-Rho** 是一個很神奇的演算法，用於在 $O(n^{\frac{1}4}) $的期望時間複雜度內計算合數 n 的某個非平凡因子（除了1和它本身以外能整除它的數）。事書上給出的複雜度是 $O(\sqrt{p})$ ， p 是 n 的某個最小因子，滿足 p 與 n/p 互質。雖然是隨機的，但 **Pollard Rho** 演算法在實際環境中執行的相當不錯，不會被卡。 > 簡單來說：**Pollard-Rho演算法**是 **John Pollard**發明的一種能 快速找到大整數的一個非1、非自身的因子 的演算法。 ## 字首知識 1. [快速乘](https://www.cnblogs.com/RioTian/p/13927813.html) 2. [Miller-Rabin 素數檢驗演算法](https://www.cnblogs.com/RioTian/p/13927952.html) ## PollardRhoPollardRho 演算法： ### 問題模型： 給一個數 $n$ ，你需要快速求出它的一個非平凡因子。 對於一個比較小的數（ $n≤10^9$ ），我們直接暴力列舉質因子就行但太大了我們就必須考慮一下隨機演算法了。 對於一個非常大的合數 $n≤10^{18}$ （如果可能是質數，我們可以用**Miller Rabin**判一下）我們要求 nn 的某一個非平凡因子，如果 $n$ 的質因子很多（就是約數很多）我們也可以暴力隨機弄一弄，但如果是一些（像由兩個差不多的而且很大的質數乘得的 $n$ ）它的非平凡因子就只有兩個而且 $n$ 本身還很大，此時如果我們隨機的話複雜度 $O(n)$ ，這個太難接受了，所以我們想辦法優化一下。 ### 1.我們求 gcd 如果我們直接隨機求 $n$ 的某一個約數複雜度很高，我們可以考慮求一下 gcdgcd 。因為我們可以保證一個合數它絕對存在一個質因子小於 $√n$ ，所以在 $n$ 就存在至少 $√n$ 個數與 $n$ 有大於一的公約數。於是我們隨機的可能性就提高到了 $O(\sqrt{n})$ 。 ### 2.玄學隨機法： $x^2+c \ \ _{x=rand(),c=rand()}$ 其實網上大多都叫這種方法為生日悖論（悖論是不符合經驗但符合事實的結論），有興趣的可以去看看。但博主也不懂覺得這個方法跟我們的 **PollardRho** 的關聯。在1中我們用gcd的方案，現在我們考慮有沒有辦法能夠快速的找到這些與 $n$ 有公約數的數（或者說找到一種隨機生成方法，使隨機生成這類我們需要的數的概率變大一點）。 這個隨機生成方法最終被 $Pollard$ 研發出來了：就是標題那個式子！ 我們定義 $$(y=x\ x=x×x+c \ g=gcd(y−x,p))$$ 為一次操作，我們發現 $g>1$ （就是找到我們需要的數）的機率出奇的高，根據國際上大佬的瘋狂測試，在 $n$ 中找到一對 $y−x$ 使兩者有公約數的概率接近 $n^{\frac14}$ ，不得不說太玄學又太優秀了！ ### 判環 這種隨機生成方法雖然優秀，但也有一些需要注意之處，比如有可能會生成一個環，並不斷在這個環上生成以前生成過一次的數，所以我們必須寫點東西來判環： 1. 我們可以讓 $y$ 根據生成公式以比 $x$ 快兩倍的速度向後生成，這樣當$y$ 再次與 $x$ 相等時，$x$ 一定跑完了一個圈且沒重複跑多少！ 2. 我們可以用倍增的思想，讓 $y$ 記住 $x$ 的位置，然後 $x$ 再跑當前跑過次數的一倍的次數。這樣不斷讓 $y$ 記住 $x$ 的位置，$x$ 再往下跑，因為倍增所以當 $x$ 跑到 $y$ 時，已經跑完一個圈，並且也不會多跑太多（這個必須好好想一想，也可以看看程式碼如何實現的）（因為這種方法會比第一種用的更多，因為它在 $PollardRho$ 二次優化時還可以起到第二個作用）（這裡就給第二種的程式碼吧） ```cpp inline ll rho(ll n) { ll x, y, c, g; rg i, j; x = y = rand(); c = rand(); //初始化 i = 0, j = 1; //倍增初始化 while (++i) { x = (ksc(x, x, n) + c) % n; // x只跑一次 if (x == y) break; //跑完了一個環 g = gcd(Abs(y - x), n); //求gcd（注意絕對值） if (g > 1) return g; if (i == j) y = x, j <<= 1; //倍增的實現 } } ``` ## PollardRho 演算法的二次優化： 我們發現其實 $PollardRho$ 演算法其實複雜度不止 $O(n^{\frac14})$ ，它每一次操作都會求一次 $gcd$ ，所以複雜度會變成 $O(n^{\frac14}log)$ 我們發現這個 $log$ 實在有些拖慢速度（雖然實際上也很快了）。於是一個很奇妙的優化誕生了！ **二次優化**：我們發現我們在每一次生成操作中，乘法之後所模的數就是我們的 $n$ ，而我們要求的就是 $n$ 的某一個約數！也就是說我們現在的模數並不是一個質數！而根據取模的性質：如果模數和被模的數都含有一個公約數，那麼這次模運算的結果必然也會是這個公約數的倍數！！！所以如果我們將若干個 $(y−x)$ 乘到一起，因為中途模的是 $n$ ，所以如果我的若干個 $(y−x)$ 中有一個與 $n$ 有公約數，最後的結果定然也會含有這個公約數！所以我們完全可以多算幾次 $(y−x)$ 的乘積在來求 $gcd$ （一般連續算127次再求一次 $gcd$ （這個應該有大佬測試過了）），這樣我們的複雜度就能降一個 $log$ 級別，跑的飛快！ **對原演算法的一些影響：**任何一個優化都是要考慮其對原演算法的影響的。這個優化也會有一些影響：首先，因為我們會等大概127次後再去 $gcd$ ，然而我們有可能在生成時碰上一個環，我們有可能還沒生成127次就跳出這個環了，這樣就無法得出答案；其次，我們的 $PollardRho$ 演算法實在太玄學優秀了，可能我們跑127次之後，所有 $(y−x)$ 的乘積就變成了 $n$ 的倍數（模 $n$ 意義下得到 $0$ ），所以我們不能完全就只呆板的等127次在取模。 **那怎麼辦呢？**（在上文就已經提到過了：倍增時我們的 $i$ 和 $j$ 會在二次優化時起到另一種作用。）沒錯！就是倍增！我們可以用倍增，分別在生成（1次，2次，4次，8次，16次，32次......）之後進行 $gcd$ ！這樣就可以完全避免上述的兩個影響，而且我們是倍增來gcd的這對我們的複雜度影響不大！！！！！然後我們看程式碼吧！ ```cpp inline ll rho(ll p) { //求出p的非平凡因子 ll x, y, z, c, g; rg i, j; //先擺出來（z用來存（y-x）的乘積） while (1) { //保證一定求出一個因子來 y = x = rand() % p; //隨機初始化 z = 1; c = rand() % p; //初始化 i = 0, j = 1; //倍增初始化 while (++i) { //開始玄學生成 x = (ksc(x, x, p) + c) % p; //可能要用快速乘 z = ksc(z, Abs(y - x), p); //我們將每一次的（y-x）都累乘起來 if (x == y || !z) break; //如果跑完了環就再換一組試試（注意當z=0時，繼續下去是沒意義的） if (!(i % 127) || i == j) { //我們不僅在等127次之後gcd我們還會倍增的來gcd g = gcd(z, p); if (g > 1) return g; if (i == j) y = x, j <<= 1; //維護倍增正確性，並判環（一箭雙鵰） } } } } ``` 然後是兩道道例題，都有點難度的。 ## code1： ```cpp #include using namespace std; #define rg register int typedef long long ll; typedef unsigned long long ull; typedef long double lb; int n; ll ans, m; inline ll qr() { //快讀 char ch; while ((ch = getchar()) < '0' || ch > '9') ; ll res = ch ^ 48; while ((ch = getchar()) >= '0' && ch <= '9') res = res * 10 + (ch ^ 48); return res; } inline ll Abs(ll x) { return x < 0 ? -x : x; } //取絕對值 inline ll gcd(ll x, ll y) { //非遞迴求gcd ll z; while (y) { z = x; x = y; y = z % y; } return x; } inline ll ksc(ull x, ull y, ll p) { // O(1)快速乘（防爆long long） return (x * y - (ull)((lb)x / p * y) * p + p) % p; } inline ll ksm(ll x, ll y, ll p) { //快速冪 ll res = 1; while (y) { if (y & 1) res = ksc(res, x, p); x = ksc(x, x, p); y >>= 1; } return res; } inline bool mr(ll x, ll p) { // mille rabin判質數 if (ksm(x, p - 1, p) != 1) return 0; //費馬小定理 ll y = p - 1, z; while (!(y & 1)) { //二次探測 y >>= 1; z = ksm(x, y, p); if (z != 1 && z != p - 1) return 0; if (z == p - 1) return 1; } return 1; } inline bool prime(ll x) { if (x < 2) return 0; // mille rabin判質數 if (x == 2 || x == 3 || x == 5 || x == 7 || x == 43) return 1; return mr(2, x) && mr(3, x) && mr(5, x) && mr(7, x) && mr(43, x); } inline ll rho(ll p) { //求出p的非平凡因子 ll x, y, z, c, g; rg i, j; //先擺出來（z用來存（y-x）的乘積） while (1) { //保證一定求出一個因子來 y = x = rand() % p; //隨機初始化 z = 1; c = rand() % p; //初始化 i = 0, j = 1; //倍增初始化 while (++i) { //開始玄學生成 x = (ksc(x, x, p) + c) % p; //可能要用快速乘 z = ksc(z, Abs(y - x), p); //我們將每一次的（y-x）都累乘起來 if (x == y || !z) break; //如果跑完了環就再換一組試試（注意當z=0時，繼續下去是沒意義的） if (!(i % 127) || i == j) { //我們不僅在等127次之後gcd我們還會倍增的來gcd g = gcd(z, p); if (g > 1) return g; if (i == j) y = x, j <<= 1; //維護倍增正確性，並判環（一箭雙鵰） } } } } inline void prho(ll p) { //不斷的找他的質因子 if (p <= ans) return; //最優性剪紙 if (prime(p)) { ans = p; return; } ll pi = rho(p); //我們一次必定能求的出一個因子，所以不用while while (p % pi == 0) p /= pi; //記得要除盡 return prho(pi), prho(p); //分開繼續分解質因數 } int main() { // freopen("in.txt","r",stdin); // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0); n = qr(); srand(time(0)); //隨機數生成必備！！！ for (rg i = 1; i <= n; ++i) { ans = 1; prho((m = qr())); if (ans == m) puts("Prime"); else printf("%lld\n", ans); } return 0; } ``` ## code2： 題目連結：[杭電3864 D_num](http://acm.hdu.edu.cn/showproblem.php?pid=3864) 根據這道題的答案性質，我們知道這個數在除了1以外還有三個約數，仔細想一下不難發現這類數最多隻能有兩個質因子，而這兩個質因子分別就是它的兩個約數，而第三個約數自然就是它自己了！（當然還有可能出現質數的三次方這類特殊的數，所以我們要特判一下） 所以我們的思路就是，先求出 $n$ 的一個小於它的約數，然後我們可以用這個約數算出 $n$ 的另一個約數，然後我們判一下它是不是質數的三次方這類特殊的數（是就可以直接輸出答案了），不是的話我們判斷這兩個約數是否分別是不同的質因子，是就能輸出了，不是就說明這不是 $D\_num$ ```cpp #include using namespace std; #define ll long long #define db double #define rg register int #define cut {puts("is not a D_num");continue;} //這個可以看一看 ll n; inline ll qr() { char ch; while ((ch = getchar()) < '0' || ch > '9') ; ll res = ch ^ 48; while ((ch = getchar()) >= '0' && ch <= '9') res = res * 10 + (ch ^ 48); return res; } inline ll Abs(ll x) { return x < 0 ? -x : x; } inline ll gcd(ll x, ll y) { ll z; while (y) { z = x, x = y, y = z % y; } return x; } inline ll ksm(ll x, ll y, ll p) { ll res = 1; while (y) { if (y & 1) res = (__int128)res * x % p; x = (__int128)x * x % p; y >>= 1; } return res; } inline bool mr(ll x, ll p) { if (ksm(x, p - 1, p) != 1) return 0; ll y = p - 1, z; while (!(y & 1)) { y >>= 1; z = ksm(x, y, p); if (z != 1 && z != p - 1) return 0; if (z == p - 1) return 1; } return 1; } inline bool prime(ll p) { if (p < 2) return 0; if (p == 2 || p == 3 || p == 5 || p == 7 || p == 43) return 1; return mr(2, p) && mr(3, p) && mr(5, p) && mr(7, p) && mr(43, p); } inline ll rho(ll p) { ll x, y, z, c, g; rg i, j; while (1) { y = x = rand() % p; z = 1, c = rand() % p; i = 0, j = 1; while (++i) { x = ((__int128)x * x + c) % p; z = (__int128)z * Abs(y - x) % p; if (x == y) break; if (i % 72 || i == j) { g = gcd(z, p); if (g > 1 && g != p) return g; if (i == j) y = x, j <<= 1; } } } } int main() { // freopen("in.txt","r",stdin); // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0); srand(time(0)); ll p, q; while (scanf("%lld", &n) != EOF) { if (prime(n) || n <= 2) cut; p = rho(n); q = n / p; if ((__int128)p * p == q || (__int128)q * q == p) //這題很卡細節的 { printf("%lld %lld %lld\n", p, q, n); continue; } if (!prime(p) || !prime(q) || p == q) cut; //注意兩個相等也不行 printf("%lld %lld %lld\n", min(p, q), max(p, q), n); } return 0; } ``` ## 參考 https://blog.csdn.net/qq_39972971/article/details/82346390 https://ld246.com/article/1569220928268 https://zhuanlan.zhihu.com/p/267884783