米勒拉宾算法（素性测试）

米勒拉宾素性测试

对于一个数n，如果想要判断它是否为素数，常规的方法为试除法。即，让n依次除以2到sqrt（n）以内的整数。如果有出现除尽的情况，则为合数。
该方法的时间复杂度为O（sqrt（n））在面对n为长整型的时候有可能超出时间要求。

因此普遍采用米勒拉宾算法进行素性判定。在此之前介绍一种伪素数判定方法——小费马定理。
小费马定理为：若有素数p，则对任意的数a( a为正整数 ，且a < p)，$a^{p-1}\equiv 1(mod$ $p)m$ 。反之，若有任意的a( a为正整数 ，且a < p)使得p不满足 $a^{p-1}\equiv 1(mod$ $p)$ ，p一定为合数。
可以发现若是能够举出所有的a，都能满足上式，是不是就说明p是素数呢？其实不是因为有一类合数也可以做到这一点，这一类合数叫做 Carmichael 数。前三个这样的数是561 ，1105，1729。
这样的数真让人不爽。所以采用这种方法测出来的所谓素数是不一定的。叫做伪素数。哪怕你枚举出所有的a，也不可避免。

在小费马定理的基础上有人设计出米勒拉宾随机素数测试法。可以大大的提高检测素数的正确性，但是同样并非一定正确，错误可能性却小到可以接受。

该方法同样利用了枚举多个a的做法，以提高算法的可靠性，对于每一个a，又采用了特殊的方法处理。这基于另外一个定理：如果p是素数，x是小于p的正整数，且$x^{2}mod$ $p=1$，那么要么$x=1$，要么$x=p-1$。该定理证明如下：如果p为素数，x是小于p的正整数， 且$x^{2}mod$ $p=1$ ，说明p能够整除$（x+1）（x-1）$。但是p是素数，那么只可能是$x-1$能被$p$整除(此时$x=1$)或$x+1$q能被$p$整除(此时 $x=p-1$)。

判断一个数是不是素数光靠上面的方法是不可靠的，因为p如果是合数的话，也有可能有$x^2\equiv 1mod(p)$ 且 $x=1$或者 $x=p-1$；但是多排除几次p不为合数的话，就增大了p是素数的可能性 ，这是这个算法的核心思想。
因此例如341这个数。可知 $(2^{340})\equiv 1(mod$ $340)$; $(2^{170})\equiv 1(mod$ $340)$； 但是发现 $2^{85}mod$ $341=32$。这足以证明341是一个合数，而不是一个素数。

首先判断要判断的数n是不是2，在判断n是不是奇数。然后尽可能的在令$d=n-1$，在$d$中除去2，使得$n=d\ast (2^t)$,d为奇数，t的值并不关心。如果n是一个素数，那么或者$a^{d}mod$ $n=1$，或者存在某个i使得$a^{d\ast 2^i}mod$ $n=n-1(0<=i<r)$（注意i可以等于0，这就把$a^{d}mod$ $n=n-1$的情况统一到后面去了）。

求$a^{d}mod(n)$的算法以及求$d^2$的算法是采用的快速幂取模算法。但是在d为long long的情况下有可能乘法溢出。有更加优秀的算法存在。

#include <iostream>
#include <cstdio>
using namespace std;
typedef long long LL;

LL mulmod( LL a, LL b , LL p )
{
    LL  d =1;
    a = a%p;
    while( b>0 )
    {
        if(b&1)
            d = (d*a)%p;
        a = (a*a)%p;
        b>>=1;
    }
    return d;
}

bool witness( LL a,LL n)
{
    LL d = n-1 ;
    if( n ==2 ) return true ;
    if( !(n&1) ) return false ;
    while(!(d&1)) d = d/2;
    LL t = mulmod(a,d,n);
    while((d!=n-1) && (t!=1)&&(t!=n-1))
    {
        t = mulmod( t ,2,n);
        d=d<<1;
    }
    return (t==n-1)||(d&1);
}

bool isprime( LL n)
{
    int a[3] = {2,7,61};
    for(int i=0;i<3;i++)
        if(!witness(a[i],n))
            return false;
    return true;
}
int main()
{
    LL s;
    cin>>s;
    if(isprime(s))
        cout<<"YES";
    else
        cout<<"NO";
    return 0;
}

转自：https://blog.csdn.net/qq_37957829/article/details/77335072