素数判断算法 - 拉宾-米勒测试定理(c++实现)

在1000如此小的素数判断，在不考虑效率的情况下可以利用素数的定义来判断

printf("2 ");//2是唯一一个偶素数        
for( int a = 3; a <= 1000; a+=2) //步进为2, 因为只有奇数才有可能是素数（已排除了2）      
    {      
        bool bDivision = false;      
        int _nTemp = sqrt((float)a); //只需要检测自身开平方的数      
        for( int i = 3; i < nTemp; i+=2)      
        {      
            if( a % i == 0 )      
            {      
                bDivision = true;      
                break;      
            }      
        }      
     
        if( !bDivision )      
            printf("%d ", a );       
    }    

//========================================================

由于卡米歇尔数的存在，导致 费马小定理 无法判断一个数是否是素数。

费马小定理: 设p是素数, a是任意整数且 a!三0( mod p )， 则

                     a^(p-1)  三  1(mod p)

 //========================================================

卡米歇尔数:它是合数, 当 1<=a<=n, 都有 a^n 三 a(mod n)

 //========================================================

卡米歇尔数的考塞特判别法: 设n是合数，则n是卡米歇尔数当且仅当它是奇数，且整除n的每个素数p满足下述两个条件:

1)p^2 不整除 n

2)p-1 整除 n-1

 //========================================================                  

如果要判断相当大的素数最好使用 合数的拉宾-米勒测试定理

合数的拉宾-米勒测试定理: 设n是奇素数, 记 n-1 = 2^k * q , q 是奇数, 对不被n整除的某个a, 如果下述两个条件都成立，则n是合数.

a) a^q !三 1（mod n);

b) 对所有 i = 0, 1, 2, ...., k-1,    a^((2^i)*q) !三 -1(mod n);

//======================================================== 

这里给出了合数的拉宾-米勒测试定理  a  的取值:

if n < 1,373,653, it is enough to test a = 2 and 3.
if n < 9,080,191, it is enough to test a = 31 and 73.
if n < 4,759,123,141, it is enough to test a = 2, 7, and 61.
if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11.
//======================================================== 

// montgomery快速幂模算法 (n ^ p) % m, 与power算法极类似      
unsigned __int64 montgomery(unsigned __int64 n, unsigned __int64 p, unsigned __int64 m)      
{       
    unsigned __int64 r = n % m;      
    unsigned __int64 tmp = 1;      
    while (p > 1)      
    {      
        if ((p & 1)!=0)      
        {      
            tmp = (tmp * r) % m;      
        }      
        r = (r * r) % m;      
        p >>= 1;      
    }      
    return (r * tmp) % m;      
}      
     
//返回true:n是合数, 返回false:n是素数      
bool R_M_Help(unsigned __int64 a, unsigned __int64 k, unsigned __int64 q, unsigned __int64 n)      
{      
    if ( 1 != montgomery( a, q, n ) )      
    {      
        int e = 1;      
        for ( int i = 0; i < k; ++i )      
        {      
            if ( n - 1 == montgomery( a, q * e, n ) )       
                return false;      
                 
            e <<= 1;      
        }      
              
        return true;      
    }      
     
    return false;      
}      
     
//拉宾-米勒测试 返回true:n是合数, 返回false:n是素数        
bool R_M( unsigned __int64 n )       
{      
    if( n < 2 )      
        throw 0;      
     
    if ( n == 2 || n == 3 )      
    {      
        return false;      
    }      
     
    if( (n & 1) == 0 )      
        return true;      
     
    // 找到k和q, n = 2^k * q + 1;      
    unsigned __int64 k = 0, q = n - 1;      
    while( 0 == ( q & 1 ) )      
    {      
        q >>= 1;      
        k++;      
    }      
     
    /*if n < 1,373,653, it is enough to test a = 2 and 3.     
    if n < 9,080,191, it is enough to test a = 31 and 73.     
    if n < 4,759,123,141, it is enough to test a = 2, 7, and 61.     
    if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11.*/     
               
    if( n < 1373653 )      
    {      
        if( R_M_Help(2, k, q, n )       
         || R_M_Help(3, k, q, n ) )      
            return true;      
    }      
    else if( n < 9080191 )      
    {      
        if( R_M_Help(31, k, q, n )       
         || R_M_Help(73, k, q, n ) )      
            return true;      
    }        
    else if( n < 4759123141 )      
    {      
        if( R_M_Help(2, k, q, n )       
         || R_M_Help(3, k, q, n )      
         || R_M_Help(5, k, q, n )      
         || R_M_Help(11, k, q, n ) )      
            return true;      
    }      
    else if( n < 2152302898747 )      
    {      
        if( R_M_Help(2, k, q, n )       
         || R_M_Help(3, k, q, n )      
         || R_M_Help(5, k, q, n )      
         || R_M_Help(7, k, q, n )      
         || R_M_Help(11, k, q, n ) )      
            return true;      
    }      
    else       
    {      
        if( R_M_Help(2, k, q, n )       
         || R_M_Help(3, k, q, n )      
         || R_M_Help(5, k, q, n )      
         || R_M_Help(7, k, q, n )      
         || R_M_Help(11, k, q, n )      
         || R_M_Help(31, k, q, n )      
         || R_M_Help(61, k, q, n )      
         || R_M_Help(73, k, q, n ) )      
            return true;      
    }      
     
    return false;      
}