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