算法思想

基本思路:
对于给定的一个整数 n,显然如果 n 为素数(Miller-Rabin 算法判断),那么算法结束,返回唯一素因子 n。

否则,pollard's rho 算 ***试着找到 n 的一个因子 a(a 并不一定是素数),然后递归 Pollard_rho( a ) 和 Pollard_rho( n/a ) ,即将 n 分为了两个部分,将原问题变成了规模更小的两个子问题。

如何求因子 a?
一般情况下可能会本能的想到枚举,但这样过于耗时。寻找因子是该算法中最重要的一步,该算法中采用随机化算法来查找因子 a。假设此时 n 仅有 2 个因子 p 和 q,那么如果我们用随机数来找 n 的因子,成功率为 1/n。接下来我们思路是如何提升成功率。
寻找一个因子 a,等价于寻找是否存在 k 个数,使得其中有 图片说明 ,由于生日悖论 我们可以得知当图片说明 时,该概率是 50%,提高了成功率。但不幸的是对于 10 位的整数,k = 1e5 时,仍要做图片说明 次比较,我们仍然选取 k 个数:图片说明 但我们不再询问是否存在图片说明 可以整除n,转而询问是否有图片说明 的情况。因为如果要问有多少个数可以整除n,那么只有p,q。但是如果问有多少个数满足图片说明,就会有更多的数的组合满足了,从而也增加了成功率。而图片说明 就是一个n的因子。但是经过验证,还是需要选择图片说明 才能有相当高的成功率。而图片说明 数量仍然很大,以至于可能不方便存在内存中。为了解决无法储存太多数的问题,Pollard's rho Algorithm 只将两个数存放在内存中,具体思路是:并不随机生成 k 个数并两两比较,而是一个一个地生成并检查连续的两个数。反复执行这个步骤并希望能够得到我们想要的数。

我们使用一个函数来生成伪随机数。
图片说明

使用的时候从x1=2开始,让x2 = f(X1), x3 = f(x2) ,xn+1=f(xn),于是顺着这个策略,不断判断gcd(xn-xn-1,n) 是否大于1即可。

利用这个函数因为x1和x2再调整时最终一定会出现循环,形成一个类似希腊字母rho的形状,故因此得名。于是问题变成了“如何判断环的出现”。一种方法是记录所有出现过的数,当然这会耗费大量内存,故舍弃;另一种方法是 Floyd 发明的算法,这里可以举一个有趣的例子说明“假设 A 和 B 在一个很长的圆形轨道上走,那么我们如何判断 B 是否走完一圈呢?我们可以让 B 的速度是 A 的二倍,他们同时出发,当 B 第一次追上 A,就知道 B 至少已经走了一圈”,同样的道理运用到该算法中,框架如下:

a = 2;
b = 2;
while(b != a){
    a = f(a);    //一倍速 
    b = f(f(b));//二倍速
    p = GCD(|b-a| , n);
    if(p > 1) return "Found factor: p"; 
}
return "Failed";

算法流程:

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

//****************************************************************
// Miller_Rabin 算法进行素数测试
//速度快,而且可以判断 <2^63的数
//****************************************************************

ll prime[10] = {2,3,5,7,11,13,17,19,23,29};
ll Quick_Multiply(ll a, ll b, ll c)
{
    ll ans = 0, res = a;
    while(b)
    {
        if(b & 1)
            ans = (ans + res) % c;
        res = (res + res) % c;
        b >>= 1;
    }
    return ans;
}
ll Quick_Power(ll a, ll b, ll m)
{
    ll ans = 1;
    while(b)
    {
        if(b & 1)
        {
            ans = ans * a % m;
        }
        a = a * a % m;
        b >>= 1;
    }
    return ans;
}
bool Miller_Rabin(ll x)    //判断素数
{
    ll i, j, k;
    ll s = 0, t = x-1;
    if(x == 2) return true; //2是素数
    if(x < 2 || !(x & 1)) return false;  //如果x是偶数或者是0,1,那它不是素数
    while(!(t & 1)) //将x-1分解成(2^s)*t的样子
    {
        s++;
        t >>= 1;
    }
    for(i = 0; i < 10 && prime[i] < x; i++) //随便选一个素数进行测试
    {
        ll a = prime[i];
        ll b = Quick_Power(a, t, x);   //先算出a^t
        for(j = 1; j <= s; j++) //然后进行s次平方
        {
            k = Quick_Multiply(b,b,x);  //求b的平方
            if(k == 1 && b != 1 && b != x-1)
                return false;
            b = k;
        }
        if(b != 1) return false;    //用费马小定理判断
    }
    return true;    //如果进行多次测试都是对的,那么x就很有可能是素数
}

//************************************************
//pollard_rho 算法进行质因数分解
//************************************************

ll factor[110];  //用来存放被分解的因子(这里也可以用map来存)
int tot = 0;        //因子个数
ll gcd(ll a, ll b)
{
    /* 返回a和b的最大公约数 */
    if(a == 0) return 1;
    if(a < 0) return gcd(-a, b);
    while(b) {
        ll t = a % b;
        a = b; b = t;
    }
    return a;
}

ll pollard_rho(ll x, ll c)
{
    /* 返回 x 的一个因子或 x 本身 */
    ll i = 1, k = 2;
    ll tx = rand()%x;
    ll y = tx;
    while(true)
    {
        i++;
        tx = (Quick_Multiply(tx, tx, x) + c) % x;
        ll d = gcd(y-tx, x);
        if(d != 1 && d != x) return d;  //x是个合数,所以因子应该大于1,小于n
        if(y == tx) return x;   //如果循环一圈了,仍未找到因子,返回到findFac函数重新取c
        if(i == k) y = tx, k += k;  //保持y是1倍速,tx是二倍速
    }
}

void findFac(ll n)
{
    /* 对 n 进行质因数分解 */
    if(Miller_Rabin(n))
    {
        factor[++tot] = n;
        return;
    }
    ll p = n;
    while(p >= n) p = pollard_rho(p, rand()%(n-1)+1);   //rand()%(n-1)+1 可以保证c值范围是在1~n-1之间
    findFac(p); //递归分解
    findFac(n/p);
}

int main()
{
    ll x, T;
    scanf("%lld",&T);
    while(T--)
    {
        memset(factor, 0, sizeof(factor));
        tot = 0;
        scanf("%lld",&x);
        tot = 0;
        if(Miller_Rabin(x)) {
            printf("isprime\n");
            printf("%lld\n",x);
        }
        else {
            printf("noprime\n");
            findFac(x);
            sort(factor+1, factor+tot+1);
            for(int i = 1; i < tot; i++) {
                if(factor[i] != factor[i+1])
                printf("%lld ",factor[i]);
            }
            printf("%lld\n",factor[tot]);
        }
    }
}

参考博客:
https://hacpai.com/article/1569220928268
https://www.cnblogs.com/fzl194/p/9047710.html