大素数测试的Miller-Rabin算法

Miller-Rabin算法本质上是一种概率算法，存在误判的可能性，但是出错的概率非常小。出错的概率到底是多少，存在严格的理论推导。

费尔马小定理

• 如果$p$是质数且$(a,p)=1$，则有$a^{p-1}\equiv1({\rm mod}\,{p})$。
  当然反过来不一定成立。即当$a^{p-1}{\%}{p}=1$时，$p$未必是质数。但是这个概率比较小。所以利用费尔马小定理来检测素数，不能保证时刻都对，只能保证出错的概率比较小。
  给定正整数n，问n是否为质数（显然只需判断正奇数），最基本的做法就是计算$2^{n-1}\%n$是否为1。如果不是1，n肯定为合数；否则，n可能为质数。

有限域上的平方根定理

• 如果$p$是一个奇质数且$e\geq1$，则方程
  $$x^2 \equiv 1({\rm mod}\,{p^e})$$
  仅有两个根$x=1$或者$x=-1$，注意到在模$p$的意义下，$x=-1$等价于$x=p-1$，$\pm1$也称为1的平凡平方根
• 很容易有一个推论，如果对模n存在1的非平凡平方根，n一定是合数

Miller-Rabin算法

利用上面两个定理，就可以构造出Miller-Rabin算法。考虑到n肯定是奇数(偶数的情况自己想去)，则n一定可以表示为$n-1=2^s*d$，其中$s\geq1$且$d$是奇数。则

$$a^{n-1}=a^{2^{s}*d}=(((a^d)^2)^{...})^2$$
也就是说，$a^{n-1}$相当于$a^d$平方若干次。例如当$n=7$时，$a^{n-1}$就是$a^6$，就是$a^3$的平方。当$n=13$时，$a^{n-1}$就是$a^{12}$，就是$a^3$的平方的平方。
以$n=13$的情况进行说明（所有运算都是在模n的意义下，以下的文字说明省略了这一点），任取一个$a,1\lt{a}\lt13$，计算$a^3$，再将其平方一次得到$a^6$，注意到$a^3$是$a^6$的平方根（废话），根据平方根定理的推论，如果$a^6=1且a^3\neq\pm1$，则$n$肯定是合数。将$a^6$平方一次得到$a^{12}$，同样，如果$a^{12}=1且a^6\neq\pm1$，则$n$肯定是合数。最后，根据费尔马小定理，如果$a^{12}\neq1$，则$n$肯定是合数。否则，$n$有极大概率为质数。
为了增加得到正确判断的概率，可以将$a$重复取不同的值，对每一个$a$验证一次$a^d$到$a^{n-1}$的过程。不过，考虑到ACM的特殊性，测试数据应该不会选择伪素数特别是强伪素数。所以很多题目的AC程序实际上只来一次即可。

typedef long long llt;
//利用二进制计算a*b%mod
llt multiMod(llt a,llt b,llt mod){
    llt ret = 0LL;
    a %= mod;
    while( b ){
        if ( b & 1LL ) ret = ( ret + a ) % mod, --b;
        b >>= 1LL;
        a = ( a + a ) % mod;
    }
    return ret;
}

//计算a^b%mod
llt powerMod(llt a,llt b,llt mod){
    llt ret = 1LL;
    a %= mod;
    while( b ){
        if ( b & 1LL ) ret = multiMod(ret,a,mod),--b;
        b >>= 1LL;
        a = multiMod(a,a,mod);
    }
    return ret;
}

//Miller-Rabin测试,测试n是否为素数
bool Miller_Rabin(llt n,int repeat){
    if ( 2LL == n || 3LL == n ) return true;
    if ( !( n & 1LL ) ) return false;

    //将n分解为2^s*d
    llt d = n - 1LL;
    int s = 0;
    while( !( d & 1LL ) ) ++s, d>>=1LL;

    srand((unsigned)time(0));
    for(int i=0;i<repeat;++i){//重复repeat次
        llt a = rand() % ( n - 3 ) + 2;//取一个随机数,[2,n-1)
        llt x = powerMod(a,d,n);
        llt y = 0LL;
        for(int j=0;j<s;++j){
            y = multiMod(x,x,n);
            if ( 1LL == y && 1LL != x && n-1LL != x ) return false;
            x = y;
        }
        if ( 1LL != y ) return false;
    }
    return true;
}