目录
一、Miller-Rabin算法介绍
Miller-Rabin 算法(强伪素数检测)利用原理:若 n 为素数,则根据费马小定理:若n是一个素数,且(a,n)=1,则
且1𝑚𝑜𝑑𝑛只有 1 和-1 两个平方根。
首先,假设n为一个非2素数,则 n-1 为一个偶数,任何一个偶数可以写成的形式,其中m为一个奇数。随机选取a,满足1 ≤ 𝑎 ≤ 𝑛 − 1。将
赋值给一个新的变量b。如果,𝑏 ≡ 1(𝑚𝑜𝑑𝑛),则表明假设成立,n 为大于 2 的素数。否则,进行k次循环,每次循环判断,若𝑏≡−1(𝑚𝑜𝑑𝑛),则n为大于2的素数,否则将
赋值给b,继续循环。若最终没有得到𝑏≡−1(𝑚𝑜𝑑𝑛),则n为一个合数。
且当n为2时,n-1=1,k=0,m=1,a=1,则b=1满足𝑏≡1(𝑚𝑜𝑑𝑛),说明Miller-Rabin算法对n为2同样适用。
但该算法有概率出错:
1、在取 a 时,若 a 取到 1,则无论 n 是否为素数,一定等于 1,即𝑏 ≡ 1(𝑚𝑜𝑑𝑛)
必定成立,导致将其判断为素数。因 a 若无法取到 1,只会影响到 2 是否为素数的判定,所以,在算法实现时可令 a 无法取到 1,且在生成 a 前判断是否为2。以使算法成功判断素数。
2、当 a 取得 n-1 时,则因 m 必定为奇数,所以,
此时𝑏 ≡ −1(𝑚𝑜𝑑𝑛)成立,会将合数误判为素数。
3、因费马小定理的逆命题不成立,存在的合数。如:341=,若取a=2,
。
Miller-Rabin算法出现错误的概率约为,根据上述分析认为错误原因主要在于a的随机选取上,综上,可以通过多次随机选择a避免1,2情况的发生,以及降低发生3的概率,但多次取a会使判断时间变长,因
≈ 0.0002已经较小,所以在实现中可以取6次a分别判断。
二、快速乘与快速幂
(一)、分析
在C++中,为测试较大的数据,采用long long类型,long long类型为8字节,即可计算64bit的数据,但在幂计算时会出现溢出现象,所以介绍一种快速乘与快速幂的运算方式,利用分治的思想进行幂运算,利用加法代替乘法,以防止溢出的发生。具体分析可参考大佬的文章
(二)、快速乘代码
long long mod_mul(long long a, long long b,long long n)
{
long long result = 0;
while (b > 0)
{
if (b & 1)
result = (result + a) % n;
a = (a + a) % n;
b = b >> 1;
}
return result;
}
(三)、快速幂代码
long long mod_exp(long long a, long long b, long long n)
{
long long result = 1;
while (b > 0)
{
if (b & 1)
result = mod_mul(result, a, n);
a = mod_mul(a, a, n);
b = b >> 1;
}
return result;
}
三、完整Miller-Rabin代码
(一)、C++代码
#include <iostream>
using namespace std;
long long mod_mul(long long a, long long b, long long n)//快速乘法
{
long long result = 0;
while (b > 0)
{
if (b & 1)//判断是否为偶数
result = (result + a) % n;
a = (a + a) % n;
b = b >> 1;//除二操作
}
return result;
}
long long mod_exp(long long a, long long b, long long n)//快速幂
{
long long result = 1;
while (b > 0)
{
if ((b & 1) > 0)
result = mod_mul(result, a, n);
a = mod_mul(a, a, n);
b = b >> 1;//除二操作
}
return result;
}
bool isprime(long long n)
{
int k = 0;
long long p = n - 1;
while ((p & 1) == 0)//判断是否为奇数
{
p = p >> 1;//除二操作
k++;
}
for (int i = 0; i < 6; i++)
{
long long a = rand() % (n - 1 - 1 + 1) + 1;
long long b = mod_exp(a, p, n);
bool flag = false;
if (b == 1)
continue;
for (int j = 0; j < k; j++)
if ((b + 1) % n == 0)
{
flag = true;
break;
}
else
b = (b * b) % n;
if (flag)
continue;
return false;
}
return true;
}
int main()
{
long long N;
cin >> N;
if (isprime(N))
cout << N << " is prime!";
else
cout << N << " is composite!";
}
(二)、python代码
import random
def isprime(n):
k,p=0,n-1
while (p&1)==0:
p=p>>1
k+=1
for j in range(6):
a=random.randint(1,n-1)
b=pow(a,p,n)
flag=0
if b==1:
continue
for i in range(k):
if (b+1)%n==0:
flag=1
break
else:
b=(b*b)%n
if flag==1:
continue
else:
return False
return True
n=int(input())
if isprime(n):
print(n,"is prime”)
else:
print(n,"is composite")