Miller-Rabin素数检测算法

 

Fermat定理 若n是奇素数，a是任意正整数(1≤ a≤ n−1)，则 aⁿ⁻¹≡1pmod n。[]

Miller-Rabin算法的理论基础 如果n是一个奇素数，将n−1表示成2^s*r的形式，r是奇数，a与n是互素的任何整数，那么a^r≡1pmod n或者对某个j(0 ≤ j≤ s−1, j∈Z)等式a^2jr≡−1 pmod n成立。[]

这个理论是由Fermat定理推导而来：n是一个奇素数，则方程 x²≡1pmod n只有±1两个解。

定理 3 设x，y和n是整数，如果x² = y²pmod n, 但x≠± ypmod n，则(x−y)和n的公约数中有n的非平凡因子。

输入：一个大于3的奇整数n和一个大于等于1的安全参数t（用于确定测试轮数）。

输出：返回n是否是素数（概率意义上的，一般误判概率小于((1/2))⁸⁰)即可)。

1. 将n-1表示成2^sr

2. 对i从1到t做循环做以下操作：

   1. 选择一个随机整数a（2 ≤ a ≤ n−2）

   2. 计算y ← a^r bmod n

   3. 如果y≠1并且y≠ n−1循环做下面的操作，否则转3：

      1. j ← 1

      2. 当j ≤ s−1并且y ≠ n−1循环做下面操作，否则跳到（iv.）

      3. 计算y ← y² bmod n，如果y = 1返回“合数”，否则j ← j + 1

      4. 如果y ≠ n−1则返回“合数”

3. 返回“素数”

其中第(vi.)步是基于定理3得来的。 []

经过独立的t轮Miller-Rabin算法将一个合数误判为素数的可能性不大于4^−t。这个概率是给予Fermat定理的算法中是最好的。这个误判概率是利用下面的两个定理证明的。

定理 1：设d = gcd(k,m)那么在有限群{g1,g + 2,···,gm = 1}中（g是有限群的生成元，m是有限群的阶）有d个元素满足方程x^k = 1。

定理 2： 设p是一个素奇数，p−1 = 2^sh（h是奇数），那么在乘法群(Z/pZ)^*中满足方程 x^2rt = −1pmod p（t是奇数）的元素个数是：0，如果r≥ s；2^rgcd(h,t)如果r<s。

利用这两个定理，对算法输入n分3种情况讨论：

1. n是可以被某个素数的平方整除的时候；

2. n是两个不同的素数的乘积的时候；

3. n是两个以上不同素数乘积的时候。

这样就可以证明Miiler-Rabin算法的误判概率上界。

Miller-Rabin算法是一个概率算法，算法的计算集中于(b)步和(c)步的循环中，最坏情况是(iv.)的循环没有中途推出，则一轮Miller-Rabin算法的最坏情况复杂度为(1 + O(1))log2(n)(以模n乘法为基本操作)。如果以单精度乘法操作作为时间复杂度的衡量，则一轮优化的Miller-Rabin算法的最坏情况时间复杂度是O(log2³(n))。从时间复杂度来看Miller-Rabin算法的性能是很好的。在实际应用中，Miller-Rabin算法的实际执行速度也很快。

代码实现：

#include <iostream>

using namespace std;

typedef unsigned __int64 llong;

llong mod_pro(llong x,llong y,llong n)
{
    llong ret=0,tmp=x%n;
    while(y)
    {
        if(y&0x1)if((ret+=tmp)>n)ret-=n;
        if((tmp<<=1)>n)tmp-=n;
        y>>=1;
    }
    return ret;
}
llong mod(llong a,llong b,llong c)
{
    llong ret=1;
    while(b)
    {
        if(b&0x1)ret=mod_pro(ret,a,c);
        a=mod_pro(a,a,c);
        b>>=1;
    }
    return ret;
}
llong ran()
{
    llong ret=rand();
    return ret*rand();
}
bool is_prime(llong n,int t)
{
    if(n<2)return false;
    if(n==2)return true;
    if(!(n&0x1))return false;
    llong k=0,m,a,i;
    for(m=n-1;!(m&1);m>>=1,k++);
    while(t--)
    {
        a=mod(ran()%(n-2)+2,m,n);
        if(a!=1)
        {
            for(i=0;i<k&&a!=n-1;i++)
                a=mod_pro(a,a,n);

            if(i>=k)return false;
        }
    }
    return true;
}
int main()
{
    llong n;
    while(scanf("%I64u",&n)!=EOF)
        if(is_prime(n,3))
            cout<<"YES\n";
        else
            cout<<"NO\n";
        return 0;
}