一个很有意思的算法

例题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;
}