质数进阶
质数进阶#
质数（英语：Prime number），又称素数，指在大于1的自然大于1的自然数中，除了1和此整数自
身外，无法被其他自然数整除的数（也可定义为只有1和本身两个因数的数）。
暴力枚举判断
根据素数的概念，即除了 1 和此整数自身外，无法被其他自然数整除的数
cin >> n
flag = true
for i = 2 to n-1{
if(n % i == 0) flag = false
}
cout << flag
Eratosthenes筛法
素数的倍数一定不是素数
从 2 开始，把所有 2 的倍数都标记
然后进行 3 的遍历，执行相同的操作
到 4 的时候，由于之前已经将其标记，则对下一个素数 5 进行操作
依次进行上述操作，直到最后获得所有的素数
isPrime[] = true
primeCount = 0
For i = 2 .. N
If isPrime[i] Then
primeCount = primeCount + 1
multiple = 2
While (i * multiple ≤ N)
isPrime[i * multiple] = false
multiple = multiple + 1
End While
End If
End For
此算法有很多不足之处，比如 6，在素数为 2 时处理过一次，在为3的时候也处理了一次，重复了相同的操作。下面我们将介绍一个改进算法，欧拉筛法。
Eular质数筛法
规定每个合数只用最小的一个质因数去筛选，比如6有2和3两个因数，只用2进行筛选和置位操作，3的情况通过条件跳过。保证了每个合数只被他的最小素因子筛到一次，欧拉筛法的思想很巧妙，特别是应用在求一些积性函数的时候会比普通筛法更快，比如欧拉函数。
isPrime[] = true
primeList = []
primeCount = 0
For i = 2 .. N
If isPrime[i] Then
primeCount = primeCount + 1
primeList[ primeCount ] = i
End If
For j = 1 .. primeCount
If (i * primeList[j] > N) Then
Break
End If
isPrime[ i * primeList[j] ] = false
If (i % primeList[j] == 0) Then
Break
End If
End If
End For
Miller_Rabin质数检测
这种质数算法是基于费马小定理的一个扩展，首先我们要知道什么是费马小定理：
对于质数p和任意整数a，有a^p ≡ a(mod p)(同余)。
反之，若满足a^p ≡ a(mod p)，p也有很大概率为质数。
将两边同时约去一个a，则有a^(p-1) ≡ 1(mod p)
也即是说：假设我们要测试n是否为质数。我们可以随机选取一个数a，然后计算a^(n-1) mod n，如果结果不为1，我们可以100%断定n不是质数。
否则我们再随机选取一个新的数a进行测试。如此反复多次，如果每次结果都是1，我们就假定n是质数。
该测试被称为Fermat测试。
Fermat测试不一定是准确的，有可能出现把合数误判为质数的情况。
Miller和Rabin在Fermat测试上，建立了Miller-Rabin质数测试算法。
与Fermat测试相比，增加了一个二次探测定理：
如果p是奇素数，则 x^2 ≡ 1(mod p)的解为 x ≡ 1 或 x ≡ p - 1(mod p)
将这两条定理合起来，也就是最常见的Miller-Rabin测试。
但一次MR测试仍然有一定的错误率。为了使我们的结果尽可能的正确，我们需要进行多次MR测试，这样可以把错误率降低。
Miller-Rabin(n):
If (n <= 2) Then
If (n == 2) Then
Return True
End If
Return False
End If
If (n mod 2 == 0) Then
// n为非2的偶数，直接返回合数
Return False
End If
// 我们先找到的最小的a^u，再逐步扩大到a^(n-1)
u = n - 1; // u表示指数
while (u % 2 == 0)
u = u / 2
End While // 提取因子2
For i = 1 .. S // S为设定的测试次数
a = rand_Number(2, n - 1) // 随机获取一个2~n-1的数a
x = a^u % n
While (u < n)
// 依次次检查每一个相邻的 a^u, a^2u, a^4u, ... a^(2^k*u)是否满足二次探测定理
y = x^2 % n
If (y == 1 and x != 1 and x != n - 1) // 二次探测定理
// 若y = x^2 ≡ 1(mod n)
// 但是 x != 1 且 x != n-1
Return False
End If
x = y
u = u * 2
End While
If (x != 1) Then // Fermat测试
Return False
End If
End For
Return True
Wiki中的伪代码比上文中的简洁一些，并且有介绍了一些小技巧：比如如果n<2^64，只用选取a=2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37做测试即可
10^30的Miller-Rabin质数检测
如果用Java 的大数相关类型的话，就比较简单了，然而用C++，就。。。
#include <bits/stdc++.h>
#define MAXL 4
#define M10 1000000000
#define Z10 9
using namespace std;
const int zero[MAXL - 1] = {0};
struct bnum{
int data[MAXL]; // 断成每截9个长度
// 读取字符串并转存
void read(){
memset(data, 0, sizeof(data));
string buf;
cin >> buf;
int len = buf.length();
int i = 0, k;
while (len >= Z10){
for (k = len - Z10; k < len; ++k){
data[i] = data[i] * 10 + buf[k] - '0';
}
++i;
len -= Z10;
}
if (len > 0){
for (k = 0; k < len; ++k){
data[i] = data[i] * 10 + buf[k] - '0';
}
}
}
bool operator == (const bnum &x){
return memcmp(data, x.data, sizeof(data)) == 0;
}
bnum & operator = (const int x){
memset(data, 0, sizeof(data));
data[0] = x;
return *this;
}
bnum operator + (const bnum &x){
int i, carry = 0;
bnum ans;
for (i = 0; i < MAXL; ++i){
ans.data[i] = data[i] + x.data[i] + carry;
carry = ans.data[i] / M10;
ans.data[i] %= M10;
}
return ans;
}
bnum operator - (const bnum &x){
int i, carry = 0;
bnum ans;
for (i = 0; i < MAXL; ++i){
ans.data[i] = data[i] - x.data[i] - carry;
if (ans.data[i] < 0){
ans.data[i] += M10;
carry = 1;
}
else{
carry = 0;
}
}
return ans;
}
// assume *this < x * 2
bnum operator % (const bnum &x){
int i;
for (i = MAXL - 1; i >= 0; --i){
if (data[i] < x.data[i]){
return *this;
}
else if (data[i] > x.data[i]){
break;
}
}
return ((*this) - x);
}
bnum & div2(){
int i, carry = 0, tmp;
for (i = MAXL - 1; i >= 0; --i){
tmp = data[i] & 1;
data[i] = (data[i] + carry) >> 1;
carry = tmp * M10;
}
return *this;
}
bool is_odd(){
return (data[0] & 1) == 1;
}
bool is_zero(){
for (int i = 0; i < MAXL; ++i){
if (data[i]){
return false;
}
}
return true;
}
};
void mulmod(bnum &a0, bnum &b0, bnum &p, bnum &ans){
bnum tmp = a0, b = b0;
ans = 0;
while (!b.is_zero()){
if (b.is_odd()){
ans = (ans + tmp) % p;
}
tmp = (tmp + tmp) % p;
b.div2();
}
}
void powmod(bnum &a0, bnum &b0, bnum &p, bnum &ans){
bnum tmp = a0, b = b0;
ans = 1;
while (!b.is_zero()){
if (b.is_odd()){
mulmod(ans, tmp, p, ans);
}
mulmod(tmp, tmp, p, tmp);
b.div2();
}
}
bool MillerRabinTest(bnum &p, int iter){
int i, small = 0, j, d = 0;
for (i = 1; i < MAXL; ++i){
if (p.data[i]){
break;
}
}
if (i == MAXL){
// small integer test
if (p.data[0] < 2){
return false;
}
if (p.data[0] == 2){
return true;
}
small = 1;
}
if (!p.is_odd()){
return false; // even number
}
bnum a, s, m, one, pd1;
one = 1;
s = pd1 = p - one;
while (!s.is_odd()){
s.div2();
++d;
}
for (i = 0; i < iter; ++i){
a = rand();
if (small){
a.data[0] = a.data[0] % (p.data[0] - 1) + 1;
}
else{
a.data[1] = a.data[0] / M10;
a.data[0] %= M10;
}
if (a == one){
continue;
}
powmod(a, s, p, m);
for (j = 0; j < d && !(m == one) && !(m == pd1); ++j){
mulmod(m, m, p, m);
}
if (!(m == pd1) && j > 0){
return false;
}
}
return true;
}
相关参考：
Miller-Rabin算法+大数算法参考了f_zyj的博文Miller-Rabin算法+大数算法
转载请注明出处：http://www.cnblogs.com/ygdblogs
