素性测试的Miller-Rabin算法完全解析 （C语言实现、Python实现）

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

以下为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 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()