Jacobi symbol using Lehmer's algorithm.

[Torbjorn Granlund]

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