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