一、Miller-Rabin算法介绍

Miller-Rabin 算法（强伪素数检测）利用原理：若n为素数，则根据费马小定理：若n是一个素数，且(a,n)=1，则$a ^ {n-1}=1modn$。且1$mod$n只有1和-1 两个平方根。

首先，假设n为一个非2素数，则n-1为一个偶数，任何一个偶数可以写成$2 ^ k m$的形式，其中m为一个奇数。随机选取a，满足$1\le a\le n-1$。将$a^m(mod\ n)$赋值给一个新的变量b。如果，$b\equiv1(mod\ n)$，则表明假设成立，n 为大于 2 的素数。否则，进行k次循环，每次循环判断，若$b\equiv -1(mod\ n)$，则n为大于2的素数，否则将$b ^ 2(mod\ n)$赋值给b，继续循环。若最终没有得到$b\equiv -1(mod\ n)$，则n为一个合数。

且当n为2时，n-1=1，k=0，m=1，a=1，则b=1满足$b\equiv 1(mod\ n)$说明Miller-Rabin算法对n为2同样适用。

$\textcolor{red}{但该算法有概率出错}$：

1、在取a时，若a取到1，则无论n是否为素数，$a^m(mod\ m)$一定等于1，即$b\equiv1(mod\ n)$必定成立，导致将其判断为素数。因a若无法取到1，只会影响到2是否为素数的判定，所以，在算法实现时可令a无法取到1，且在生成a前判断是否为2。以使算法成功判断素数。
2、当a取得 n-1 时，则因 m 必定为奇数，所以，$a ^ m+1=(n-1) ^ m+1=n((n-1) ^ {m-1}-(n-1) ^ {m-2}\dots\dots+1)$此时$b\equiv -1(mod\ n)$成立，会将合数误判为素数。
3、因费马小定理的逆命题不成立，存在合数。如：$341=11\times31$，若取a=2，$2 ^ {340}\equiv1(mod\ 341)$。

Miller-Rabin算法出现错误的概率约为$\frac{1}{4}$，根据上述分析认为错误原因主要在于a的随机选取上，综上，可以通过多次随机选择a避免1，2情况的发生，以及降低发生3的概率，但多次取a会使判断时间变长，因$(\frac{1}{4}) ^ 6 \approx 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")

不知所谓の博客

Miller-Rabin算法素检测介绍及简单实现