miller_robin大素数判定

参考了ACdreamer大神的博客http://blog.csdn.net/acdreamers/article/details/7913786

在51nod上看了个10^30范围的素数判定，打表肯定是不行了，应该用miller-robin素数判定加java的BigInteger

首先基于fermat小定理就是gcd(a,p)=1时，a^(p-1)%p=1,所以在生成rand(2,p-1)的随机数后,如果输入的是个素数,那么必有pow(a,p-1)%p=1,这里采用快速幂取模，在一次探查后，若结果不为1，则p必然不是一个素数，否则有75%的概率是素数，那么探查10次有

1-(0.75)^10=0.94369的概率是素数（好像是这样）。当然可以多探查到20次，成功率0.996828

在这里还有一个优化，就是根据上面博客叙述的二次探查定理x^2==1(mod p)在p为素数的情况下x必为1或p-1,那么就有一个优化，先滤出p-1的2次幂，假设p-1=m*2^sqr，先计算x=pow(a,m)%p，之后每次平方看是否存在y=x*x%p==1但是x又不为1，p-1的，这就说明p必不是素数

代码如下：

import java.io.*;
import java.util.*;
import java.math.BigInteger;
public class minller{
    public static final int iter=10;
    public static BigInteger fast_fac(BigInteger a,BigInteger n,BigInteger mod){
        BigInteger TWO=BigInteger.valueOf(2),ans;
        a=a.mod(mod);
        ans=BigInteger.ONE;
        while (!(n.equals(BigInteger.ZERO)))
        {
            if ((n.mod(TWO)).equals(BigInteger.ONE)){
                ans=(ans.multiply(a)).mod(mod);
                n=n.subtract(BigInteger.ONE);
            }
            n=n.divide(TWO);
            a=(a.multiply(a)).mod(mod);
        }
        return ans;
    }
    public static boolean millerRobin(BigInteger p){
        BigInteger TWO=BigInteger.valueOf(2),a,x,n,y=BigInteger.ZERO,p_1;
        Random rnd=new Random();

        p_1=p.subtract(BigInteger.ONE);
        if (p.equals(BigInteger.ONE)) return false;
        if (p.equals(TWO)) return true;
        if ((p.mod(TWO)).equals(BigInteger.ZERO)) return false;
        
        int MAX,sqr=0;
        //滤出2的次方
        n=p.subtract(BigInteger.ONE);
        while ((n.mod(TWO)).equals(BigInteger.ZERO)){
            ++sqr;
            n=n.divide(TWO);
        }
        
        if (n.compareTo(BigInteger.valueOf(10000))==1)
            MAX=10000;
        else MAX=n.intValue();
        
        
        for (int i=0;i<iter;++i){
            a=BigInteger.valueOf(rnd.nextInt(MAX-2)+2);
            x=fast_fac(a,n,p);
            if (x.equals(BigInteger.ONE)) continue;
            for (int j=0;j<sqr;++j){
                y=(x.multiply(x)).mod(p);
                if (y.equals(BigInteger.ONE) && !x.equals(BigInteger.ONE) && !x.equals(p_1))
                    return false;
                x=y;
            }
            if (y.equals(BigInteger.ONE)==false)
                return false;
        }
        return true;
    }
    public static void main(String[] argv){
        Scanner cin=new Scanner(System.in);
        while (cin.hasNextBigInteger()){
            BigInteger n=cin.nextBigInteger();
            if (millerRobin(n))
                System.out.println("Yes");
            else System.out.println("No");
        }
    }
}