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

 

以下为C语言源代码：

#include <stdio.h>
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; i<s; i++) { //进行s次二次探测
        newX=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; i<N+10; i++) {
        n=i;
        //check(n);        
        //result=isPrime(n);
        result=isPrimeByMR(n);
        if(result==TRUE)
            printf("%llu\n",n);

    }
    return 0;
}

 

以下为Python语言源代码：

import math
import randomdef isPrime(n): #这是传统的方法，用于与Miller-Rabin算法的结果进行对比。if (n & 1) == 0 or n % 5 == 0:return Falsebnd = int(math.sqrt(n))for i in range(2, bnd + 1):if n % i == 0:return Falsereturn True# 实现(a*b)%c,用类似快速幂的方式实现。
# 将模运算下的乘法用加法实现的思想是时间换空间。
# 这可以作为典型的时间换空间的例子。
def quickMult(a, b, c):result = 0while b > 0:if b & 1:result = (result + a) % ca = (a + a) % cb >>= 1return result#因为Python支持大整数运算，所以在此也可以不用快速乘法，而直接使用乘法*。
def quickPower(a, b, c):result = 1while b > 0:if (b & 1):# result=result*a%c #此为直接乘法result = quickMult(result, a, c)# a=a*a%c #此为直接乘法a = quickMult(a, a, c)b >>= 1return 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 - 1while (d & 1) == 0: #将d中因子2全部提取出来。s += 1d >>= 1x = 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 = newXif x != 1:  # 用费马小定理的逆否命题判断，此时x=a^(n-1) (mod n)，那么n确定为合数。return Falsereturn True  # 用费马小定理的逆命题判断。能经受住考验至此的数，大概率为素数。# 经过连续特定次数的Miller-Rabin测试后，
# 如果返回值为TRUE表示n为素数，返回值为FALSE表示n为合数。
def isPrimeByMR(n):if ((n & 1) == 0 or n % 5 == 0):return Falsefor i in range(100):if MillerRabinPrimeTest(n) == False:return Falsereturn 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 = 1000000000000000000for 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()