From 674cd11ec63b56561980249cb19a0db54bacfa86 Mon Sep 17 00:00:00 2001
From: mdw <mdw>
Date: Mon, 13 Dec 1999 15:35:27 +0000
Subject: [PATCH] More changes.  Still embryonic.

---
 manual/catacomb.tex |   9 ++-
 manual/mp-mod.tex   | 156 ++++++++++++++++++++++++++++++++++++----------------
 manual/mp-mp.tex    | 144 ++++++++++++++++++++++++++++++++++++++++--------
 manual/mp-mpx.tex   |  57 ++++++++++---------
 4 files changed, 266 insertions(+), 100 deletions(-)

diff --git a/manual/catacomb.tex b/manual/catacomb.tex
index 2a2d3507..7cd84231 100644
--- a/manual/catacomb.tex
+++ b/manual/catacomb.tex
@@ -1,6 +1,6 @@
 %%% -*-latex-*-
 %%%
-%%% $Id: catacomb.tex,v 1.1 1999/12/10 23:27:11 mdw Exp $
+%%% $Id: catacomb.tex,v 1.2 1999/12/13 15:35:27 mdw Exp $
 %%%
 %%% Catacomb manual
 %%%
@@ -29,6 +29,9 @@
 %%%----- Revision history ---------------------------------------------------
 %%%
 %%% $Log: catacomb.tex,v $
+%%% Revision 1.2  1999/12/13 15:35:27  mdw
+%%% More changes.  Still embryonic.
+%%%
 %%% Revision 1.1  1999/12/10 23:27:11  mdw
 %%% Embryonic library reference manual.
 %%%
@@ -39,6 +42,7 @@
 \usepackage{mdwlist}
 \usepackage{mdwtab}
 \usepackage{mdwmath}
+\usepackage{mathenv}
 \usepackage{sverb}
 \usepackage{syntax}
 \usepackage{url}
@@ -71,10 +75,13 @@
   \par\nobreak\@nobreaktrue%
 }
 
+\show\tab@colstack
 \newcolumntype\!{!{\vline \qquad \qquad \vline}}
 
 \def\mp{\mathop{\operator@font MP}}
 
+\def\jacobi#1#2{{{#1}\overwithdelims()#2}}
+
 \def\mont#1{\tilde{#1}}
 
 \def\mm{\mathit{mm}}
diff --git a/manual/mp-mod.tex b/manual/mp-mod.tex
index 6c04ec36..bb3a9c77 100644
--- a/manual/mp-mod.tex
+++ b/manual/mp-mod.tex
@@ -12,8 +12,8 @@ performing modular arithmetic.
 \label{sec:mp-barrett}
 
 P.~Barrett invented an elegant method for computing modular reductions.  It
-requires a small amount of precomputation, depending only on the modulus $m$.
-Thereafter, it can compute residues mod $m$ using only multiplication and
+requires a small amount of precomputation, depending only on the modulus~$m$.
+Thereafter, it can compute residues mod~$m$ using only multiplication and
 arithmetic shifts.
 
 Montgomery's method (section~\ref{sec:mp-mont}, page~\pageref{sec:mp-mont})
@@ -26,7 +26,7 @@ The declarations required to use Barrett reduction are in
 initialized by calling \code{mpbarrett_create}
 (page~\pageref{fn:mpbarrett-create}) and disposed of by
 \code{mpbarrett_destroy} (page~\pageref{fn:mpbarrett-destroy}).  A number
-less then $m^2$ may then be reduced modulo $m$ by passing it to
+less then~$m^2$ may then be reduced modulo~$m$ by passing it to
 \code{mpbarrett_reduce} (page~\pageref{fn:mpbarrett-reduce}) along with the
 appropriate reduction context.
 
@@ -35,6 +35,78 @@ function \code{mpbarrett_exp} (page~\pageref{fn:mpbarrett-exp}) is provided
 which uses a simple square-and-multiply method combined with Barrett
 reduction.
 
+\subsubsection{How it works}
+
+For brevity, let $b = 2^{\code{MPW_BITS}}$ be the base we're working in.  Let
+$m$ be the modulus, and let $k$ be an integer such that $b^{k - 1} \le m <
+b^k$.
+
+The precomputation simply involves calculating $\mu = \lfloor b^{2k} / m
+\rfloor$.  Hence, we have the bound on~$\mu$
+\begin{equation}
+  \frac{b^{2k}}{m} - 1 < \mu \le \frac{b^{2k}}{m}.
+  \label{eq:mpbarrett-mu-bound}
+\end{equation}
+
+Let $x$ be an integer such that $0 \le x < b^{2k}$.  Firstly, let
+$q_0 = \lfloor x / b^{k - 1} \rfloor$.  Thus,
+\begin{equation}
+  \frac{x}{b^{k - 1}} - 1 < q_0 \le \frac{x}{b^{k - 1}}.
+  \label{eq:mpbarrett-q0-bound}
+\end{equation}
+This division an be done simply by ignoring the least significant $k - 1$
+words in the representation of~$x$, requiring no computation.
+
+Next, calculate $q = \lfloor q_0 \mu / b^{k + 1} \rfloor$.  The division can
+again be performed simply by ignoring the least significant $k + 1$ words of
+the result.
+
+Determining the bounds on $q$ is fairly easy, given
+equations \ref{eq:mpbarrett-mu-bound} and~\ref{eq:mpbarrett-q0-bound}.
+\[
+  \biggl\lfloor
+    \frac{(x / b^{k - 1} - 1)(b^{2k} \!/ m - 1) }{b^{k + 1}}
+  \biggr\rfloor
+  \le q \le
+  \biggl\lfloor \frac{(x / b^{k - 1})(b^{2k} \!/ m)}{b^{k + 1}} \biggr\rfloor
+\]
+Thus,
+\[
+  \biggl\lfloor
+   \frac{x}{m} - \frac{x}{b^{2k}} - \frac{b^{k - 1}}{m} + \frac{1}{b^{k + 1}}
+  \biggr\rfloor
+  \le q \le
+  \biggl\lfloor \frac{x}{m} \biggr\rfloor
+\]
+Since $b^{k - 1} \le m$ and $x < b^{2k}$, this can be reduced to simply
+\begin{equation}
+  \biggl\lfloor \frac{x}{m} \biggr\rfloor - 2
+  \le q \le
+  \biggl\lfloor \frac{x}{m} \biggr\rfloor
+  \label{eq:mpbarrett-q-bound}
+\end{equation}
+
+Finally, let $r = x - qm$.
+
+Firstly, because of the bounds on $q$ given in
+equation~\ref{eq:mpbarrett-q-bound}, we have
+\begin{equation}
+  x \bmod m \le r \le x \bmod m + 2m
+  \label{eq:mpbarrett-r-bound}
+\end{equation}
+Since $q$ is an integer, $r$ is equal to $x \bmod m + am$ where $a$ is $0$,
+$1$ or~$2$.
+
+Obviously, $r \le x \bmod m + 2m < 3m$.  If $b > 3$, then
+$3 m < 3 b^k < b^{k + 1}$, because $m < b^k$.  Hence the computation of $r$
+may be performed modulo $b^{k + 1}$, and in particular only the $k + 1$ least
+significant words of the product~$qm$ need be computed.
+
+The actual result $x \bmod m$ may be computed from $r$ by repeatedly
+subtracting $m$ while the result is still greater than $m$.  The bound in
+equation~\ref{eq:mpbarrett-r-bound} shows that this need be done at most
+twice.
+
 \subsubsection{The \code{mpbarrett_create} function}
 \label{fn:mpbarrett-create}
 
@@ -48,7 +120,7 @@ reduction.
 \fsec{Description}
 
 Initializes a Barrett reduction context $b$ suitable for performing
-reductions modulo $m$.  The memory used by the context block must be provided
+reductions modulo~$m$.  The memory used by the context block must be provided
 by the caller, and must not be discarded until the context is destroyed.
 
 The computation performed is as follows.  Let $k = 2 \cdot |MP_LEN(|m|)|$.
@@ -84,22 +156,12 @@ discarded.
 \fsec{Description}
 
 Calculates $x \bmod m$, where $m$ is the modulus configured into the Barrett
-reduction context $b$.  The argument $d$ is a suggested destination.
+reduction context~$b$.  The argument~$d$ is a suggested destination.
 
 Note that not all values of $x$ are suitable for reduction: in particular,
 negative values are improperly handled.  Let $k = 2 \cdot |MP_LEN(|m|)|$;
 then $x$ must be in the interval $[\,0, 2^{\code{MPW_BITS} \cdot k})$.  It's
-probably safer (and eqsier) to ensure that $x < m^2$.
-
-For reference, the calculation works as follows:
-\begin{itemize}
-\item To simplify the notation, let $b = 2^{\code{MPW_BITS}}$.
-\item Let $q_0 = \lfloor x / b^{k - 1} \rfloor$.  Let
-  $q = \lfloor q_0 \mu / b^{k + 1} \rfloor$.
-\item Let $r = (x - q m) \bmod b^{k + 1}$.
-\item At this point, $r \equiv x \pmod m$, but $r$ may be too large.  It can
-  be reduced by repeated subtraction until it's in range.
-\end{itemize}
+probably safer (and easier) to ensure that $x < m^2$.
 
 \fsec{Returns}
 
@@ -118,7 +180,7 @@ A value $d$ such that $0 \le d < m$ and $d \equiv x \pmod m$.
 \fsec{Description}
 
 Computes $x^e \bmod m$, where $m$ is the modulus configured into the Barrett
-reduction context $b$.  The argument $d$ is a suggested destination.
+reduction context~$b$.  The argument~$d$ is a suggested destination.
 
 \fsec{Returns}
 
@@ -133,8 +195,8 @@ reductions.  The method requires a number of values to be precomputed.
 \begin{itemize}
 \item Let $m$ be the modulus.  Let $b$ be the radix used to represent
   multiprecision integers.  (In Catacomb, then, $b = 2^{\code{MPW_BITS}}$.)
-  For Montgomery reduction to work, $m$ and $b$ must have no common factors.
-\item Let $R = b^n > m$ be the smallest power of $b$ which exceeds $m$.
+  For Montgomery reduction to work, $m$ and~$b$ must have no common factors.
+\item Let $R = b^n > m$ be the smallest power of $b$ which exceeds~$m$.
 \item Precompute $m'$ such that $m m' \equiv -1 \pmod b$.
 \item Precompute $R \bmod m$ and $R^2 \bmod m$, since these are generally
   useful when using Montgomery's method.
@@ -142,7 +204,7 @@ reductions.  The method requires a number of values to be precomputed.
 Given $x < R^2$ and the above values, Montgomery's algorithm computes $\mont
 x = x R^{-1} \bmod m$ efficiently, using only multiplications, additions and
 word-sized shifts.  The quantity $\mont x$ is called the \emph{Montgomery
-reduction of $x$}.
+reduction of~$x$}.
 
 At first sight, this isn't particularly useful.  However, an interesting
 thing happens if you run an extra factor of $R$ all the way through your
@@ -179,13 +241,12 @@ members:
 \begin{description}
   \def\makelabel#1{\code{#1}\hfil}
 \item[mp *m] The modulus $m$ with which the reduction context is working.
-\item[mp *mi] The quantity $m'$ where $m m' \equiv -1
-  \pmod{ 2^{\code{MPW_BITS}}}$.
-\item[size_t shift] The smallest integer $n$ such that
-  $R = 2^{\code{MPW_BITS} \cdot n} > m$.
-\item[mp *r] The quantity $R \bmod m$.  The Montgomery multiplicative
+\item[mp *mi] The integer $m'$ where $m m' \equiv -1 \pmod{R}$.
+\item[size_t n] The smallest integer $n$ such that $R = 2^{\code{MPW_BITS}
+    \cdot n} > m$.
+\item[mp *r] The integer $R \bmod m$.  The Montgomery multiplicative
   identity.
-\item[mp *r2] The quantity $R^2 \bmod m$.
+\item[mp *r2] The integer $R^2 \bmod m$.
 \end{description}
 All of these are set up by \code{mpmont_create} and should not be changed.
 
@@ -217,8 +278,8 @@ initialized.  Suppose $m$ is the modulus we need to work with:
 mpmont mm;
 mpmont_create(&mm, m);
 \end{listing}
-The next step is usually to multiply all of the inputs to the calculation by
-$R$.  This is usually best done like this:
+The next step is usually to multiply all of the inputs to the calculation
+by~$R$.  This is usually best done like this:
 \begin{listing}
 xr = mpmont_mul(&mm, MP_NEW, x, mm.r2);
 \end{listing}
@@ -279,17 +340,20 @@ page~\pageref{sec:mp-barrett}).
 
 \subsubsection{How it works}
 
-Let $x_0$ be the integer we wish to reduce.  Compute $u_1 = x_0 m' \bmod b$,
-and $v_1 = x_0 + u_1 m$.
+Let $x$ be the integer we wish to reduce.  Compute $u = x m' \bmod R$, and
+$v = x + u m$.
 
-Obviously $u_1 = x_0 m' + bk$ for some integer $k$.  So $v_1 = x_0 + u_1 m =
-x_0 + x_0 m m' + b k m$.  Note that $m m' \equiv -1 \pmod b$; thus $m m' =
-b k' - 1$ for some $k'$, so $v_1$ then becomes $x_0(1 + bk' - 1) + b k m =
-b (k' x_0 + k m)$.  Let $x_1 = v_1 / b = k' x_0 + k m$.
+Examining $v \bmod R$, we note that $v \equiv x + m m' x \equiv 0 \pmod R$,
+since $m m' \equiv -1 \pmod R$.  Thus, $v$ is a multiple of $R$.  Let $x' = v
+/ R$.
 
-Now note that $b k' = m m' + 1$, so $b k' \equiv 1 \pmod m$; hence, $k'
-\equiv b^{-1} \pmod m$.  Thus $x_1 \equiv x_0 b^{-1} \pmod m$.
+Examining $x' \bmod m$ is slightly trickier.  By definition, $u = x m' - k R$
+for some integer $k$.  Then $v = x + m(x m' - k R)$, and $x' \equiv v R^{-1}
+\equiv x R^{-1} \pmod m$.
 
+Finally, if $x < m R$ to start with, $u < R$ and so $v < 2 m R$.  Hence,
+$x' < 2 m$, so at most one subtraction of $m$ is enough to compute $x
+R^{-1} \bmod m$ from $x'$.
 
 
 \subsubsection{The \code{mpmont_create} function}
@@ -304,9 +368,9 @@ Now note that $b k' = m m' + 1$, so $b k' \equiv 1 \pmod m$; hence, $k'
 
 \fsec{Description}
 
-Initializes a Montgomery reduction context $\mm$, using the given modulus
-$m$.  The caller must provide memory for the context, and must ensure that
-the memory remains available until the context is destroyed.
+Initializes a Montgomery reduction context $\mm$, using the given
+modulus~$m$.  The caller must provide memory for the context, and must ensure
+that the memory remains available until the context is destroyed.
 
 If Catacomb is compiled with debugging support enabled, it will abort if $m$
 is negative or even.
@@ -338,14 +402,14 @@ had acquired.
 
 \fsec{Description}
 
-Performs a Montgomery reduction operation on the integer $x$, modulo $m$.
-The integer $d$ is a suggested destination.
+Performs a Montgomery reduction operation on the integer~$x$, modulo~$m$.
+The integer~$d$ is a suggested destination.
 
-For Montgomery reduction to work, $x$ must be less than $R^2$.  In practice,
+For Montgomery reduction to work, $x$ must be less than~$mR$.  In practice,
 it's probably safest to ensure that $x < m^2$.
 
-This function is particularly efficient if $d = x$, i.e., if you overwrite
-$x$ with its Montgomery reduction.
+This function can be particularly efficient if $d = x$, i.e., if you
+overwrite $x$ with its Montgomery reduction.
 
 \fsec{Returns}
 
@@ -363,10 +427,10 @@ An integer $d$ such that $0 \le d < m$ and $d R \equiv x \pmod m$.
 
 \fsec{Description}
 
-Multiplies $x$ and $y$, and returns the Montgomery reduction of the product
+Multiplies $x$ and~$y$, and returns the Montgomery reduction of the product
 modulo $m$.  The integer $d$ is a suggested destination.
 
-Both $x$ and $y$ must be less than $R$, and it's probably safer in practice
+Both $x$ and~$y$ must be less than $R$, and it's probably safer in practice
 if you ensure that $x, y < m$.
 
 \fsec{Returns}
diff --git a/manual/mp-mp.tex b/manual/mp-mp.tex
index 4f3aa670..45d829fa 100644
--- a/manual/mp-mp.tex
+++ b/manual/mp-mp.tex
@@ -64,8 +64,8 @@ more than once, which can otherwise lead to confusing bugs.
 
 \fsec{Returns}
 
-The number of words used in the representation of the integer $m$.  Note that
-the argument $m$ may be evaluated multiple times by this macro.
+The number of words used in the representation of the integer~$m$.  Note that
+the argument~$m$ may be evaluated multiple times by this macro.
 
 
 \subsubsection{The \code{mp_burn} function}
@@ -80,8 +80,8 @@ the argument $m$ may be evaluated multiple times by this macro.
 
 \fsec{Description}
 
-Sets the \code{MP_BURN} flag on the integer $m$.  Whenever memory associated
-with $m$, or results derived from it, is deallocated, it is overwritten with
+Sets the \code{MP_BURN} flag on the integer~$m$.  Whenever memory associated
+with~$m$, or results derived from it, is deallocated, it is overwritten with
 zeros to prevent traces being found in swap files or core dumps.
 
 
@@ -176,7 +176,7 @@ A pointer to a newly allocated multiprecision integer.
 
 \fsec{Description}
 
-Destroys the multiprecision integer $m$.  All resouces allocated to the
+Destroys the multiprecision integer~$m$.  All resouces allocated to the
 integer are freed.  It's most unusual for this to be the function you want:
 most of the time \code{mp_drop} (page~\pageref{fn:mp-drop}) is much more
 useful.
@@ -228,9 +228,9 @@ copying any data.\footnote{%
 
 \fsec{Description}
 
-The function \code{mp_copy} adds a reference to the multiprecision integer
-$m$.  The integer will not be destroyed until all references to it have been
-dropped.
+The function \code{mp_copy} adds a reference to the multiprecision
+integer~$m$.  The integer will not be destroyed until all references to it
+have been dropped.
 
 The macro \code{MP_COPY} performs exactly the same action as the
 \code{mp_copy} function, except that it uses inline code rather than a
@@ -239,7 +239,7 @@ function call.  It's probably smaller and faster than the function call too.
 \fsec{Returns}
 
 The address of the new reference.  This will usually be the same as the
-original argument $m$.
+original argument~$m$.
 
 \subsubsection{The \code{mp_drop} function and \code{MP_DROP} macro}
 \label{fn:mp-drop}
@@ -255,7 +255,7 @@ original argument $m$.
 \fsec{Description}
 
 The function \code{mp_drop} removes (`drops') a reference to the
-multiprecision integer $m$.  If there are no more references to the integer,
+multiprecision integer~$m$.  If there are no more references to the integer,
 it is destroyed, as if by \code{mp_destroy} (page~\pageref{fn:mp-destroy}).
 
 The macro \code{MP_DROP} performs exactly the same action as the
@@ -275,7 +275,7 @@ function call, so it's slightly faster but uses a little more code.
 
 \fsec{Description}
 
-If the integer $m$ has only one reference, the \code{mp_split} function
+If the integer~$m$ has only one reference, the \code{mp_split} function
 returns $m$ unchanged; otherwise, a copy is made and returned, and the
 reference count of $m$ is decremented.  Thus, either way, \code{mp_split}
 returns a pointer to an integer with the same value as $m$ and a single
@@ -284,7 +284,7 @@ refernce.  Thus, the result of \code{mp_split} is therefore safe to modify.
 The macro \code{MP_SPLIT} performs the same action as the \code{mp_split}
 function except that it uses inline code rather than a function call, and it
 modifies its argument in place rather than returning the new reference.  The
-macro \code{MP_SPLIT} may evaluate its argument $m$ multiple times.
+macro \code{MP_SPLIT} may evaluate its argument~$m$ multiple times.
 
 \fsec{Returns}
 
@@ -310,32 +310,32 @@ The actual behaviour is complicated:
 
 \begin{itemize}
 \item If $m$ is equal to \code{MP_NEW}, or the \code{MP_CONST} flag is set on
-  $m$, then allocate a new integer $d$ of size $\sz$ by calling
+  $m$, then allocate a new integer~$d$ of size $\sz$ by calling
   \code{mp_create} (page~\pageref{fn:mp-create}).
 \item Otherwise, if $m$ has more than one reference, a reference is removed
   asif by \code{mp_drop} (page~\pageref{fn:mp-drop}) and a new integer $d$ of
   size $\sz$ is allocated as if by \code{mp_create}.
-\item Otherwise, the integer $m$ is extended to $\sz$ words, as if by calling
-  \code{mp_ensure} (page~\pageref{fn:mp-ensure}), and $d$ is set equal to
-  $m$.
+\item Otherwise, the integer~$m$ is extended to $\sz$ words, as if by calling
+  \code{mp_ensure} (page~\pageref{fn:mp-ensure}), and $d$ is set equal
+  to~$m$.
 \end{itemize}
 
 The resulting integer $d$ is returned from the function.
 
 The upshot of all of this is that $d$ has exactly one reference, space for
-$\sz$ words of data, and no particular value set.  Any other references to
-the original integer $m$ (whether counted or not) continue to refer to the
-same value.  The total number of references to $d$ and $m$ after
+$\sz$~words of data, and no particular value set.  Any other references to
+the original integer~$m$ (whether counted or not) continue to refer to the
+same value.  The total number of references to $d$ and~$m$ after
 \code{mp_modify} returns is equal to the initial number of references to $m$
 before the call.
 
 The macro \code{MP_MODIFY} does the same job using inline code; it sets $m$
-to the new integer $d$ when it's finished.  The macro expands to a fairly
+to the new integer~$d$ when it's finished.  The macro expands to a fairly
 large chunk of code.
 
 \fsec{Returns}
 
-The function \code{mp_modify} returns the destination integer $d$ as
+The function \code{mp_modify} returns the destination integer~$d$ as
 described above.
 
 
@@ -359,9 +359,9 @@ of a multiprecision integer.
 \fsec{Description}
 
 The function \code{mp_resize} resizes the array used to hold the value of the
-integer $m$ such that it has space for $\sz$ words.  No check is made to see
+integer~$m$ such that it has space for $\sz$ words.  No check is made to see
 whether $m$'s array is already the correct size.  The integer's array limit
-pointer \code{vl} is set to point $\sz$ words higher than the base pointer
+pointer \code{vl} is set to point $\sz$~words higher than the base pointer
 \code{v}.
 
 The macro \code{MP_RESIZE} performs exactly the same action as the
@@ -427,7 +427,7 @@ The macro \code{MP_SHRINK} performs exactly the same operation as
 
 \fsec{Description}
 
-Minimizes the amount of memory used by the represenation of the integer $m$.
+Minimizes the amount of memory used by the represenation of the integer~$m$.
 This isn't usually useful unless $m$ is going to remain constant for a long
 time, or has become extremely large for some reason.
 
@@ -444,6 +444,102 @@ time, or has become extremely large for some reason.
 \subsection{More advanced algorithms}
 \label{sec:mp-adv}
 
+\subsubsection{The \code{mp_gcd} function}
+\label{fn:mp-gcd}
+
+\fsec{Synopsis}
+
+\begin{listinglist}
+|#include <catacomb/mp.h>| \\
+|void mp_gcd(mp **|$g$|, mp **|$a$|, mp **|$b$|, mp *|$x$|, mp *|$y$|);|
+\end{listinglist}
+
+\fsec{Description}
+
+The function \code{mp_gcd} implements an `extended GCD' algorithm.  Given
+integers $x$ and~$y$, it computes integers $g$, $a$, and~$b$ such that $g =
+\gcd(x, y)$ is the greatest common divisor of $x$ and $y$, and $ax + by = g$.
+
+If the arguments $g$, $a$ or~$b$ are null pointers, the appropriate result is
+discarded; otherwise, the current value of $|*|g$, $|*|a$ or $|*|b$ is used
+as a suggested destination for the result, and replaced by the actual result.
+
+If both $a$ and~$b$ are null, the results $a$ and $b$ are not computed.
+Hence \code{mp_gcd} can efficiently perform a simple greatest-common-divisor
+computation.
+
+A major use for \code{mp_gcd} is the calculation of modular inverses.  If
+$\gcd(x, n) = 1$ then \code{mp_gcd} computes integers $a$ and~$b$ such that
+$ax + bn = 1$, i.e., $ax \equiv 1 \pmod n$, so $a$ is the inverse of $x$,
+written $a \equiv x^{-1} \pmod n$.
+
+In order to better support this usage, \code{mp_gcd} implements a \emph{sign
+preference} system.  If only $a$ is requested, it will either be zero or
+have the same sign as~$y$; otherwise, $b$ will be zero or have the same sign
+as~$x$.  Thus, $x^{-1} \bmod n$ may be conveniently computed using the code
+\begin{listing}
+mp *xinv = MP_NEW;
+mp_gcd(0, 0, &xinv, n, x);
+\end{listing}
+(It's conventional to ensure that $x > y$, even though it doesn't actually
+make any difference.)
+
+The function \code{mp_gcd} works correctly for all integers $x$ and $y$.  The
+following properties of the results are guaranteed:
+\begin{itemize}
+  \unverb\|
+\item $\gcd(x, y) = \gcd(y, x)$ for all integers $x$ and~$y$.
+\item $\gcd(x, 0) = \gcd(0, x) = x$ for all integers~$x$.
+\item $|a| < |y|$, unless $y = 0$; similarly, $|b| < |x|$ if $x \ne 0$.
+\item If $x = y = 0$, then $a = b = 0$.  Otherwise, if $x = 0$, then $a = 0$
+  and $b = 1$; similarly, if $y = 0$, then $a = 1$ and $b = 0$.
+\end{itemize}
+
+\subsubsection{The \code{mp_jacobi} function}
+\label{fn:mp-jacobi}
+
+\fsec{Synopsis}
+
+\begin{listinglist}
+|#include <catacomb/mp.h>| \\
+|int mp_jacobi(mp *|$a$|, mp *|$n$|);|
+\end{listinglist}
+
+\fsec{Returns}
+
+The value of the  Jacobi symbol $\jacobi{a}{n}$.\footnote{%
+  Schneier writes the Jacobi symbol as $J(a, n)$.  The notation
+  $(\!\frac{a}{n}\!)$ seems more usual.} %
+The Jacobi symbol $\jacobi{a}{n}$ is defined for all integers~$a$ and odd
+integers~$n$ in terms of the Legendre symbol $\jacobi{a}{p}$ for prime~$p$,
+which is:
+\[
+  \jacobi{a}{p} = \begin{cases}
+     0 & if $a$ is a multiple of $p$ \\
+    -1 & if $a$ is a quadratic nonresidue mod $p$ \\
+     1 & if $a$ is a quadratic residue mod $p$
+   \end{cases}
+\]
+(An integer $x$ is a quadratic residue mod $n$ if and only if there exists an
+integer $y$ such that $y^2 \equiv x \pmod n$.  Note that $\jacobi{a}{p}
+\equiv a^{(p - 1)/2} \pmod p$.)
+
+The Jacobi symbol is defined as follows: let $p_0^{e_0} p_1^{e_1} \ldots
+p_k^{e_k} = n$ be the prime factors of $n$.  Then
+\[
+  \jacobi{a}{n} = \jacobi{a}{p_0}^{e_0} \jacobi{a}{p_1}^{e_1} \cdots
+    \jacobi{a}{p_k}^{e_k}
+\]
+Fortunately, the Jacobi symbol can be computed efficiently without factoring
+$n$ or performing modular exponentiations.
+
+The Jacobi symbol is used in the Solovay-Strassen primality test (not
+implemented in Catacomb -- the Rabin-Miller test is better), and in some
+subliminal channels in digital signature algorithms.
+
+The RSA public-key encryption algorithm leaks the Jacobi symbol of its
+plaintext.
+
 %%% Local Variables: 
 %%% mode: latex
 %%% TeX-master: "catacomb"
diff --git a/manual/mp-mpx.tex b/manual/mp-mpx.tex
index f775f6b6..28c14243 100644
--- a/manual/mp-mpx.tex
+++ b/manual/mp-mpx.tex
@@ -54,7 +54,7 @@ frequently used abbreviation for |(mpw)(|$x$| & MPW_MAX)|.
 
 \fsec{Returns}
 
-The number of bytes occupied by an array of $n$ \code{mpw}s (i.e., $n \cdot
+The number of bytes occupied by an array of $n$~\code{mpw}s (i.e., $n \cdot
 |sizeof(mpw)|$).
 
 \subsubsection{The \code{MPW_RQ} macro}
@@ -70,7 +70,7 @@ The number of bytes occupied by an array of $n$ \code{mpw}s (i.e., $n \cdot
 \fsec{Returns}
 
 The number of \code{mpw}s required to represent a multiprecision integer
-occupying $\sz$ octets in an external representation.
+occupying $\sz$~octets in an external representation.
 
 
 \subsection{Low-level multiprecision integer representation}
@@ -85,7 +85,7 @@ pointers to the array's \emph{base} (i.e., its first element) and
 \emph{limit} (i.e., the location immediately following its last element).
 
 Let $v$ be the base and $\vl$ the limit of a multiprecision integer array.
-The array itself is notated as $v .. \vl$.  The array's size in words may be
+The array itself is notated as~$v .. \vl$.  The array's size in words may be
 computed as $\vl - v$ in the normal C fashion.  The integer represented by
 the array, denoted $\mp(v .. \vl)$, is defined to be
 \[
@@ -100,8 +100,8 @@ the arithmetic algorithms don't need to consider high-order words which make
 no difference to the final result anyway.
 
 Whenever a result is too large to be represented in the memory allocated for
-it, high-order bits are discarded.  Thus, a result written to an array of $k$
-words is reduced modulo $2^{\code{MPW_BITS} \cdot k}$.
+it, high-order bits are discarded.  Thus, a result written to an array of
+$k$~words is reduced modulo $2^{\code{MPW_BITS} \cdot k}$.
 
 
 \subsection{Low-level macros}
@@ -122,9 +122,8 @@ manipulating multiprecision integers at the MPX level.
 
 \fsec{Description}
 
-The \code{MPX_SHRINK} macro reduces the limit $\vl$ of a multiprecision
-integer array so that either $v = \vl$ or $\vl[-1] \ne 0$.  The $\vl$
-argument must be an lvalue, since the macro writes the result back when it
+The \code{MPX_SHRINK} macro reduces the limit~$\vl$ of a multiprecision
+integer array so that either $v = \vl$ or~$\vl[-1] \ne 0$.  The argument~$\vl$ must be an lvalue, since the macro writes the result back when it
 finishes.
 
 \subsubsection{The \code{MPX_BITS} macro}
@@ -140,7 +139,7 @@ finishes.
 \fsec{Description}
 
 Determines the smallest number of bits which could represent the number
-$\mp(v .. \vl)$, and writes the answer to $b$, which must therefore be an
+$\mp(v .. \vl)$, and writes the answer to~$b$, which must therefore be an
 lvalue.  The result is zero if the number is zero; otherwise $b$ is the
 largest integer such that
 \[ \mp(v .. \vl) \ge 2^{(b - 1) \bmod \code{MPW_BITS}}.\]
@@ -159,11 +158,11 @@ largest integer such that
 \fsec{Description}
 
 Determines the smallest number of octets which could represent the number
-$\mp(v .. \vl)$, and writes the answer to $o$, which must therefore be an
+$\mp(v .. \vl)$, and writes the answer to~$o$, which must therefore be an
 lvalue.  This is useful for determining appropriate buffer sizes for the
 results of \code{mpx_storel} and \code{mpx_storeb}.
 
-The result $o$ can be calculated from the number of bits $b$ reported by
+The result $o$ can be calculated from the number of bits~$b$ reported by
 \code{MPX_BITS} as $o = \lceil b / 8 \rceil$; the algorithm used by
 \code{MPX_OCTETS} is more efficient than this, however.
 
@@ -240,8 +239,8 @@ in an octet array of a particular size.
 
 \fsec{Description}
 
-Stores the number held in the array $v .. \vl$ to the array of $\sz$ octets
-starting at address $p$ in little-endian byte order (i.e., least significant
+Stores the number held in the array $v .. \vl$ to the array of $\sz$~octets
+starting at address~$p$ in little-endian byte order (i.e., least significant
 byte first).
 
 \subsubsection{The \code{mpx_loadl} function}
@@ -257,8 +256,8 @@ byte first).
 
 \fsec{Description}
 
-Loads into the array $v .. \vl$ the number represented in the array of $\sz$
-octets starting at address $p$ in little-endian byte order (i.e., least
+Loads into the array $v .. \vl$ the number represented in the array of
+$\sz$~octets starting at address~$p$ in little-endian byte order (i.e., least
 significant byte first).
 
 \subsubsection{The \code{mpx_storeb} function}
@@ -274,8 +273,8 @@ significant byte first).
 
 \fsec{Description}
 
-Stores the number held in the array $v .. \vl$ to the array of $\sz$ octets
-starting at address $p$ in big-endian byte order (i.e., least significant
+Stores the number held in the array $v .. \vl$ to the array of $\sz$~octets
+starting at address~$p$ in big-endian byte order (i.e., least significant
 byte last).
 
 \subsubsection{The \code{mpx_loadb} function}
@@ -291,8 +290,8 @@ byte last).
 
 \fsec{Description}
 
-Loads into the array $v .. \vl$ the number represented in the array of $\sz$
-octets starting at address $p$ in big-endian byte order (i.e., least
+Loads into the array $v .. \vl$ the number represented in the array of
+$\sz$~octets starting at address $p$ in big-endian byte order (i.e., least
 significant byte last).
 
 
@@ -326,8 +325,8 @@ integer.
 
 \fsec{Description}
 
-Stores in $\dv .. \dvl$ the result of shifting $\mp(\av .. \avl)$ left by $n$
-bits (i.e., multiplying it by $2^n$).
+Stores in $\dv .. \dvl$ the result of shifting $\mp(\av .. \avl)$ left by
+$n$~bits (i.e., multiplying it by~$2^n$).
 
 \subsubsection{The \code{mpx_lsr} function}
 \label{fn:mpx-lsr}
@@ -345,10 +344,10 @@ bits (i.e., multiplying it by $2^n$).
 \fsec{Description}
 
 Stores in $\dv .. \dvl$ the result of shifting $\mp(\av .. \avl)$ right by
-$n$ bits (i.e., dividing it by $2^n$, rounding towards zero).
+$n$~bits (i.e., dividing it by~$2^n$, rounding towards zero).
 
 
-\subsection{Low level arithmetic}
+\subsection{Low-level arithmetic}
 
 The remaining functions perform standard arithmetic operations on large
 integers.  The form for the arguments is fairly well-standardized:
@@ -395,7 +394,7 @@ $\dv .. \dvl$.  The two arrays $v .. \vl$ and $\dv .. \dvl$ may be the same.
 \fsec{Description}
 
 The function \code{mpx_ucmp} performs an unsigned comparison of two unsigned
-multiprecision integers $a$ and $b$, passed in the arrays $\av .. \avl$ and
+multiprecision integers $a$ and~$b$, passed in the arrays $\av .. \avl$ and
 $\bv .. \bvl$ respectively.
 
 The macro \code{MPX_UCMP} provides a slightly more readable interface for
@@ -407,7 +406,7 @@ $a \mathrel{\synt{rel-op}} b$.
 
 The function \code{mpx_ucmp} returns a value less then, equal to, or
 greater than zero depending on whether $a$ is less than, equal to or greater
-than $b$.
+than~$b$.
 
 The macro \code{MPX_UCMP} returns a nonzero result if $a
 \mathrel{\synt{rel-op}} b$ is true, and zero if false.
@@ -447,7 +446,7 @@ destination array may be equal to either or both source arrays.\footnote{%
 
 \fsec{Description}
 
-The function \code{mpx_uaddn} adds a small integer $n$ (expressed as a single
+The function \code{mpx_uaddn} adds a small integer~$n$ (expressed as a single
 \code{mpw}) to the multiprecision integer held in $\dv .. \dvl$.
 
 The macro \code{MPX_UADDN} performs exactly the same operation, but uses
@@ -494,7 +493,7 @@ obtained.
 
 \fsec{Description}
 
-The function \code{mpx_usubn} subtracts a small integer $n$ (expressed as a
+The function \code{mpx_usubn} subtracts a small integer~$n$ (expressed as a
 single \code{mpw}) from the multiprecision integer held in $\dv .. \dvl$.
 
 The macro \code{MPX_USUBN} performs exactly the same operation, but uses
@@ -538,7 +537,7 @@ destination array may not be equal to either source array.
 \fsec{Description}
 
 The function \code{mpx_umuln} multiplies the multiprecision integer passed in
-$\av .. \avl$ by a small integer $n$ (expressed as a single \code{mpw}),
+$\av .. \avl$ by a small integer~$n$ (expressed as a single \code{mpw}),
 writing the product $n \cdot \mp(\av .. \avl)$ to the destination array $\dv
 .. \dvl$.  The destination array may be equal to the source array $\av
 .. \avl$.
@@ -564,7 +563,7 @@ inline code rather than calling a function.
 \fsec{Description}
 
 The function \code{mpx_umlan} multiplies the multiprecision integer passed in
-$\av .. \avl$ by a small integer $n$ (expressed as a single \code{mpw}), and
+$\av .. \avl$ by a small integer~$n$ (expressed as a single \code{mpw}), and
 adds it to the value already stored in the destination array $\dv .. \dvl$.
 The destination array may be equal to the source array $\av .. \avl$,
 although this isn't very useful.
-- 
2.11.0

