# Jacobi symbol using Lehmer's algorithm.

Torbjorn Granlund tg at gmplib.org
Sun Jan 24 23:12:24 CET 2010

```nisse at lysator.liu.se (Niels Möller) writes:

> I'm looking into left-to-right Jacobi now, although I'm not sure how
> much time I will have for it. I think it should run at almost the same
> speed as gcd (either Lehmer or sub-quadratic). Just one additional table
> lookup per quotient.

Below is my current code. So far still quadratic (based on Lehmer's
algorithm). In my intitial benchmark, it seems to beat mpz_jacobi in
current GMP for n >= 3.

There seem to be a correlation between you saying "I don't have time for
X" and your immediate quality implementation of X.  :-)

The jacobi logic is implemented using a state of 6 bits (current
exponent for (-1)^e, two least significant bits of the current
remainders, and one bit saying which remainder we're currently
subtracting from the other.

Each time a multiple of one remainder in the sequence subtracted from
the other (that would be the quotient or some number smaller than the
quotient), the state is updated using a table lookup based on the
current state, the least two bits of the subtracted multiple/quotient,
and one bit saying which remainder is subtracted from which. I.e, the
table consists of 2^9 values of 6 bits each. I have special code for the
small cases n = 1 and 2 (which it might be better to rewrite using the
binary algorithm), and then the lehmer code is basically a copy of the
corresponding gcd code with appropriate calls to jacobi_update added,

There is some small cases code in the old jacobi code too.

Should be a fairly easy exercise to extend to a subquadratic algorithm,
using hgcd rather than hgcd2 for inputs above some threshold.

I hope you will also not have time for that, and expect the code to be

I timed your new code on shell.gmplib.org (Core i7, 2.67 GHz):

1   0.245   0.365   1.488
2   1.380   1.042   0.755
3   2.484   1.763   0.710
4   3.538   2.489   0.703
5   5.006   3.255   0.650
6   6.478   4.104   0.634
7   8.217   5.256   0.640
8   2.572   1.559   0.606
9  11.109   7.245   0.652
10   4.490   2.647   0.589
11   2.287   1.296   0.567
12  20.528   9.734   0.474
13   1.757   0.771   0.439
14  27.364  12.768   0.467
15  14.887   6.715   0.451
16  33.141  14.934   0.451
17  37.011  15.538   0.420
18  39.828  16.537   0.415
19   1.921   0.836   0.435
20  43.984  18.940   0.431
21  46.885  19.962   0.426
22  52.022  21.438   0.412
23  54.145  22.357   0.413
24  58.302  23.661   0.406
25  59.316  25.056   0.422
26   4.960   1.960   0.395
27  68.991  26.894   0.390
28  69.322  28.267   0.408
29  76.479  28.897   0.378
30  80.736  30.723   0.381
31  85.257  31.622   0.371
32  90.433  32.818   0.363
33  18.586   6.904   0.371
34  48.594  17.673   0.364
35  11.271   4.083   0.362

Really impressive speedups!  For huge operands, the new code is about
10x faster than the old code.

--
Torbjörn
```