大数素性检验

大数素性检验：

  1 #include <iostream>
  2 #include <stdlib.h>
  3 //#include "big.h"
  4
  5 //#pragma comment(lib, "miracl.lib")
  6 using namespace std;
  7
  8 // Montgomery 快速幂模算法(n^p)%m
  9 unsigned __int64 Montgomery(unsigned __int64 n, unsigned __int64 p, unsigned __int64 m)
 10 {
 11     unsigned __int64 r = n % m;
 12     unsigned __int64 k = 1;
 13     while (p > 1)
 14     {
 15         if ((p & 1) != 0)
 16         {
 17             k = (k * r) % m;
 18         }
 19         r = (r * r) % m;
 20         p >>= 1;
 21     }
 22     return (r * k) % m;
 23 }
 24
 25 // 返回true:n是合数, 返回false:n是素数
 26 bool R_M_Help(unsigned __int64 a, unsigned __int64 k, unsigned __int64 q, unsigned __int64 n)
 27 {
 28     if (1 != Montgomery(a, q, n))
 29     {
 30         int e = 1;
 31         for (int i = 0; i < k; ++i)
 32         {
 33             if (n - 1 == Montgomery(a, q * e, n))
 34                 return false;
 35             e <<= 1;
 36         }
 37         return true;
 38     }
 39     return false;
 40 }
 41
 42 //拉宾-米勒测试 返回true:n是合数, 返回false:n是素数
 43 bool R_M(unsigned __int64 n)
 44 {
 45     if (n < 2)
 46         throw 0;
 47     if (n == 2 || n == 3)
 48     {
 49         return false;
 50     }
 51     if ((n & 1) == 0)
 52         return true;
 53
 54     // 找到k和q, n = 2^k * q + 1;
 55     unsigned __int64 k = 0, q = n - 1;
 56     while (0 == (q & 1))
 57     {
 58         q >>= 1;
 59         k++;
 60     }
 61     /*if n < 1,373,653, it is enough to test a = 2 and 3.
 62     if n < 9,080,191, it is enough to test a = 31 and 73.
 63     if n < 4,759,123,141, it is enough to test a = 2, 7, and 61.
 64     if n < 2,152,302,898,747, it is enough to test a = 2, 3, 5, 7, and 11.*/
 65     if (n < 1373653)
 66     {
 67         if (R_M_Help(2, k, q, n) || R_M_Help(3, k, q, n))
 68             return true;
 69     }
 70     else if (n < 9080191)
 71     {
 72         if (R_M_Help(31, k, q, n) || R_M_Help(73, k, q, n))
 73             return true;
 74     }
 75     else if (n < 4759123141)
 76     {
 77         if (R_M_Help(2, k, q, n)
 78             || R_M_Help(3, k, q, n)
 79             || R_M_Help(5, k, q, n)
 80             || R_M_Help(11, k, q, n))
 81             return true;
 82     }
 83     else if (n < 2152302898747)
 84     {
 85         if (R_M_Help(2, k, q, n)
 86             || R_M_Help(3, k, q, n)
 87             || R_M_Help(5, k, q, n)
 88             || R_M_Help(7, k, q, n)
 89             || R_M_Help(11, k, q, n))
 90             return true;
 91     }
 92     else
 93     {
 94         if (R_M_Help(2, k, q, n)
 95             || R_M_Help(3, k, q, n)
 96             || R_M_Help(5, k, q, n)
 97             || R_M_Help(7, k, q, n)
 98             || R_M_Help(11, k, q, n)
 99             || R_M_Help(31, k, q, n)
100             || R_M_Help(61, k, q, n)
101             || R_M_Help(73, k, q, n))
102             return true;
103     }
104     return false;
105 }
106
107 int main()
108 {
109     ////初始化一个500位10进制的大数系统
110     //Miracl precision(500, 10);
111     //char a[512];
112     //memset(a, 0, 512);
113
114     //big p = mirvar(0); //大素数p
115     //expb2(100, p);
116
117     //while (1)
118     //{
119     //    // nxprime(pd,x)找到一个x大于pd的素数，返回值为BOOL
120     //    nxprime(p, p);
121
122     //    // 将一个大数根据其进制转换成一个字符串
123     //    cotstr(p, a);
124     //    if (isprime(p))
125     //        cout << "p is a prime!" << "\n" << a << "\n";
126     //    getchar();
127
128     //}
129
130     long long num;
131     while (1)
132     {
133         scanf_s("%lld",&num);
134         if (!R_M(num))
135         {
136             cout << "Is prime num." << endl;
137         }
138         else
139         {
140             cout << "No." << endl;
141         }
142
143     }
144
145     return 0;
146 }