因数分解 Pollard rho

算法思路

随机生成两个数a,ba,ba,b，然后求gcd⁡(n,a−b)\gcd\pod{n,a-b}gcd(n,a−b)，如果其值不为111，则这个数就是nnn的一个因数。

对于a,ba,ba,b并不是随机生成的，而是按照一定的顺序生成的，一般情况下采用函数F(x)=x2+CF(x)=x^2+CF(x)=x2+C，其中CCC为随机生成的一个数。

通过不断地生成a,ba,ba,b，直到找到nnn的一个因数。

随机数生成

在使用rand()函数的时候需要注意，它生成的数的范围时在655366553665536以内的，因此如果要获得更大的随机数，需要自行手写。

当通过F(x)F(x)F(x)生成出一个相同的数之后，也就是说F(x)F(x)F(x)成环了，此时就应该重新随机一个CCC。

这个算法名称中的rho就是代表着F(x)F(x)F(x)生成的数成环。

成环判断

如果是要记录随机生成的每个数，并用此来判断新生成的一个随机数有没有出现过显然时不可行的，因为这里随机生成的a,ba,ba,b都非常大，很难高效地存储与判断。

不过根据它一定成环的性质，就可以用另一种方法来判断。

其实就是一个追及过程，追及的距离为从初始值到第一次出现重复随机数的次数。

具体如下：

每次令a=F(a)，b=F(F(b))a=F(a)，b=F(F(b))a=F(a)，b=F(F(b))，即aaa每次走一步，bbb每次走两步，比bbb多走一步。

当bbb追上aaa的时候，即b=ab=ab=a的时候，就可以判断成环。

一些优化技巧

乘法的优化

对于两long long的数相乘，并取模，貌似没有什么有效的解决方案，直接乘会溢出，而且也没有可用的更大的整数类型。

一种比较显然的做法，利用快速幂的思路，就可以写出快速乘

ll mul (ll a , ll b , ll p) //a*b%p
{ll s = 0;for ( ; b ; ){if (b & 1) s = (s + a) % p;a = (a + a) % p;b = b >> 1;}return s;
}

不难看出来这个的时间复杂度是O(log⁡b)O(\log b)O(logb)的，而且这里的a,ba,ba,b都非常大，于是每次乘法的需要循环的次数估计要有四五十。

再回过头来看一下Pollard rho算法，它需要计算乘法的地方是非常多的，于是这样的时间开销是无法接受的。

基于底层的计算原理，其实还有一种更快的乘法，称之为O(1)快速乘。

ll mul (ll a , ll b , ll p) //a*b%p
{ll s = a * b - (ll)((long double) a * b / p +0.5) * p;return s < 0 ? s + p : s;
}

两个作差，要么都溢出，要么都不溢出，最终的结果都是一样了。

这个也可以参见2009年国家集训队论文，骆可强同学的论文。

最小公约数的优化

先要知道最小公约数的一个性质，对于一个正整数bbb，则有gcd⁡(a,n)=gcd⁡(a×b,n)\gcd(a,n)=\gcd(a\times b,n)gcd(a,n)=gcd(a×b,n)。

考虑一下每次求公因数的时间复杂度，也是log⁡\loglog级别的，

如果频繁调用，其产生的时间开销也是无法接受的。

因此可以利用上面的性质，先累乘起来，再求gcd⁡\gcdgcd，减少调用的次数。

Code

//#pragma GCC optimize (2)
//#pragma G++ optimize (2)
#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <iostream>
#include <vector>
#include <queue>
#define G getchar
#define ll long long
using namespace std;ll read()
{char ch;for(ch = G();(ch < '0' || ch > '9') && ch != '-';ch = G());ll n = 0 , w;if (ch == '-'){w = -1;ch = G();} else w = 1;for(;'0' <= ch && ch <= '9';ch = G())n = (n<<1)+(n<<3)+ch-48;return n * w;
}const int N = 100003;
ll n , ans;
ll pri[N] , k[N] , m;
int ss[16] = {2 , 3 , 5 , 7 , 11 ,13 , 17 , 19 , 23 , 29 , 325 , 9375 , 28187 , 450775 , 9870504 , 1795265022};ll mul (ll a , ll b , ll p) //a*b%p
{/*ll s = 0;for ( ; b ; ){if (b & 1) s = (s + a) % p;a = (a + a) % p;b = b >> 1;}return s;*/ll s = a * b - (ll)((long double) a * b / p +0.5) * p;return s < 0 ? s + p : s;
}ll ksm (ll x , ll y , ll p)
{ll s = 1;for ( ; y ; ){if (y & 1) s = mul(s , x , p);x = mul (x , x , p);y = y >> 1;}return s;
}bool MR_detect (ll n , ll a)
{if (n == a) return 1;if (a % n == 0) return 1;if (ksm(a , n - 1 , n) != 1) return 0;ll p = n - 1 , lst = 1;for ( ; ((p & 1) ^ 1) && lst == 1 ; ){p = p >> 1;lst = ksm(a , p , n);if (lst != 1 && lst != n - 1) return 0;}return 1;
}bool MR(ll  n)
{if (n < 2) return 0;for (int i = 0 ; i < 16 ; i ++)if (! MR_detect(n , ss[i])) return 0;return 1;
}ll F(ll x , ll C , ll p)
{return (mul(x , x , p) + C) % p;
}ll gcd(ll a , ll b)
{return a % b ? gcd(b , a % b) : b;
}ll Rand()
{return (ll)rand() + ((ll)rand() << 15) + ((ll)rand() << 30) + ((ll)rand() << 45);
}ll PR(ll n)
{if (MR(n)) return n;if (n == 4) return 2;for ( ; ; ){ll C = Rand() % (n - 1) + 1;ll p1 = 0 , p2 = 0;ll s = 1 , tmp;for ( ; ; ){for (int i = 0 ; i < 128 ; i++){p2 = F(F(p2 , C , n) , C , n);p1 = F(p1 , C , n);tmp = mul(s , abs(p1 - p2) , n);if (tmp == 0 || p1 == p2) break;s = tmp;}tmp = gcd(n , s);if (tmp > 1) return tmp;if (p1 == p2) break;}}
}int main()
{//freopen("1.txt","r",stdin);//freopen("2.txt","w",stdout);srand((unsigned)0);//printf("350\n");//for (int i = 1 ; i <= 350 ; i++)printf("%lld\n", Rand());for (int T = read() ; T ;T--){n = read();if (MR(n)) printf("Prime\n"); else{m = 0;for ( ; n != 1 ; ){ll tmp = PR(n);for ( ; !MR(tmp) ; )tmp = PR(tmp);m++;pri[m] = tmp;k[m] = 0;for ( ; n % tmp == 0 ; ){k[m]++;n = n / tmp;}}//for (int i = 1 ; i <= m ; i++)//    printf("%lld^%lld\n", pri[i] , k[i]);ans = 0;for (int i = 1 ; i <= m ; i++)ans = max(ans , pri[i]);printf("%lld\n", ans);}}   return 0;
}