| Time Limit: 6000MS | Memory Limit: 65536K | |
| Total Submissions: 36448 | Accepted: 9737 | |
| Case Time Limit: 4000MS | ||
Description
Input
Output
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是素数，那么只可能是x-1能被p整除(此时x=1) 或 x+1能被p整除(此时x=p-1)。
如果我们经过测试得到a^(p-1) ≡1（mod p），如果p - 1为偶数，那么可以把p - 1表示成d × 2x这样的形式，我们就可以从ad×2^(x-1) = √(ap-1)开始，判断这个数mod p是不是等于1, 如果不是或者指数与2互质时，则判断这个数mod p是不是等于1或p-1，如果不是则无法通过测试(用上面二次探测定理的第二种描述来理解这个过程更好)。这种测试对于这个题目的数据范围已经毫无压力了。
判断素数的问题解决了，接下来要找一个合数的最小质因子。对于这么大的数据范围，我们应该采用一种叫做Pollard-rho的大数分解算法，然后找到一个数的最小质因子（实际上，这个算法理论上可以在O(N1/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 位的整数，我们只需要选取?? = 105个随机数而不是1010个。 不幸的是，这没有节省我们的开销，对于??( ?? = 105 )个人来说， 我们仍然要做?? 2 = 1010次比较和除法。但幸运的是，这里有一个更好的办法。
因为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 个数并两两进行比较，而是一个一个地生成并检查连续的两个数。
反复执行这个步骤并希望能够得到我们想要的数。
我们使用一个函数来生成伪随机数。
换句话说，我们不断地使用函数 ?? 来生成（看上去或者说感觉上像的）随机数。
并不是所有的函都能够这样做，但是有一个神奇的函数可以。它就是
??(??) = ( ??2 + ?? ) mod ??
( 我们可以自己指定??，也可以用 rand()随即函数来生成，这不是重点 )
我们从??1 = 2 或者其他数开始，让??2 = ??(??1), ??3 = ??(??2), … … 生成规则为??n+1 = ??(??n).
可以发现对于大部分的数据这个算法能够正常运行，但是对于某些数据，它将会进入无限循环。为什么呢？这是因为存在 ?? 环的原因。
那么，如何探测环的出现呢？
1. 一种方法是记录当前产生过的所有的数??1, ??2, … … ????,并检测是否存在???? = ????(?? < ??)。 在实际过程中，当 ?? 增长到一定大小时，可能会造成的内存不够用的情况。
2. 我们可以记录第一个随机数，从范围1开始，将当前数与记录的数比较是否相同，检测完一个范围没有发现相同则更新记录的数，并把范围扩大两倍，则要么找到因数退出循环，要么肯定会出现记录的数更新到循环节中，并且查找范围大于等于循环节长度的时候，这样便能检测出循环节退出循环，不会出现死循环。
下面为此题的AC代码+注释：
1 #include <iostream> 2 #include <cstdlib> 3 using namespace std; 4 typedef long long LL; 5 LL minfactor, p[11] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29}; 6 LL gcd(LL a, LL b) { 7 return b ? gcd(b, a % b) : a; 8 } 9 LL qmult(LL a, LL b, LL mod) { // 快速乘模 10 LL sum = 0; 11 while (b) { 12 if (b & 1) { 13 sum += a; 14 if (sum >= mod) sum -= mod; // 此处无需用%，%运算比减法慢很多 15 } 16 b >>= 1, a <<= 1; 17 if (a >= mod) a -= mod; 18 } 19 return sum; 20 } 21 LL qpow(LL a, LL b, LL mod) { // 快速幂模 22 LL res = 1; 23 while (b) { 24 if (b & 1) res = qmult(res, a, mod); 25 b >>= 1; 26 a = qmult(a, a, mod); 27 } 28 return res; 29 } 30 bool prime_test(LL n, LL a) { // 对整数n，底数a进行测试，返回true表示通过测试 31 LL p = qpow(a, n - 1, n); 32 if (p != 1) return false; 33 else { // 二次探测 34 LL s = n - 1; 35 while (!(s & 1) && p == 1) { 36 s >>= 1; 37 p = qpow(a, s, n); 38 } 39 if (p == 1 || p == n - 1) return true; 40 else return false; 41 } 42 } 43 bool Miller_Rabin(LL n) { // 对整数n进行Miller_Rabin素数测试，返回true表示通过测试 44 if (n <= 29) { // if这一块其实可以不用 45 for (int i = 0; i < 10; i++) { 46 if (n == p[i]) return true; 47 } 48 return false; 49 } 50 for (int i = 0; i < 10; i++) { // 利用前10个素数作为底数测试的正确率已经非常高 51 if (gcd(n, p[i]) == 1 && !prime_test(n, p[i])) return false; 52 } 53 return true; 54 } 55 LL randf(LL x, LL n, LL c) { // 满足要求的产生伪随机数函数 56 return (qmult(x, x, n) + c) % n; 57 } 58 LL pollard_rho(LL n, LL c) { // 查找n的因数，c为上面函数要用的随机数，c也可自己指定（但要有变化） 59 LL x = rand() % n, y = x, i = 1, k = 2, p; // 随机生成随机数的初始值，也可自己指定 60 while (true) { 61 i++; 62 x = randf(x, n, c); 63 p = gcd(y - x + n, n); 64 if (p > 1 && p < n) return p; 65 if (y == x) return n; // 判圈，返回n表示查找失败，要更新随机种子重新查找 66 if (i == k) { y = x; k <<= 1; } // 更新范围和记录的数 67 } 68 } 69 void find_factor(LL n) { // 查找所有因数 70 if (Miller_Rabin(n)) { 71 minfactor = min(minfactor, n); 72 return ; 73 } 74 LL p = n; 75 while (p == n) p = pollard_rho(n, rand() % (n - 1) + 1); // 查找失败则更新随机种子重新查找，直到找到因子 76 find_factor(p); // 递归查找更小因子 77 find_factor(n / p); 78 } 79 80 int main() { 81 int t; cin >> t; 82 while (t--) { 83 LL N; cin >> N; 84 if (Miller_Rabin(N)) cout << "Prime" << endl; 85 else { 86 minfactor = N; 87 find_factor(N); 88 cout << minfactor << endl; 89 } 90 } 91 return 0; 92 }
