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
ready by tomorrow. ;-)
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
More information about the gmp-devel
mailing list