[Gmp-commit] /var/hg/gmp: Jacobi algorithm documentation, from Seth Troisi

[mercurial at gmplib.org]

details:   /var/hg/gmp/rev/578d5dbebf07
changeset: 17781:578d5dbebf07
user:      Niels M?ller <nisse at lysator.liu.se>
date:      Tue Jul 30 07:23:50 2019 +0200
description:
Jacobi algorithm documentation, from Seth Troisi

* doc/gmp.texi (Jacobi Symbol): Update algorithm documentation.
* tests/mpz/t-jac.c: Comment update.

diffstat:

 ChangeLog         |   6 ++++
 doc/gmp.texi      |  80 ++++++++++++++++++++++++++++++++++++++++++------------
 tests/mpz/t-jac.c |   3 --
 3 files changed, 68 insertions(+), 21 deletions(-)

diffs (125 lines):

diff -r 36b5993aefae -r 578d5dbebf07 ChangeLog
--- a/ChangeLog	Sat Jul 13 22:25:31 2019 +0200
+++ b/ChangeLog	Tue Jul 30 07:23:50 2019 +0200
@@ -1,3 +1,9 @@
+2019-07-30  Niels Möller  <nisse at lysator.liu.se>
+
+	From Seth Troisi:
+	* doc/gmp.texi (Jacobi Symbol): Update algorithm documentation.
+	* tests/mpz/t-jac.c: Comment update.
+
 2019-06-17  Torbjörn Granlund  <tg at gmplib.org>
 
 	* config.guess: Work around upstream configfsf.guess's regression wrt
diff -r 36b5993aefae -r 578d5dbebf07 doc/gmp.texi
--- a/doc/gmp.texi	Sat Jul 13 22:25:31 2019 +0200
+++ b/doc/gmp.texi	Tue Jul 30 07:23:50 2019 +0200
@@ -8934,23 +8934,64 @@
 @subsection Jacobi Symbol
 @cindex Jacobi symbol algorithm
 
-[This section is obsolete.  The current Jacobi code actually uses a very
-efficient algorithm.]
-
- at code{mpz_jacobi} and @code{mpz_kronecker} are currently implemented with a
-simple binary algorithm similar to that described for the GCDs (@pxref{Binary
-GCD}).  They're not very fast when both inputs are large.  Lehmer's multi-step
-improvement or a binary based multi-step algorithm is likely to be better.
-
-When one operand fits a single limb, and that includes @code{mpz_kronecker_ui}
-and friends, an initial reduction is done with either @code{mpn_mod_1} or
- at code{mpn_modexact_1_odd}, followed by the binary algorithm on a single limb.
-The binary algorithm is well suited to a single limb, and the whole
-calculation in this case is quite efficient.
-
-In all the routines sign changes for the result are accumulated using some bit
-twiddling, avoiding table lookups or conditional jumps.
-
+ at c Editor Note: I don't see other people defining the inputs, it would be nice
+ at c here because the code uses (a/b) where other references use (n/k)
+
+Jacobi symbol @m{\left(a \over b\right), (@var{a}/@var{b})}
+
+Initially if either operand fits in a single limb, a reduction is done with
+either @code{mpn_mod_1} or @code{mpn_modexact_1_odd}, followed by the binary
+algorithm on a single limb.  The binary algorithm is well suited to a single limb,
+and the whole calculation in this case is quite efficient.
+
+For inputs larger than @code{GCD_DC_THRESHOLD}, @code{mpz_jacobi},
+ at code{mpz_legendre} and @code{mpz_kronecker} are computed via the HGCD (Half
+GCD) function, as a generalization to Lehmer's algorithm.
+
+Most GCD algorithms reduce @math{a} and @math{b} by repeatatily computing the
+quotient @m{q = \lfloor a/b \rfloor, q = floor(a/b)} and iteratively replacing
+
+ at c Couldn't figure out macros with commas.
+ at tex
+$$ a, b = b, a - q * b$$
+ at end tex
+ at ifnottex
+ at math{a, b = b, a - q * b}
+ at end ifnottex
+
+Different algorithms use different methods for calculating q, but the core
+algorithm is the same if we use @ref{Lehmer's Algorithm} or
+ at ref{Subquadratic GCD, HGCD}.
+
+At each step it is possible to compute if the reduction inverts the Jacobi
+symbol based on the two least significant bits of @var{a} and @var{b}.  For
+more details see ``Efficient computation of the Jacobi symbol'' by
+M@"oller (@pxref{References}).
+
+A small set of bits is thus used to track state
+ at itemize
+ at item
+current sign of result (1 bit)
+
+ at item
+two least significant bits of @var{a} and @var{b} (4 bits)
+
+ at item
+a pointer to which input is currently the denominator (1 bit)
+ at end itemize
+
+In all the routines sign changes for the result are accumulated using fast bit
+twiddling which avoids conditional jumps.
+
+The final result is calculated after verifying the inputs are coprime (GCD = 1)
+by raising @m{(-1)^e,(-1)^e}
+
+Much of the HGCD code is shared directly with the HGCD implementations, such
+as the 2x2 matrix calculation, @xref{Lehmer's Algorithm} basecase and
+ at code{GCD_DC_THRESHOLD}.
+
+The asymptotic running time is @m{O(M(N)\log N),O(M(N)*log(N))}, where
+ at math{M(N)} is the time for multiplying two @math{N}-limb numbers.
 
 @need 1000
 @node Powering Algorithms, Root Extraction Algorithms, Greatest Common Divisor Algorithms, Algorithms
@@ -10910,9 +10951,12 @@
 Symposium on Computer Arithmetic, 1993, pp.@: 260 to 271.  Reprinted as ``More
 on Multiplying and Squaring Large Integers'', IEEE Transactions on Computers,
 volume 43, number 8, August 1994, pp.@: 899-908.
+
+ at item
+Niels M@"oller, ``Efficient computation of the Jacobi symbol'', @texlinebreak{}
+ at uref{https://arxiv.org/abs/1907.07795}
 @end itemize
 
-
 @node GNU Free Documentation License, Concept Index, References, Top
 @appendix GNU Free Documentation License
 @cindex GNU Free Documentation License
diff -r 36b5993aefae -r 578d5dbebf07 tests/mpz/t-jac.c
--- a/tests/mpz/t-jac.c	Sat Jul 13 22:25:31 2019 +0200
+++ b/tests/mpz/t-jac.c	Tue Jul 30 07:23:50 2019 +0200
@@ -25,9 +25,6 @@
 
 	   t-jac -p | gp -q
 
-   It takes a while because the output from "t-jac -p" is big.
-
-
    Enhancements:
 
    More big test cases than those given by check_squares_zi would be good.  */