Pollard-Rho算法
作者：互联网
一个很有意思的算法
要求的是一个数的最大质因子，常规做法肯定是过不了的这个数据范围的
MillerRabin概率测素数法
根据两个定理
费马小定理:若\(p\)为质数,\(a^{p-1}\equiv 1 \pmod{p}\)
二次检查定理:若\(p\)为质数,\(x^2 \equiv 0 \pmod{p}\)的解为\(x_1=1,x_2=p-1\)
有了这两个定理,就可以开始MillerRabin算法了
通过费马小定理可知,如果一个数满足费马小定理,那么这个数很可能为素数,并不是一定.
对于一个素数 \(p , p-1\) 一定是偶数,那么不妨设 \(p-1=m \times 2^q\)
那么我们就可以通过一个数 \(a^m\) 不断自乘得到下面这个序列
\(a^{m \times 2},a^{m \times 2^2},...,a^{m \times 2^q}\)
若该序列的数全部通过二次检查,那么这个数就很大概率是素数了
通过前人智慧,我们知道每通过一次检查,那么\(p\)不是素数的概率为\(\frac{1}{4}\),
那么经过这一系列的检查,\(p\)不是素数的概率就近乎为0了
通过MillerRabin算法,我们就可以O(1)判断一个数是不是素数
Pollard-Rho算法
接下来才是重头戏
我们不妨考虑构造一个伪随机数序列,然后取相邻的两项来求gcd.为了生成一串优秀的随机数,Pollard设计了这样一个函数:
\(f(x)=(x^2+c)\mod N\)其中c是一个随机的常数.
我们随便取一个\(x_1\),令\(x_2=f(x_1),x_3=f(x_2),...,x_i=f(x_{i-1})\).
在一定的范围内,这个数列是基本随机的,可以取相邻两项作差求gcd.
不过之所以叫伪随机数,是因为这个数列里面会包含有"死循环".
{1,8,11,8,11,8,11,...}这个数列是\(c=7,N=20,x_1=1\)时得到的随机数列.
这个数列会最终在8和11之间不断循环.循环节的产生不难理解:在模N意义下,整个函数的值域只能是{0,1,2,...,N-1}
总有一天,函数在迭代的过程中会带入之前的值,得到相同的结果.
生成数列的轨迹很像一个希腊字母\(\rho\),所以整个算法的名字叫做Pollard-Rho.
为了判断环的存在,可以用一个简单的Floyd判圈算法,也就是"龟兔赛跑".
假设乌龟为\(t\),兔子为\(r\),初始时\(t=r=1\).假设兔子的速度是乌龟的一倍.
过了时间\(i\)后,\(t=i,r=2i\).此时两者得到的数列值\(x_t=x_i,x_r=x_{2i}\).
假设环的长度为\(c\),在环内恒有:\(x_i=x_{i+c}\).
如果龟兔"相遇",此时有:\(x_r=x_t\),也就是\(x_i=x_{2i}=x_{i+kc}\).
此时两者路径之差正好是环长度的整数倍。
这样以来,我们得到了一套基于Floyd判圈算法的Pollard Rho 算法.
int f(int x,int c,int n)
{
return (x*x+c)%n;
}
int pollard_rho(int N)
{
int c = rand() % (N - 1) + 1;
int t=f(0, c, N),r=f(f(0, c, N), c, N);//两倍速跑
while(t != r)
{
int d = gcd(abs(t-r), N);
if (d > 1)
return d;
t=f(t, c, N),r=f(f(r, c, N), c, N);
}
return N;//没有找到,重新调整参数c
}
由于求\(gcd\)的时间复杂度是\(O(\log{N})\)的,频繁地调用这个函数会导致算法变得异常慢.
我们可以通过乘法累积来减少求\(gcd\)的次数.
显然的,如果\(\gcd(a,b)>1\),则有\(\gcd(ac,b)>1\),c是某个正整数.
更近一步的,由欧几里得算法,我们有\(gcd(a,b) = gcd(ac\mod b,b) >1\).
我们可以把所有测试样本\(|t-r|\)在模\(N\)意义下乘起来,再做一次\(gcd\).
选取适当的相乘个数很重要.
我们不妨考虑倍增的思想:每次在路径\(\rho\)上倍增地取一段\([2^{k-1},2^k]\)的区间.
将\(2^{k-1}\)记为\(l\),\(2^k\)记为\(r\).
取而代之的,我们每次取的\(gcd\)测试样本为\(|x_i-x_l|\),其中\(l \le i \le r\).
我们每次积累的样本个数就是\(2^{k-1}\)个,是倍增的了.
这样可以在一定程度上避免在环上停留过久,或者取gcd的次数过繁的弊端.
当然，如果倍增次数过多，算法需要积累的样本就越多。我们可以规定一个倍增的上界。
\(127=2^7-1\)，据推测应该是限制了倍增的上界。
一旦样本的长度超过了127，我们就结束这次积累，并求一次\(\gcd\)。
127这个数应该是有由来的。
128似乎作为“不超过某个数的质数个数”出现在时间复杂度中。
不过你也可以理解为，将"迭代7次"作为上界是实验得出的较优方案。
这样,我们就得到了一套完整的,基于路径倍增的Pollard Rho 算法.
ll Pollard_Rho(ll x)
{
ll s = 0, t = 0, c = (ll)rand() % (x - 1) + 1;
int step = 0, goal = 1;
ll val = 1;
for (goal = 1; ; goal <<= 1, s = t, val = 1)
{
for (step = 1; step <= goal; ++step)
{
t = f(t, c, x);
val = (lll)val * abs(t - s) % x;
if (step % 127 == 0)
{
ll d = gcd(val, x);
if (d > 1)
return d;
}
}
ll d = gcd(val, x);
if (d > 1)
return d;
}
}
完整代码
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<ctime>
#define lll __int128
using namespace std;
typedef long long ll;
ll prime[8] = {2, 3, 5, 7, 11, 13, 17, 19};
ll max_factor;
ll qmi(ll a, ll b, ll c)
{
ll res = 1, base = a;
while (b)
{
if (b & 1)
res = (lll)res * base % c;
base = (lll)base * base % c;
b >>= 1;
}
return res;
}
bool MillerRabin(ll n)
{
if (n == 2)
return true;
if (n < 2 || !(n & 1))
return false;
ll t = n - 1, k = 0;
while (!(t & 1))
k++, t >>= 1;
for (int i = 0; i < 8; ++i)
{
if (n == prime[i])
return true;
ll a = qmi(prime[i], t, n), temp = a;
for (int j = 1; j <= k; ++j)
{
temp = (lll)a * a % n;
if (temp == 1 && a != 1 && a != n - 1)
return false;
a = temp;
}
if (a != 1)
return false;
}
return true;
}
ll gcd(ll a, ll b)
{
if (b == 0)
return a;
return gcd(b, a % b);
}
ll f(ll x, ll c, ll n)
{
return ((lll)x * x + c) % n;
}
ll Pollard_Rho(ll x)
{
ll s = 0, t = 0, c = (ll)rand() % (x - 1) + 1;
int step = 0, goal = 1;
ll val = 1;
for (goal = 1; ; goal <<= 1, s = t, val = 1)
{
for (step = 1; step <= goal; ++step)
{
t = f(t, c, x);
val = (lll)val * abs(t - s) % x;
if (step % 127 == 0)
{
ll d = gcd(val, x);
if (d > 1)
return d;
}
}
ll d = gcd(val, x);
if (d > 1)
return d;
}
}
void find(ll x)
{
if (x <= max_factor || x < 2)
return;
if (MillerRabin(x))
{
max_factor = max_factor > x ? max_factor : x;
return;
}
ll p = x;
while (p >= x)
p = Pollard_Rho(x);
while (x % p == 0)
x /= p;
find(x), find(p);
}
int main()
{
int T;
scanf("%d", &T);
while (T--)
{
srand((unsigned)time(0));
ll n;
scanf("%lld", &n);
max_factor = 0;
find(n);
if (max_factor == n)
printf("Prime\n");
else
printf("%lld\n", max_factor);
}
return 0;
}
标签：return,gcd,int,ll,Pollard,算法,Rho 来源： https://www.cnblogs.com/A2484337545/p/14436317.html