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

分手后的思念是犯贱 2022-04-24 07:48 170阅读 0赞

# 前言 #

这是一个伤心的故事：*我已经好久没有写博客了*。。。  
当然这次跟大家分享的是一个很玄学的东西——**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 1 1 <= T <=  20 20 20 ), then the following T lines each contains an integer number N ( 2 2 2 <= N <  2 54 2^\{54\} 254).

#### 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类型。  
因此，我们就必须引入今天的 **玄学主题** 了： **P o l l a r d r h o Pollard rho Pollardrho**.

### *Pollardrho* ###

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

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

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

当然，其实也不是一直循环，是在满足一个条件后就可以退出了:  
![在这里插入图片描述][aHR0cHM6Ly9pLmxvbGkubmV0LzIwMTkvMDQvMjQvNWNiZmY4YWM5MTgzYi5qcGc]  
因为在大量数据中表明，x的变化会满足上图，也就是希腊字母  ρ \\rho ρ。  
因此在那个圆圈的时候，就说明已经出现循环了，因此此次随机的值也就没有任何意义了，就可以退出了。

### 注意事项 ###

1.  因为这个数据十分的巨大，要用  l o n g l o n g longlong longlong 来存。 当然用 l o n g l o n g longlong longlong也是有可能爆空间的，因此要**边更新边模**。
2.  当然，在x变化的地方，也是要用到**快速乘**（类似于快速幂）以防爆空间的。 在 m i l l e r r a b i n miller\_rabin millerrabin判素数中会运用到**小费马定理**和**二次探测定理**，在这之中也要用到快速幂来计算，否则依然会爆空间。
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);
        }
    }

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MzkyNTIzOQ_size_16_color_FFFFFF_t_70]: /images/20220217/fcff4bdfbbd94a488c948605e358cfdc.png
[Link 1]: https://baike.baidu.com/item/%E7%94%9F%E6%97%A5%E6%82%96%E8%AE%BA/2715290?fr=aladdin
[aHR0cHM6Ly9pLmxvbGkubmV0LzIwMTkvMDQvMjQvNWNiZmY4YWM5MTgzYi5qcGc]: /images/20220217/62c1d3f16a914414a92e66ba0aecbd18.png