src/numtheory/cl_nt_jacobi.cc

   1 // jacobi().
   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 "cln/integer.h"
  13 #include "integer/cl_I.h"
  14 #include "cln/exception.h"
  15 #include "base/cl_xmacros.h"
  16
  17 namespace cln {
  18
  19 int jacobi (const cl_I& a, const cl_I& b)
  20 {
  21         // Check b > 0, b odd.
  22         if (!(b > 0))
  23                 throw runtime_exception();
  24         if (!oddp(b))
  25                 throw runtime_exception();
  26  {      Mutable(cl_I,a);
  27         Mutable(cl_I,b);
  28         // Ensure 0 <= a < b.
  29         a = mod(a,b);
  30         // If a and b are fixnums, choose faster routine.
  31         if (fixnump(b))
  32                 return jacobi(FN_to_V(a),FN_to_V(b));
  33         var int v = 1;
  34         for (;;) {
  35                 // (a/b) * v is invariant.
  36                 if (b == 1)
  37                         // b=1 implies (a/b) = 1.
  38                         return v;
  39                 if (a == 0)
  40                         // b>1 and a=0 imply (a/b) = 0.
  41                         return 0;
  42                 if (a > (b >> 1)) {
  43                         // a > b/2, so (a/b) = (-1/b) * ((b-a)/b),
  44                         // and (-1/b) = -1 if b==3 mod 4.
  45                         a = b-a;
  46                         if (FN_to_V(logand(b,3)) == 3)
  47                                 v = -v;
  48                         continue;
  49                 }
  50                 if ((a & 1) == 0) {
  51                         // b>1 and a=2a', so (a/b) = (2/b) * (a'/b),
  52                         // and (2/b) = -1 if b==3,5 mod 8.
  53                         a = a>>1;
  54                         switch (FN_to_V(logand(b,7))) {
  55                                 case 3: case 5: v = -v; break;
  56                         }
  57                         continue;
  58                 }
  59                 // a and b odd, 0 < a < b/2 < b, so apply quadratic reciprocity
  60                 // law  (a/b) = (-1)^((a-1)/2)((b-1)/2) * (b/a).
  61                 if (FN_to_V(logand(logand(a,b),3)) == 3)
  62                         v = -v;
  63                 swap(cl_I, a,b);
  64                 // Now a > 2*b, set a := a mod b.
  65                 if ((a >> 3) >= b)
  66                         a = mod(a,b);
  67                 else
  68                         { a = a-b; do { a = a-b; } while (a >= b); }
  69         }
  70 }}
  71
  72 }  // namespace cln