POJ 1811 -- Prime Test (Miller-Rabin素数测试 + Pollard-rho因子分解)

题目链接：http://poj.org/problem?id=1811

Prime Test

Time Limit: 6000MS		Memory Limit: 65536K
Total Submissions: 36448		Accepted: 9737
Case Time Limit: 4000MS

Description

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⁵⁴).

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个整数N（2≤N≤254），判断它们是不是素数，如果不是，输出它的最小质因子。

解析：通过Miller-Rabin素数测试判断是否素数，通过Pollard-rho因子分解算法得到最小因子。

费马小定理：假如p是质数，且GCD(a,p)=1，那么 a^(p-1) ≡1（mod p）（假如p是质数，且a,p互质，那么 a的(p-1)次方除以p的余数恒等于1）。
用费马小定理的逆命题可以来判定素数，但是其逆命题不一定成立，即费马小定理是一个数为素数的必要条件，不是充分条件，所以有一定的概率出错，需要多次测试来减小出错的概率。
每次测试随机取一个正整数a,要保证 GCD(a,p)=1; 用快速幂来计算a^(p-1)mod p看是否为 1，如果为1则通过这次测试。
但是这样也可能出错，因为存在称做Carmichael数的合数，使得对这些合数进行多次测试都能通过，但是它们不是素数。前3个Carmichael数是561,1105,1729。Carmichael数是非常少的。在1~100000000范围内的整数中，只有255个Carmichael数。利用二次探测定理可以对上面的素数判定算法作进一步改进，以避免将Carmichael数当作素数。实际上就是添加了必要条件，这样我们就对Miller_Rabin素数测试进行了强化，称为改进Miller_Rabin素数测试。注意改进后的算法仍然是不确定算法，但非常实用。通常认为，Miller_Rabin素数测试的正确率可以令人接受，随机选取k个底数进行测试算法的失误率大概为4^(-k)。

二次探测定理：如果p是一个素数，且0<x<p,则方程x^2≡1(mod p)的解必为 x=1或p-1。
或者：如果p是一个素数，则x^2≡1(mod p)的解满足 x≡±1(mod p)。

这是显然的，因为相当于p能整除，也即p能整除(x+1)(x-1)。

 

由于p是素数，那么只可能是x-1能被p整除(此时x=1) 或 x+1能被p整除(此时x=p-1)。

 

 

 

 

如果我们经过测试得到a^(p-1) ≡1（mod p），如果p - 1为偶数，那么可以把p - 1表示成d × 2^x这样的形式，我们就可以从a^{d×2^(x-1)} = √(a^p-1)开始，判断这个数mod p是不是等于1, 如果不是或者指数与2互质时，则判断这个数mod p是不是等于1或p-1，如果不是则无法通过测试(用上面二次探测定理的第二种描述来理解这个过程更好)。这种测试对于这个题目的数据范围已经毫无压力了。

 

 

 

 

判断素数的问题解决了，接下来要找一个合数的最小质因子。对于这么大的数据范围，我们应该采用一种叫做Pollard-rho的大数分解算法，然后找到一个数的最小质因子（实际上，这个算法理论上可以在O(N^1/4)间内找到一个数的所有质因子）。

 

 

Pollard-rho算法：

问题的设置：

我们假设 𝑵 是一个能被分解为 𝒑 ∗ 𝒒 的数 ( 𝑵 = 𝒑 ∗ 𝒒 且 𝒑 ≠ 𝒒 ) ，我们的目标是找到 𝑵 的其中一个因子 𝒑 或 𝒒 (另外一个因子可以通过除 𝑁 来得到 )。

 

现在，我们随机从[2,N]中选出n个数，刚好刷出p或q的概率很小，但是如果它们两两作差呢，差为p或q的概率会如何？

如果我们在[1,1000]中选取 𝑘 个数 𝑥1,… … , 𝑥𝑘， 取得的𝑘个数中满足𝑥𝑖 − 𝑥𝑗 = 42 的概率是多少呢？这个概率和 𝑘 有什么关系呢？

 

大约在 𝑘 = 30 的时候，成功的概率已经超过一半。 我们有一个大区间[1,1000]，但是仅仅生成约 30 个随机数就让我们成功的概率超过了一半， 大约在 𝑘 = 100 的时候，我们几乎可以确保我们是成功的。这是一个非常重要的观察， 它被成为 生日问题 (𝑏𝑖𝑟𝑡ℎ𝑑𝑎𝑦 𝑝𝑟𝑜𝑏𝑙𝑒𝑚) 或者 生日悖论 (𝑏𝑖𝑟𝑡ℎ𝑑𝑎𝑦 𝑝𝑎𝑟𝑎𝑑𝑜𝑥)。

(注：当k=sqrt(N)时，概率约为50%。因为sqrt(N)个数任意组队，有N/2种方案。)

我们已经知道对于𝑘 = √𝑁, 我们可以将可能性提高到50%， 所以，可以粗略地说，我们可以将可能性从1/𝑁 提高到 √ (1/𝑁)。

因此，对于一个 10 位的整数，我们只需要选取𝑘 = 10⁵个随机数而不是10¹⁰个。 不幸的是，这没有节省我们的开销，对于𝑘( 𝑘 = 10⁵ )个人来说， 我们仍然要做𝑘 2 = 10¹⁰次比较和除法。但幸运的是，这里有一个更好的办法。 

 

 

因为N=p*q，p和q都很稀有，但是p和q的倍数便是一波大军。

我们可以选取 𝑘 个数 𝑥1, … … , 𝑥𝑘，不再问是否存在𝑥𝑖 − 𝑥𝑗能够整除𝑁， 转而询问是否存在 𝑔𝑐𝑑(|𝑥_i − 𝑥_j| , 𝑁) > 1的情况。若存在，则𝑔𝑐𝑑(|𝑥𝑖 − 𝑥j| , 𝑁) 为N的一个因子，通过递归找出更小的因子，直到因子为素数。

 

所以，一个简单的策略如下：

1.  在区间[𝟐, 𝑵 − 𝟏]中随即选取 𝒌 个数，𝒙𝟏, … … , 𝒙𝒌                                                                                                                                                                  

 

 

2.  判断是否存在𝒈𝒄𝒅(|𝒙_i − 𝒙_j| , 𝑵) > 𝟏, 若存在，𝒈𝒄𝒅(|𝒙_i − 𝒙_j| ,𝑵) 是𝑵的一个因子

 

但是很早就出现了一个问题，我们需要选取大约 𝑁^1/4 个数， 这个数量太大了，以至于我们并不能将其存放在内存中 .

Pollard′s rho algorithm 解决了数太多无法储存的问题，只将两个数存放在内存中。

我们并不随机生成 k 个数并两两进行比较，而是一个一个地生成并检查连续的两个数。

反复执行这个步骤并希望能够得到我们想要的数。

 

 

我们使用一个函数来生成伪随机数。

换句话说，我们不断地使用函数 𝑓 来生成（看上去或者说感觉上像的）随机数。

并不是所有的函都能够这样做，但是有一个神奇的函数可以。它就是

𝑓(𝑥) = ( 𝑥² + 𝑎 ) mod 𝑁

 

( 我们可以自己指定𝑎，也可以用 rand()随即函数来生成，这不是重点 )

 

我们从𝑥₁ = 2 或者其他数开始，让𝑥₂ = 𝑓(𝑥₁), 𝑥₃ = 𝑓(𝑥₂), … … 生成规则为𝑥_n+1 = 𝑓(𝑥_n).

 

可以发现对于大部分的数据这个算法能够正常运行，但是对于某些数据，它将会进入无限循环。为什么呢？这是因为存在 𝑓 环的原因。

 

 

那么，如何探测环的出现呢？

1. 一种方法是记录当前产生过的所有的数𝑥1, 𝑥2, … … 𝑥𝑛,并检测是否存在𝑥𝑙 = 𝑥𝑛(𝑙 < 𝑛)。 在实际过程中，当 𝑛 增长到一定大小时，可能会造成的内存不够用的情况。

 

2. 我们可以记录第一个随机数，从范围1开始，将当前数与记录的数比较是否相同，检测完一个范围没有发现相同则更新记录的数，并把范围扩大两倍，则要么找到因数退出循环，要么肯定会出现记录的数更新到循环节中，并且查找范围大于等于循环节长度的时候，这样便能检测出循环节退出循环，不会出现死循环。

 

下面为此题的AC代码+注释：

 

 

#include <iostream>
#include <cstdlib>
using namespace std;
typedef long long LL;
LL minfactor, p[11] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
LL gcd(LL a, LL b) {
    return b ? gcd(b, a % b) : a;
}
LL qmult(LL a, LL b, LL mod) {  // 快速乘模
    LL sum = 0;
    while (b) {
        if (b & 1) {
            sum += a;
            if (sum >= mod) sum -= mod;  // 此处无需用%，%运算比减法慢很多
        }
        b >>= 1, a <<= 1;
        if (a >= mod) a -= mod;
    }
    return sum;
}
LL qpow(LL a, LL b, LL mod) {  // 快速幂模
    LL res = 1;
    while (b) {
        if (b & 1) res = qmult(res, a, mod);
        b >>= 1;
        a = qmult(a, a, mod);
    }
    return res;
}
bool prime_test(LL n, LL a) {  // 对整数n，底数a进行测试，返回true表示通过测试
    LL p = qpow(a, n - 1, n);
    if (p != 1) return false;
    else {  // 二次探测
        LL s = n - 1;
        while (!(s & 1) && p == 1) {
            s >>= 1;
            p = qpow(a, s, n);
        }
        if (p == 1 || p == n - 1) return true;
        else return false;
    }
}
bool Miller_Rabin(LL n) {  // 对整数n进行Miller_Rabin素数测试，返回true表示通过测试
    if (n <= 29) {    // if这一块其实可以不用
        for (int i = 0; i < 10; i++) {
            if (n == p[i]) return true;
        }
        return false;
    }
    for (int i = 0; i < 10; i++) {  // 利用前10个素数作为底数测试的正确率已经非常高
        if (gcd(n, p[i]) == 1 && !prime_test(n, p[i])) return false;
    }
    return true;
}
LL randf(LL x, LL n, LL c) {  // 满足要求的产生伪随机数函数
    return (qmult(x, x, n) + c) % n;
}
LL pollard_rho(LL n, LL c) {  // 查找n的因数，c为上面函数要用的随机数，c也可自己指定（但要有变化）
    LL x = rand() % n, y = x, i = 1, k = 2, p;  // 随机生成随机数的初始值，也可自己指定
    while (true) {
        i++;
        x = randf(x, n, c);
        p = gcd(y - x + n, n);
        if (p > 1 && p < n) return p;
        if (y == x) return n;   // 判圈，返回n表示查找失败，要更新随机种子重新查找
        if (i == k) { y = x; k <<= 1; }  // 更新范围和记录的数
    }
}
void find_factor(LL n) {  // 查找所有因数
    if (Miller_Rabin(n)) {
        minfactor = min(minfactor, n);
        return ;
    }
    LL p = n;
    while (p == n) p = pollard_rho(n, rand() % (n - 1) + 1);  // 查找失败则更新随机种子重新查找，直到找到因子
    find_factor(p);     // 递归查找更小因子
    find_factor(n / p);
}

int main() {
    int t; cin >> t;
    while (t--) {
        LL N; cin >> N;
        if (Miller_Rabin(N)) cout << "Prime" << endl;
        else {
            minfactor = N;
            find_factor(N);
            cout << minfactor << endl;
        }
    }
    return 0;
}