python连续质数计算代码分析,素性测试的Miller-Rabin算法完全解析 （C语言实现、Python实现）...

因为文中存在公式，只能用图片方式上传了！

以下为C语言源代码：

#include typedef long long unsigned LLU;

typedef int BOOL;

#define TRUE 1

#define FALSE 0

BOOL isPrime(LLU n) {  //这是传统的方法，用于与Miller-Rabin算法的结果进行对比。

if( (n&1)==0 || n%5==0)

return FALSE;

LLU i,bnd;

bnd=sqrt(n);

for(i=2; i<=bnd; i++)

if(n%i==0)

return FALSE;

return TRUE;

}

//实现(a*b)%c,用类似快速幂的方式实现。

//将模运算下的乘法用加法实现的思想是时间换空间。

//这可以作为典型的时间换空间的例子。

LLU quickMult(LLU a,LLU b,LLU c){

LLU result=0;

while(b>0) {

if(b&1)

result=(result+a)%c;

a=(a+a)%c;

b>>=1;

}

return result;

}

//请注意，计算快速幂时，因为变量的数据类型为long long,是64bits长度的整数。

//a的值不能大于(2的32次方-1)≈4000000000

//否则，计算a*a时会越界，超过64bits整数的表示范围。

//例如a,b,c分别为：7087881096594 10000000000036 10000000000037时，

//正确结果 (a^b)%c=1，但是，在此计算的结果不为1。

//因此，可以改进为用“快速乘”代替*乘法运算符，能更加充分第利用64的long long unsigned型的存储空间。

//这个例子可以作为算法改进从O(n)到O(lg(n))的进步的例子。

LLU quickPower(LLU a,LLU b,LLU c) {

LLU result=1;

while(b>0) {

if(b&1)

//result=result*a%c; //此为直接乘法

result=quickMult(result,a,c);

//a=a*a%c; //此为直接乘法

a=quickMult(a,a,c);

b>>=1;

}

return result;

}

//如果返回值为TRUE表示n为素数，返回值为FALSE表示n为合数。

BOOL MillerRabinPrimeTest(LLU n) {

LLU d,x,newX,a=1;

int i;

for (i=0; i<4; i++)

a*=rand();

a=a%(n-3)+2;//随机第选取一个a∈[2,n-2]

//printf("随机选取的a=%lld\n",a);

int s=0;//s为d中的因子2的幂次数。

d=n-1;  //d取初值n-1

while( (d&1)==0) {//将d中因子2全部提取出来。

s++;

d>>=1;

}

x=quickPower(a,d,n);

for(i=0; inewX=quickPower(x,2,n);

if(newX==1 && x!=1 && x!=n-1)

return FALSE; //用二次定理的逆否命题，此时n确定为合数。

x=newX;

}

if(x!=1)

return FALSE;   //用费马小定理的逆否命题判断，此时x=a^(n-1) (mod n)，那么n确定为合数。

return TRUE; //用费马小定理的逆命题判断。能经受住考验至此的数，大概率为素数。

}

//经过连续特定次数的Miller-Rabin测试后，

//如果返回值为TRUE表示n为素数，返回值为FALSE表示n为合数。

BOOL isPrimeByMR(LLU n) {

if((n&1)==0 || n%5==0)

return FALSE;

int i;

for (i=0; i<100; i++)

if(MillerRabinPrimeTest(n)==FALSE)

return FALSE;

return TRUE;

}

//对比传统方法和Miller-Rabin算法的结果

void check(LLU n) {

char isRight;

BOOL resultA,resultB;

resultA=isPrime(n);

resultB=isPrimeByMR(n);

if(resultA==resultB) {

isRight='V';

printf("%c,%llu %d %d\n",isRight,n,resultA,resultB);

} else {

isRight='X';

printf("%c,%llu %d %d\n",isRight,n,resultA,resultB);

}

}

//测试任务：在本人笔记本电脑上测试，N=10的18次方至N+10之间的质数为：

//1000000000000000003

//1000000000000000009

//Miller-Rabin算法的速度：0.22秒

//常规算法：22秒

int main() {

srand(time(NULL));

LLU i,n,N;

N=1000000000000000000;

BOOL result;

for(i=N; in=i;

//check(n);

//result=isPrime(n);

result=isPrimeByMR(n);

if(result==TRUE)

printf("%llu\n",n);

}

return 0;

}

以下为Python语言源代码：

import math

import random

def isPrime(n): #这是传统的方法，用于与Miller-Rabin算法的结果进行对比。

if (n & 1) == 0 or n % 5 == 0:

return False

bnd = int(math.sqrt(n))

for i in range(2, bnd + 1):

if n % i == 0:

return False

return True

# 实现(a*b)%c,用类似快速幂的方式实现。

# 将模运算下的乘法用加法实现的思想是时间换空间。

# 这可以作为典型的时间换空间的例子。

def quickMult(a, b, c):

result = 0

while b > 0:

if b & 1:

result = (result + a) % c

a = (a + a) % c

b >>= 1

return result

#因为Python支持大整数运算，所以在此也可以不用快速乘法，而直接使用乘法*。

def quickPower(a, b, c):

result = 1

while b > 0:

if (b & 1):

# result=result*a%c #此为直接乘法

result = quickMult(result, a, c)

# a=a*a%c #此为直接乘法

a = quickMult(a, a, c)

b >>= 1

return result

#如果返回值为TRUE表示n为素数，返回值为FALSE表示n为合数。

def MillerRabinPrimeTest(n):

a = random.randint(2,n-2) #随机第选取一个a∈[2,n-2]

# print("随机选取的a=%lld\n"%a)

s = 0 #s为d中的因子2的幂次数。

d = n - 1

while (d & 1) == 0: #将d中因子2全部提取出来。

s += 1

d >>= 1

x = quickPower(a, d, n)

for i in range(s): #进行s次二次探测

newX = quickPower(x, 2, n)

if newX == 1 and x != 1 and x != n - 1:

return False #用二次定理的逆否命题，此时n确定为合数。

x = newX

if x != 1: # 用费马小定理的逆否命题判断，此时x=a^(n-1) (mod n)，那么n确定为合数。

return False

return True # 用费马小定理的逆命题判断。能经受住考验至此的数，大概率为素数。

# 经过连续特定次数的Miller-Rabin测试后，

# 如果返回值为TRUE表示n为素数，返回值为FALSE表示n为合数。

def isPrimeByMR(n):

if ((n & 1) == 0 or n % 5 == 0):

return False

for i in range(100):

if MillerRabinPrimeTest(n) == False:

return False

return True

#对比传统方法和Miller-Rabin算法的结果

def check(n):

resultA = isPrime(n)

resultB = isPrimeByMR(n)

if resultA == resultB:

isRight = 'V'

print("%c,%u %d %d" % (isRight, n, resultA, resultB))

else:

isRight = 'X'

print("%c,%u %d %d" % (isRight, n, resultA, resultB))

# 测试任务：在本人笔记本电脑上测试，N=10的18次方至N+10之间的质数为：

# 1000000000000000003

# 1000000000000000009

# Miller-Rabin算法的速度：0.22秒

# 常规算法：22秒

# python的常规算法，耗时更长。

def main():

# freopen("result.txt","w",stdout)

random.seed()

# res=quickPower(2,10,7)

# print("%u"%res)

N = 1000000000000000000

for i in range(N, N + 10):

n = i

#check(n)

# n=int(input())

#result=isPrime(n)

result = isPrimeByMR(n)

if result == True :

print("%u" % n)

return 0

# print(isPrime(1000000000000000003))

# a, b, c = [int(e) for e in input().split(" ")]

# print(quickPower(a,b,c))

# # 7087881096594 10000000000036 10000000000037 =1

main()