src/numtheory/cl_nt_jacobi_low.cc

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