PageRenderTime 40ms CodeModel.GetById 16ms RepoModel.GetById 0ms app.codeStats 0ms

/cln-1.3.2/src/numtheory/cl_nt_jacobi_low.cc

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