Pollard's Rho algorithm

**Pollard-Rho** It's a magic algorithm , Used in $O(n^{\frac{1}4}) $ In the expected time complexity of n Some nontrivial factor of （ except 1 And the number that can be divisible by itself ）. The complexity given in the story book is $O(\sqrt{p})$ , p yes n Some of the smallest factors of , Satisfaction p And n/p Mutuality . Although it's random , but **Pollard Rho** The algorithm performs quite well in the actual environment , It won't get stuck . > In short ：**Pollard-Rho Algorithm ** yes **John Pollard** An invention that can Find a non of large integer quickly 1、 Factors that are not themselves The algorithm of . ## Prefix knowledge 1. [ Fast ride ](https://www.cnblogs.com/RioTian/p/13927813.html) 2. [Miller-Rabin Prime number test algorithm ](https://www.cnblogs.com/RioTian/p/13927952.html) ## PollardRhoPollardRho Algorithm ： ### Problem model ： Give me a number $n$ , You need to quickly find a nontrivial factor of it . For a smaller number （ $n≤10^9$ ）, Let's just enumerate the prime factor directly, but if it's too big, we have to consider random algorithms . For a very large sum $n≤10^{18}$ （ If it could be prime numbers , We can use **Miller Rabin** Judge ） We ask for nn A nontrivial factor of , If $n$ There are many qualitative factors of （ It's a lot of approximations ） We can also do random violence , But if it's some （ As if multiplied by two similar and large primes $n$ ） It has only two nontrivial factors and $n$ It's still big , At this point, if we're random, the complexity is $O(n)$ , This is too hard to accept , So we try to optimize . ### 1. We ask gcd If we just do it at random $n$ One of the divisors of is very complex , We can consider asking for gcdgcd . Because we can guarantee that a composite number has an absolute prime factor less than $√n$ , So in $n$ There is, at least $√n$ Number and $n$ There is a common divisor greater than one . So the probability of our random increase to $O(\sqrt{n})$ . ### 2. The random method of metaphysics ： $x^2+c \ \ _{x=rand(),c=rand()}$ In fact, most people on the Internet call this method birthday paradox （ A paradox is a conclusion that does not conform to experience but conforms to facts ）, Those who are interested can go and have a look . But bloggers don't understand. They think this method is similar to ours **PollardRho** The connection of . stay 1 We use gcd The plan , Now let's consider whether there is a way to quickly find these and $n$ A number that has a common divisor （ Or find a random generation method , So that the probability of randomly generating the numbers we need is a little bit larger ）. This random generation method was eventually $Pollard$ Developed ： It's the formula of the title ！ We define $$(y=x\ x=x×x+c \ g=gcd(y−x,p))$$ For an operation , We found that $g>1$ （ Is to find the number we need ） The odds are surprisingly high , According to the crazy test of the world's bigwigs , stay $n$ We found a couple of $y−x$ So that the probability that they have a common divisor is close to $n^{\frac14}$ , I have to say that it's too metaphysical and so excellent ！ ### Judge the ring This random generation method is excellent , But there are also some things to be aware of , For example, it is possible to generate a ring , And generate numbers on this ring that have been generated once before , So we have to write something to judge the ring ： 1. We can make $y$ According to the generating formula, compare $x$ Double the speed of backward generation , It's like this $y$ Again with $x$ When equal ,$x$ Must have run a lap and didn't repeat much ！ 2. We can use the idea of multiplication , Let $y$ Remember $x$ The location of , And then $x$ Run twice as many times as you've run . This keeps making $y$ Remember $x$ The location of ,$x$ Run down again , Because of multiplication, so when $x$ Run to $y$ When , I've run a lap , And don't run too much （ You have to think about it , You can also see how the code implements ）（ Because this method will be used more than the first one , Because it's in $PollardRho$ Second optimization can also play a second role ）（ Here is the code of the second kind ） ```cpp inline ll rho(ll n) { ll x, y, c, g; rg i, j; x = y = rand(); c = rand(); // initialization i = 0, j = 1; // Multiply initialization while (++i) { x = (ksc(x, x, n) + c) % n; // x Only once if (x == y) break; // Finish a loop g = gcd(Abs(y - x), n); // seek gcd（ Pay attention to the absolute value ） if (g > 1) return g; if (i == j) y = x, j <<= 1; // The realization of multiplication } } ``` ## PollardRho Second optimization of the algorithm ： We found out that $PollardRho$ The complexity of the algorithm is more than $O(n^{\frac14})$ , Every time it operates, it asks for $gcd$ , So complexity becomes $O(n^{\frac14}log)$ We found this $log$ It's a bit slow （ Although it's actually fast ）. So a wonderful optimization was born ！ ** Second optimization **： We found that in every build operation , The number of the module after multiplication is our $n$ , And what we're asking for is $n$ A divisor of ！ In other words, our modulus is not a prime number ！ And according to the nature of the module ： If both modulus and modulus contain a common divisor , Then the result of this modular operation must also be a multiple of the common divisor ！！！ So if we take a number of $(y−x)$ Ride together , Because the half way model is $n$ , So if some of my $(y−x)$ One of them is related to $n$ There is a common divisor , The final result will certainly contain this common divisor ！ So we can count more times $(y−x)$ The product of is to find $gcd$ （ Generally, continuous calculation 127 Ask again and again $gcd$ （ This should have been tested by the big guy ））, So we can reduce our complexity by one $log$ Level , Run fast ！ ** Some effects on the original algorithm ：** Any optimization must consider its impact on the original algorithm . This optimization will also have some impact ： First , Because we'll wait about 127 I'll go again $gcd$ , However, it is possible for us to run into a ring during generation , We probably haven't generated 127 I jumped out of the ring , So there's no answer ; secondly , Our $PollardRho$ The algorithm is so metaphysical , Maybe we run 127 After that , all $(y−x)$ And the product of that becomes $n$ Multiple of （ model $n$ In the sense of getting $0$ ）, So we can't just wait for that 127 The second is taking the mold . ** So what to do ？**（ As mentioned above ： Doubling our $i$ and $j$ It will play another role in secondary optimization .） That's right ！ It's multiplication ！ We can use multiplication , They are generating （1 Time ,2 Time ,4 Time ,8 Time ,16 Time ,32 Time ......） After that $gcd$ ！ In this way, the two effects mentioned above can be completely avoided , And we're doubling it gcd This has little effect on our complexity ！！！！！ Then let's look at the code ！ ```cpp inline ll rho(ll p) { // Find out p The nontrivial factor of ll x, y, z, c, g; rg i, j; // Put it out first （z To save （y-x） The product of ） while (1) { // Make sure you find a factor y = x = rand() % p; // Random initialization z = 1; c = rand() % p; // initialization i = 0, j = 1; // Multiply initialization while (++i) { // The beginning of metaphysics x = (ksc(x, x, p) + c) % p; // Maybe you need to take the express ride z = ksc(z, Abs(y - x), p); // We will each time （y-x） I'm tired. I'm tired if (x == y || !z) break; // If you're done, try another group （ Pay attention to when z=0 When , There's no point in going on ） if (!(i % 127) || i == j) { // We're not just waiting for 127 After that gcd We're going to double that gcd g = gcd(z, p); if (g > 1) return g; if (i == j) y = x, j <<= 1; // Maintenance of multiplication accuracy , And judge the ring （ Kill two birds with one stone ） } } } } ``` Then there are two examples , It's a little difficult . ## 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() { // Read 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; } // Take the absolute value inline ll gcd(ll x, ll y) { // Non recursive seeking 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) Fast ride （ Explosion proof long long） return (x * y - (ull)((lb)x / p * y) * p + p) % p; } inline ll ksm(ll x, ll y, ll p) { // Fast power 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 Judge prime number if (ksm(x, p - 1, p) != 1) return 0; // Fermat's small theorem ll y = p - 1, z; while (!(y & 1)) { // Second detection 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 Judge prime number 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) { // Find out p The nontrivial factor of ll x, y, z, c, g; rg i, j; // Put it out first （z To save （y-x） The product of ） while (1) { // Make sure you find a factor y = x = rand() % p; // Random initialization z = 1; c = rand() % p; // initialization i = 0, j = 1; // Multiply initialization while (++i) { // The beginning of metaphysics x = (ksc(x, x, p) + c) % p; // Maybe you need to take the express ride z = ksc(z, Abs(y - x), p); // We will each time （y-x） I'm tired. I'm tired if (x == y || !z) break; // If you're done, try another group （ Pay attention to when z=0 When , There's no point in going on ） if (!(i % 127) || i == j) { // We're not just waiting for 127 After that gcd We're going to double that gcd g = gcd(z, p); if (g > 1) return g; if (i == j) y = x, j <<= 1; // Maintenance of multiplication accuracy , And judge the ring （ Kill two birds with one stone ） } } } } inline void prho(ll p) { // Constantly looking for his quality factor if (p <= ans) return; // Optimal paper cutting if (prime(p)) { ans = p; return; } ll pi = rho(p); // We can find one factor at a time , So there's no need to while while (p % pi == 0) p /= pi; // Remember to eliminate everything return prho(pi), prho(p); // Separate and continue to decompose the prime factor } int main() { // freopen("in.txt","r",stdin); // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0); n = qr(); srand(time(0)); // Random number generation is necessary ！！！ 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： Title Link ：[ Hang Dian 3864 D_num](http://acm.hdu.edu.cn/showproblem.php?pid=3864) According to the nature of the answer to this question , We know that this number is in addition to 1 There are three divisors besides , If you think about it carefully, it is not difficult to find that this kind of number can only have two prime factors at most , And these two prime factors are its two divisors respectively , And the third divisor, of course, is itself ！（ Of course, there may be special numbers such as the cubic of prime numbers , So we have to judge ） So our idea is , Find out first $n$ Is less than its divisor , Then we can use this divisor to work out $n$ Another divisor of , Then we'll decide whether it's a special number like the cubic of prime numbers （ Yes, you can output the answer directly ）, If not, we can judge whether these two divisors are different prime factors , Yes, you can output , It doesn't mean it's not $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;} // Take a look at this 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) // This question is too laggy { printf("%lld %lld %lld\n", p, q, n); continue; } if (!prime(p) || !prime(q) || p == q) cut; // Pay attention to the fact that two are equal printf("%lld %lld %lld\n", min(p, q), max(p, q), n); } return 0; } ``` ## Refer to https://blog.csdn.net/qq_39972971/article/details/82346390 https://ld246.com/article/1569220928268 https://zhuanlan.zhihu.com/p/267884783