src/numtheory/cl_IF_millerrabin.cc

   1 // cl_miller_rabin_test().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_IF.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/modinteger.h"
  13
  14 namespace cln {
  15
  16 cl_boolean cl_miller_rabin_test (const cl_I& n, int count, cl_I* factor)
  17 {
  18         // [Cohen], section 8.2, algorithm 8.2.2.
  19         var cl_modint_ring R = find_modint_ring(n); // Z/nZ
  20         var cl_I m = n-1;
  21         var uintL e = ord2(m);
  22         m = m>>e;
  23         // n-1 = 2^e*m
  24         var cl_MI one = R->one();
  25         var cl_MI minusone = R->uminus(one);
  26         for (int i = 0; i < count; i++) {
  27                 // Choosing aa small makes the expt_pos faster.
  28                 var cl_I aa = (i == 0
  29                                ? (cl_I) 2
  30                                : i <= cl_small_prime_table_size
  31                                ? (cl_I) (unsigned int) cl_small_prime_table[i-1] // small prime
  32                                : 2+random_I(n-2)); // or random >=2, <n
  33                 if (aa >= n)
  34                         break;
  35                 // Now 1 < aa < n.
  36                 var cl_MI a = R->canonhom(aa);
  37                 var cl_MI b = R->expt_pos(a,m); // b = a^m
  38                 if (b == one)
  39                         goto passed;
  40                 for (uintL s = e; s > 0; s--) {
  41                         if (b == minusone)
  42                                 goto passed;
  43                         var cl_MI new_b = R->square(b);
  44                         if (new_b == one) {
  45                                 // (b-1)*(b+1) == 0 mod n, hence n not prime.
  46                                 if (factor)
  47                                         *factor = gcd(R->retract(b)-1,n);
  48                                 return cl_false;
  49                         }
  50                         b = new_b;
  51                 }
  52                 // b = a^(2^e*m) = a^(n-1), but b != -1 mod n, hence n not prime.
  53                 if (factor) {
  54                         var cl_I g = gcd(aa,n);
  55                         if (g > 1)
  56                                 *factor = g;
  57                         else
  58                                 *factor = 0;
  59                 }
  60                 return cl_false;
  61             passed:
  62                 ;
  63         }
  64         return cl_true;
  65 }
  66
  67 }  // namespace cln