因数分解 Pollard rho
算法思路
随机生成两个数 a , b a,b a,b,然后求 gcd ( n , a − b ) \gcd\pod{n,a-b} gcd(n,a−b),如果其值不为 1 1 1,则这个数就是 n n n的一个因数。
对于 a , b a,b a,b并不是随机生成的,而是按照一定的顺序生成的,一般情况下采用函数 F ( x ) = x 2 + C F(x)=x^2+C F(x)=x2+C,其中 C C C为随机生成的一个数。
通过不断地生成 a , b a,b a,b,直到找到 n n n的一个因数。
随机数生成
在使用rand()
函数的时候需要注意,它生成的数的范围时在
65536
65536
65536以内的,因此如果要获得更大的随机数,需要自行手写。
当通过 F ( x ) F(x) F(x)生成出一个相同的数之后,也就是说 F ( x ) F(x) F(x)成环了,此时就应该重新随机一个 C C C。
这个算法名称中的rho
就是代表着
F
(
x
)
F(x)
F(x)生成的数成环。
成环判断
如果是要记录随机生成的每个数,并用此来判断新生成的一个随机数有没有出现过显然时不可行的,因为这里随机生成的 a , b a,b a,b都非常大,很难高效地存储与判断。
不过根据它一定成环的性质,就可以用另一种方法来判断。
其实就是一个追及过程,追及的距离为从初始值到第一次出现重复随机数的次数。
具体如下:
每次令 a = F ( a ) , b = F ( F ( b ) ) a=F(a),b=F(F(b)) a=F(a),b=F(F(b)),即 a a a每次走一步, b b b每次走两步,比 b b b多走一步。
当 b b b追上 a a a的时候,即 b = a b=a b=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 , b a,b a,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年国家集训队论文,骆可强同学的论文。
最小公约数的优化
先要知道最小公约数的一个性质,对于一个正整数 b b b,则有 gcd ( a , n ) = gcd ( a × b , n ) \gcd(a,n)=\gcd(a\times b,n) gcd(a,n)=gcd(a×b,n)。
考虑一下每次求公因数的时间复杂度,也是 log \log log级别的,
如果频繁调用,其产生的时间开销也是无法接受的。
因此可以利用上面的性质,先累乘起来,再求 gcd \gcd gcd,减少调用的次数。
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;
}