src/numtheory/cl_nt_isprobprime.cc

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