/cln-1.3.2/src/numtheory/cl_nt_jacobi_low.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/exception.h"
13#include "base/cl_xmacros.h"
14
15namespace cln {
16
17// Assume 0 <= a < b.
18inline 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
64int 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