src/numtheory/cl_nt_isprobprime.cc

   1 // isprobprime().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/numtheory.h"
   8
   9
  10 // Implementation.
  11
  12 #include "numtheory/cl_IF.h"
  13 #include "cln/integer_io.h"
  14 #include "cln/exception.h"
  15 #include <sstream>
  16
  17 namespace cln {
  18
  19 bool isprobprime (const cl_I& n)
  20 {
  21         if (!(n > 0)) {
  22                 std::ostringstream buf;
  23                 fprint(buf, n);
  24                 fprint(buf, " is not a positive integer.");
  25                 throw runtime_exception(buf.str());
  26         }
  27         // With a Miller-Rabin count = 50 the final error probability is
  28         // 4^-50 < 10^-30.
  29         var int count = 50;
  30         // Step 1: Trial division (rules out 87% of all numbers quickly).
  31         const uint32 trialdivide_limit = 70;
  32         var uintC l = integer_length(n);
  33         if (l <= 32) {
  34                 var uint32 nn = cl_I_to_UL(n);
  35                 if (nn <= cl_small_prime_table_limit) {
  36                         // Table lookup.
  37                         var uintL i = cl_small_prime_table_search(nn);
  38                         if (i < cl_small_prime_table_size
  39                             && ((unsigned int) cl_small_prime_table[i] == nn
  40                                 || nn == 2))
  41                                 return true;
  42                         else
  43                                 return false;
  44                 }
  45                 if ((nn % 2) == 0 || cl_trialdivision(nn,1,trialdivide_limit))
  46                         return false;
  47                 // For small n, only few Miller-Rabin tests are needed.
  48                 if (nn < 2000U) count = 1; // {2}
  49                 else if (nn < 1300000U) count = 2; // {2,3}
  50                 else if (nn < 25000000U) count = 3; // {2,3,5}
  51                 else if (nn < 3200000000U) count = 4; // {2,3,5,7}
  52         } else if (l <= 64) {
  53                 var uint32 nhi = cl_I_to_UL(ldb(n,cl_byte(32,32)));
  54                 var uint32 nlo = cl_I_to_UL(ldb(n,cl_byte(32,0)));
  55                 if ((nlo % 2) == 0 || cl_trialdivision(nhi,nlo,1,trialdivide_limit))
  56                         return false;
  57         } else {
  58                 if (evenp(n) || cl_trialdivision(n,1,trialdivide_limit))
  59                         return false;
  60         }
  61         // Step 2: Miller-Rabin test.
  62         return cl_miller_rabin_test(n,count,NULL);
  63 }
  64
  65 }  // namespace cln