src/numtheory/cl_nt_jacobi_low.cc

   1 // jacobi().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/numtheory.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/exception.h"
  13 #include "cl_xmacros.h"
  14
  15 namespace cln {
  16
  17 // Assume 0 <= a < b.
  18 inline int jacobi_aux (uintV a, uintV b)
  19 {
  20         var int v = 1;
  21         for (;;) {
  22                 // (a/b) * v is invariant.
  23                 if (b == 1)
  24                         // b=1 implies (a/b) = 1.
  25                         return v;
  26                 if (a == 0)
  27                         // b>1 and a=0 imply (a/b) = 0.
  28                         return 0;
  29                 if (a > (b >> 1)) {
  30                         // a > b/2, so (a/b) = (-1/b) * ((b-a)/b),
  31                         // and (-1/b) = -1 if b==3 mod 4.
  32                         a = b-a;
  33                         switch (b % 4) {
  34                                 case 1: break;
  35                                 case 3: v = -v; break;
  36                                 default: throw runtime_exception();
  37                         }
  38                         continue;
  39                 }
  40                 if ((a & 1) == 0) {
  41                         // b>1 and a=2a', so (a/b) = (2/b) * (a'/b),
  42                         // and (2/b) = -1 if b==3,5 mod 8.
  43                         a = a>>1;
  44                         switch (b % 8) {
  45                                 case 1: case 7: break;
  46                                 case 3: case 5: v = -v; break;
  47                                 default: throw runtime_exception();
  48                         }
  49                         continue;
  50                 }
  51                 // a and b odd, 0 < a < b/2 < b, so apply quadratic reciprocity
  52                 // law  (a/b) = (-1)^((a-1)/2)((b-1)/2) * (b/a).
  53                 if ((a & b & 3) == 3)
  54                         v = -v;
  55                 swap(uintV, a,b);
  56                 // Now a > 2*b, set a := a mod b.
  57                 if ((a >> 3) >= b)
  58                         a = a % b;
  59                 else
  60                         { a = a-b; do { a = a-b; } while (a >= b); }
  61         }
  62 }
  63
  64 int jacobi (sintV a, sintV b)
  65 {
  66         // Check b > 0, b odd.
  67         if (!(b > 0))
  68                 throw runtime_exception();
  69         if ((b & 1) == 0)
  70                 throw runtime_exception();
  71         // Ensure 0 <= a < b.
  72         if (a >= 0)
  73                 a = (uintV)a % (uintV)b;
  74         else
  75                 a = b-1-((uintV)(~a) % (uintV)b);
  76         return jacobi_aux(a,b);
  77 }
  78
  79 }  // namespace cln