[玄学]——数论高级之分解质因数（Pollard_rho）（POJ 1811）

前言

这是一个伤心的故事：我已经好久没有写博客了。。。
当然这次跟大家分享的是一个很玄学的东西——Pollard_rho
为什么说他很玄学呢？大家可以在等下的代码中看出来（全是rand。。。）
当然，这也是我的第一篇用markdown写的博客哈~~ 值得庆祝值得庆祝

正文

题目

全英版

Given a big integer number, you are required to find out whether it’s a prime number.

Input

The first line contains the number of test cases T ($1$ <= T <= $20$ ), then the following T lines each contains an integer number N ($2$ <= N < $2^{54}$).

Output

For each test case, if N is a prime number, output a line containing the word “Prime”, otherwise, output a line containing the smallest prime factor of N.

Sample Input

2
5
10

Sample Output

Prime
2

中翻版

（什么翻译啊这个）
我觉得这个翻译什么的，能自己翻还是自己动手吧。。。
其实题目描述差不多就是这样的：
给定t组，如果说输入的是一个质数，就输出prime，否则输出最小的质数。

解析

看起来这个题目似乎很简单啊，典型的 暴力打表 嘛，但是如果说这个数据很大怎么办？？关键还不只是数据大而已，是非常大诶——longlong类型。
因此，我们就必须引入今天的 玄学主题 了：$Pollardrho$.

Pollardrho

Pollard_rho是一种基于随机的算法，它的思路是先用 miller_rabin 来判断当前数是否已经是素数了，如果是的话记录并返回。如果不是，我们设要分解的数为n,那么我们考虑去找一个当前数的因数p，找到之后再分别对p和n/p分解质因数，有一点类似分治但是不是分治。
那么，我们要如何快速地找到当前数的因数呢？
这里，我们采用以下方法：先随机生成两个1~n的数，利用这两个数差的绝对值和n比较，来判断这两个数差的绝对值是不是n的一个因数。

具体如图：

在这里只到所有的p都为质数为止。
可是为什么要采用这样的方法来找n的因数呢？那么这里就要引入一个东西了——生日悖论
当然在这里就不再赘述了哈，大家可以自己看下。
其实简单来说，就是选一个概率更高的 “算法” 来进行推因数而已。因为的话，就算一个数的因子再多，你也不一定很快就能随机一个出来对吧 （毕竟RP摆在那里的呢）
但是如果说是选两个数的话，就比较容易了。设p为n的一个因数，那么这两个数就是x和x+p而已，这样随机的可能性就大得多。

当然，这两个数也不能随便随机吧，因此又来了一个玄学的方法了：
设这两个数为x和y，x是随机的一个数，令y 等于 x；
设一个随机数c，令 $x=x\ast x+c(\mathrm{m}\mathrm{o}\mathrm{d}n)$ ；
这个时候就能进行减法了；
然而每进行 $2^k$ 次运算就要重新将y赋为当前的x；
“一直循环”，直到找到n的因数为止

当然，其实也不是一直循环，是在满足一个条件后就可以退出了:

因为在大量数据中表明，x的变化会满足上图，也就是希腊字母 $\rho$。
因此在那个圆圈的时候，就说明已经出现循环了，因此此次随机的值也就没有任何意义了，就可以退出了。

注意事项

1. 因为这个数据十分的巨大，要用 $longlong$ 来存。 当然用$longlong$也是有可能爆空间的，因此要边更新边模。

2. 当然，在x变化的地方，也是要用到快速乘（类似于快速幂）以防爆空间的。 在$mille{r}_{r}abin$判素数中会运用到小费马定理和二次探测定理，在这之中也要用到快速幂来计算，否则依然会爆空间。

3. 最后一点也十分重要——wa过几次的血泪教训——在快速幂中的乘法也要用快速乘！！！！！

参考代码

#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <ctime>
using namespace std;
#define LL long long
#define min(a, b) a < b ? a : b
#define max(a, b) a > b ? a : b

template <typename T>
void re (T &x){
    x = 0; char c = getchar ();
    while (c < '0' || c > '9') c = getchar ();
    while (c >= '0' && c <= '9'){
        x = (x << 1) + (x << 3) + c - 48;
        c = getchar ();
    }
}

template <typename T>
void pr (T x){
    if (x < 0){
        putchar ('-');
        x = ~x + 1;
    }
    if (x / 10) pr (x / 10);
    putchar (x % 10 + 48);
}

int t;
LL n, minn = 0x3f3f3f3f;

LL Gcd (LL x, LL y){
    if (x == 0) return 1;
    if (x < 0) return Gcd (-x, y);
    /*if (x < y){
        x ^= y;
        y ^= x;
        x ^= y;
    }*/
    return y == 0 ? x : Gcd (y, x % y);
}

LL qkmul (LL x, LL b, LL mod){
    LL tot = 0;
    while (b){
        if (b & 1) tot = tot + x % mod;
        x = (x * 2) % mod;
        b >>= 1;
    }
    return tot % mod;
}

LL qkpow (LL x, LL b, LL mod){
    LL tot = 1;
    while (b){
        if (b & 1) tot = qkmul (tot, x, mod);
        x = qkmul (x, x, mod);
        b >>= 1;
    }
    return tot % mod;
}

bool miller_rabin (LL n){//判素数，利用二次探测和费马小定理，双重判断
    if (n == 2) return 1;
    if (n < 2 || n % 2 == 0)
        return 0;
    LL u = n - 1, pre, x;
    int k = 0;
    while ((u & 1) == 0){//可以将a^u分为(a^v)^(2^k)
        u >>= 1;
        k ++;
    }
    for (int i = 1; i <= 10; i++){
        x = rand () % (n - 1) + 1;
        x = qkpow (x, u, n);
        pre = x;
        for (int j = 1; j <= k; j++){
            x = qkmul (x, x, n);
            if (x == 1 && pre != 1 && pre != n - 1) return 0;//二次探测
            pre = x;
        }
        if (x != 1) return 0;//费马小定理
    }
    return 1;
}

//先随机生成两个数x和c，每次都使x=x*x+c(mod n) ，另外再设置一个数y，y的初始值为x，每进行2^k次操作，就让y记为现在的x
LL pollard_rho (LL n, LL c){//随机生成两个1~n的数，利用这两个数差的绝对值和n比较，来判断这两个数差的绝对值是不是n的因数
    LL i = 1, k = 2;
    LL x = rand () % n, y = x;
    while (1){
        x = (qkmul (x, x, n) + c) % n;
        LL gcd = Gcd (x - y, n);
        if (gcd > 1 && gcd < n)
            return gcd;
        if (x == y) return n;//如果说x=y了就说明随机失败，返回n重新来过
        if (i == k){
            k <<= 1;
            y = x;
        }
        i ++;
    }
}

void Find (LL n){//找出n的因数，并一直分下去，分成p和n/p，直到所有都变成素数，返回
    if (miller_rabin (n)){
        minn = min (minn, n);
        return ;
    }
    LL p = n;
    while (p == n)//理论上来说不可能大于n
        p = pollard_rho (p, rand () % (n - 1) + 1);
    Find (p);
    Find (n / p);
}

int main (){
    srand (time (0));
    re (t);
    while (t --){
        re (n);
        if (miller_rabin (n)){//素数测试
            printf ("Prime\n");
            continue;
        }
        minn = 0x3f3f3f3f;
        Find (n);
        pr (minn);
        putchar (10);
    }
}