miller_rabbin 算法

miller_rabbin算法

这是一个质数判断算法，它采用随机算法能够在很短的时间里判断一个数是否是质数。
首先，根据费马定理，若$p$为质数，取a<p,则$a^{p-1}\%p=1$.但反过来则不成立。
因为有一类数——卡迈克尔数，它不是质数，但也满足$a^{p-1}\%p=1$。
如何将卡迈克尔数甄别出来。我们要用到二次探测方法。
设$p$为我们要判断的数，将$p-1$分解为$2^r\ast k$，则有$a^{p-1}=(a^k{)}^{2^r}$
则我们可以先求出$a^k$,然后将其不断平方，并模p。如果第一次出现模$p$的值为1，那么检测上一次的值是否为$-1$,如果是,则没有毛病；否则可以断定不是质数。
证明:
令$$y^2\%p=1$$
则等价于：$$y^2-1=k\ast p$$
等价于：$$(y-1)\ast (y+1)=k\ast p$$
如果$p$是质数，则$$\{\begin{array}{l}y-1=0\\ y+1=p\end{array}$$
即y的值只能为$1$或$p-1$。
所以，如果$y\ne 1$且$y\ne p-1$,则可以断定$p$不是质数。
而如果做r次平方，模$p$的值都不为1，则$p$显然不是质数了。
那如果在平方r次的过程中，最后的结果为1了。而且第一次结果为1时，之前的结果为-1。是不是就能保证$p$是质数呢?也不能百分百地肯定。但是相比于仅通过费马定理来判断，它更严格一些了。这种方法称为二次检测。
于是，我们可以多次选择$a$，如果选了许多$a$来做上述的操作，都不能判断出$p$是合数，则$p$多半就是一个质数了。
miller_rabbin算法是一个随机算法，并不能完全保证结果准确，但是出错的概率非常小。我们取a的次数越多，出错的概率越小。一般情况下，取20~30次a，正确率已经可以达到99.9999999%了。

FZU1649
给出若干整数，判断它是否是素数。
参考代码如下：

#include<iostream>
#include<cstring>
#include<ctime>
#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;
#define LL unsigned long long int
LL n,k;
int t,r;
bool notprime;
LL ksc(LL a,LL b,LL p)
{
    LL res=0;
    while(b)
    {
        if(b&1)res=res+a;
        a=a+a;
        if(res>=p)res-=p;
        if(a>=p)a-=p;
        b>>=1;
    }
    return res;
}
LL ksm(LL a,LL b,LL p)
{
    LL res=1;
    while(b)
    {
        if(b&1)res=ksc(res,a,p);
        a=ksc(a,a,p);
        b>>=1;
    }
    return res;
}
int main()
{
    srand(time(NULL));
    while(scanf("%I64d",&n)!=-1)
    {
     if(n==2||n==3)
    {
         printf("It is a prime number.\n");
         continue;
    }
    else
    {
        if(n==1||n&1==0)
        {
            printf("It is not a prime number.\n");
            continue;
        }   
    }
    
    r=0;
    k=n-1;
    while(!(k&1))
    {
        r++;
        k>>=1;
    }
    int t=20;
    LL res=1,lastres=0;
    notprime=0;
    for(int id=0;id<=10;id++)
    {
        res=rand()%(n-1)+1;
        res=ksm(res,k,n);
        if(res==1)
        continue;
        lastres=res;
        for(int i=0;i<r;i++)
        {
            res=ksc(res,res,n);
            if(res==1)
            {
                if(lastres!=n-1) notprime=1;
               break;
            }
            lastres=res;
        }
        if(notprime==1)break;
        if(res!=1)
        {
            notprime=1;
            break;
        }
    }
    if(notprime==1)
        printf("It is not a prime number.\n");
    else 
        printf("It is a prime number.\n");
    }
    return 0;
}