Miller-Rabin素性测试

原文地址：https://www.cnblogs.com/Norlan/p/5350243.html

素数：若一个数x的约数仅仅只有1和他本身，则称之为素数，注意，1不是素数。

素数的性质：

1.素数无穷多，不存在最大素数

证明：假设最大素数存在且为p，我们可以用旧素数去构造新素数，设x=2*3*5*7*11*......*p+1，x不会被用来构造他的任何一个素数整除，当然由于素数之间也是互质的，没有被用来构造他的素数也不能整除他，只不过用来构造他的素数整除他的余数为1罢了。

2.存在任意长的一段连续区间，区间内的数全部都是合数

证明：当0<a<=n时，n!+a%a=0。而n!+2,n!+3,n!+4,......,n!+n都是合数（可分别被2,3,4,......,n整除)，我们构造出了一段n-1的连续合数区间。

注：相邻素数之间的间隔可以是任意的可能是有问题的说法，通过上述的构造方法，我们可以构造出一个连续合数区间，使得它的长度>=n-1，但是，具体这个长度是多少是未知的。以N=5为例，N!+2,N!+3,N!+4,N!+5，都是合数，但是N!+1以及N!+6都是合数，不再往两边继续举例，因此这个办法无法证明相邻素数间隔任意。

3.所有大于2的素数都可以唯一表示成两个平方数之差

证明：首先我们知道，大于2的素数都是奇数，那么我们假设这个数为2*n+1。(n+1)^2=n^2+2*n+1，而2*n+1=(n+1)^2-n^2。

唯一性：若素数p=a^2-b^2，那么p=(a+b)*(a-b)，由于p是素数，那么只能是a+b=p且a-b=1。

4.当n为大于2的整数时，2^n+1和2^n-1这两个数，至少有一个是合数

证明：2^n%3!=0。分情况讨论：

若2^n%3=1，那么2^n-1%3=0，为合数。

若2^n%3=2，那么2^n+1%3=0，为合数。

5.费马小定理：如果p是素数，a是小于p的正整数，那么a^(p-1)%p=1

证明：1.如果p是素数，那么对于任意一个小于p的正整数a，有a、2a、3a、4a、......、(p-1)a，这些数分别取余p后的结果是一个1~p-1的全排列。举个例子：5是素数，3, 6, 9, 12除以5的余数分别为3, 1, 4, 2。

反证法：如果结论不成立，我们可以知道有两个小于p的数m和n使得m*a%p=n*a%p。假设m>n，我们可以构造一个

(m-n)*a%p=0，由于p是素数，只可能(m-n)%p=0或者a%p=0，但显然的是由于m-n和a都是小于p的，因此不成立。

那么由上述结论我们可以得到a*2a*......*(p-1)a%p=(p-1)!，两边同时除去(p-1)!，即可得到a^(p-1)%p=1。

下面进入素性测试部分：

由费马小定理我们可以知道：如果p是素数，a是小于p的正整数，那么a^(p-1)%p=1。

费马素性测试：如果我们反过来考虑，如果a是小于p的正整数，a^(p-1)%p=1，是否可以得到p是素数呢。

以a=2时为例，2^4%5=1，2^340%341=1，5为素数但是341=11*31。显而易见的是，如果a^(p-1)%p！=1，那么可以直接肯定p不可能为素数，而能通过这种测试的我们称之为以2为底的伪素数。

那么自然会联想到，如果用所有小于p的数去测试p的素性呢。然而数据表明，有些合数强大到可以通过所有的这种测试。前10亿个自然数中有600个之多。

Miller-Rabin素性测试：

原理：如果p是素数，x是一个小于p的正整数，若x^2%p=1，我们可以推出(x-1)(x+1)%p=0，那么x的取值只能是1或者p-1。

具体测试方法：以341为例：

A.由于，我们假设341为素数，那么我们可以得到结论：或者。

B.进行求证，如果满足（正确），继续按照分解的原理求证或者。

如果，我们是不可以继续往下分解的，那么我们就可以判定341是以2为底的强伪素数。

C.我们发现，所以可以得到结论341并不是素数。

这种测试方法的出错概率比费马素性测试要低得多，数据表明：第一个以2为底的强伪素数为2047。第一个以2和3为底的强伪素数则大到1373653。

具体代码（核心部分）：

#include <iostream>
using namespace std ;
typedef long long ll;

ll pow_mod(ll a,ll x,ll n){//快速幂取模
    ll r=1;
    while(x){
        if(x&1){
            r=r*a%n;
        }
        a=a*a%n;
        x>>=1;
    }
    return r;
}

bool test(ll a,ll x,ll n){//a^x%n  x=n-1
    if(!(n&1))return false;//大于2的偶数不可能是素数
    while(!(x&1))x>>=1;//将x的所有2全部提取
    ll res=pow_mod(a,x,n);
    while(x!=n-1&&res!=n-1&&res!=1){
        res=res*res%n;
        x<<=1;
    }//3*3%5==4  
    return res==n-1||(x&1)==1;
    //从小求到大可以减少运算
    //如果说res==1，并且n要为素数，那么必须在x还为奇数也就是第一次进行素性判断的时候，res==1才行
    //由梅氏素性测试的原理，若t*t%n==1,则t==1或t==n-1
    //梅氏素性测试要求如果t*t%n=1的时候  才能去判定是否由t==1或者t==n-1
    //如果遇到了t*t%n=n-1的情况，此时是无法向下继续判断的，因此按照素性测试的要求，我们直接返回是强伪素数即可
    //t%n=n-1时，t*t%n=1，往后都是1
    //代码是从小判到大的，强伪素数要求要么一开始t就是1要么t在变换中途变成了n-1
}

bool isPrime(ll n){
    int a[]={2,3,5,7};
    for(int i=0;i<=3;i++){
        if(n==a[i])return true;
        if(!test(a[i],n-1,n))return false;//未通过以2 3 5 7四个数中任意一个为底时的素性测试
    }
    return true;
}

int main(){
    ll ans;
    cin>>ans;
    cout<<(isPrime(ans)?"Yes":"No");
}

防溢出的处理（利用快速乘法，加直接相乘变为利用二进制的加法）：

举个例子：3*6

A.将6转换为二进制110

B.那么我们可以按照二进制对结果进行分解：

C.加法的过程中取模就可以防止溢出

先对快速幂的部分做如下调整：

ll mul_mod(ll a,ll b,ll n){
    ll r=0;
    while(b){
        if(b&1){
            r=(r+a)%n;    
        }
        a=(a+a)%n;
        b>>=1;
    }
    return r;    
}

ll pow_mod(ll a,ll x,ll n){//快速幂取模
    ll r=1;
    while(x){
        if(x&1){
            r=mul_mod(r,a,n);
        }
        a=mul_mod(a,a,n);
        x>>=1;
    }
    return r;
}