大数素性检验

大数素性检验：

#include <iostream>
#include <stdlib.h>
//#include "big.h"

//#pragma comment(lib, "miracl.lib")
using namespace std;

// Montgomery 快速幂模算法(n^p)%m
unsigned __int64 Montgomery(unsigned __int64 n, unsigned __int64 p, unsigned __int64 m)
{
    unsigned __int64 r = n % m;
    unsigned __int64 k = 1;
    while (p > 1)
    {
        if ((p & 1) != 0)
        {
            k = (k * r) % m; 
        }
        r = (r * r) % m;    
        p >>= 1;
    }
    return (r * k) % 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;
}

int main()
{
    ////初始化一个500位10进制的大数系统
    //Miracl precision(500, 10);
    //char a[512];
    //memset(a, 0, 512);

    //big p = mirvar(0); //大素数p
    //expb2(100, p);

    //while (1)
    //{
    //    // nxprime(pd,x)找到一个x大于pd的素数，返回值为BOOL
    //    nxprime(p, p);

    //    // 将一个大数根据其进制转换成一个字符串
    //    cotstr(p, a);
    //    if (isprime(p)) 
    //        cout << "p is a prime!" << "
" << a << "
";
    //    getchar();

    //}
    
    long long num;
    while (1)
    {
        scanf_s("%lld",&num);
        if (!R_M(num))
        {
            cout << "Is prime num." << endl;
        }
        else
        {
            cout << "No." << endl;
        }

    }

    return 0;
}