More changes. Still embryonic.

[u/mdw/catacomb] / manual / mp-mod.tex

%%% -*-latex-*-

\section{Modular arithmetic}

Many public key cryptographic algorithms require modular arithmetic of
various kinds.  Modular exponentiation (calculation of $g^x \bmod n$) is
particularly important.  Catacomb provides two efficient methods for
performing modular arithmetic.


\subsection{Barrett reduction}
\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
arithmetic shifts.

Montgomery's method (section~\ref{sec:mp-mont}, page~\pageref{sec:mp-mont})
is similar, and more efficient for heavy use, but only works with odd moduli.
Barrett reduction works with an arbitrary modulus.

The declarations required to use Barrett reduction are in
\hdr{<catacomb/mpbarrett.h>}.  The precomuted values required are stored in a
\emph{Barrett reduction context} of type \code{mpbarrett}.  A context can be
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
\code{mpbarrett_reduce} (page~\pageref{fn:mpbarrett-reduce}) along with the
appropriate reduction context.

Since modular exponentiation is such a common operation in cryptography, a
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}

\fsec{Synopsis}

\begin{listinglist}
|#include <catacomb/mpbarrett.h>| \\
|void mpbarrett_create(mpbarrett *|$b$|, mp *|$m$|);|
\end{listinglist}

\fsec{Description}

Initializes a Barrett reduction context $b$ suitable for performing
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|)|$.
Calculate $\mu = \lfloor 2^{\code{MPW_BITS} \cdot k} / m \rfloor$.  The
values $k$, $\mu$ and $m$ are stored in the context block for later use.

\subsubsection{The \code{mpbarrett_destroy} function}
\label{fn:mpbarrett-destroy}

\fsec{Synopsis}

\begin{listinglist}
|#include <catacomb/mpbarrett.h>| \\
|void mpbarrett_destroy(mpbarrett *|$b$|);|
\end{listinglist}

\fsec{Description}

Disposes of a Barrett reduction context.  Resources associated with the
context block are freed.  The memory used to hold the context may safely be
discarded.

\subsubsection{The \code{mpbarrett_reduce} function}
\label{fn:mpbarrett-reduce}

\fsec{Synopsis}

\begin{listinglist}
|#include <catacomb/mpbarrett.h>| \\
|mp *mpbarrett_reduce(mpbarrett *|$b$|, mp *|$d$|, mp *|$x$|);|
\end{listinglist}

\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.

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 easier) to ensure that $x < m^2$.

\fsec{Returns}

A value $d$ such that $0 \le d < m$ and $d \equiv x \pmod m$.

\subsubsection{The \code{mpbarrett_exp} function}
\label{fn:mpbarrett-exp}

\fsec{Synopsis}

\begin{listinglist}
|#include <catacomb/mpbarrett.h>| \\
|mp *mpbarrett_exp(mpbarrett *|$b$|, mp *|$d$|, mp *|$x$|, mp *|$e$|);|
\end{listinglist}

\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.

\fsec{Returns}

A value $d$ such that $0 \le d < m$ and $d \equiv x^e \pmod m$.


\subsection{Montgomery reduction}
\label{sec:mp-mont}

Peter Montgomery has invented an ingenious method for computing modular
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$.
\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.
\end{itemize}
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$}.

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
calculation.  Note for example that the Montgomery reduction of $x R \cdot y
R$ is $xy R^2 R^{-1} \bmod m = x y R \bmod m$ -- there's still only a single
factor of $R$ in the result.  This can be removed finally by an additional
reduction step.

Catacomb provides a function for performing a simple Montgomery reduction,
and for calculating Montgomery multiplication: the Montgomery reduction of
the product of two integers.

Multiplying in the extra factor of $R$ can be done efficiently by multiplying
by $R^2 \bmod m$ and reducing.\footnote{%
  Some of the Catacomb comments refer to values with a factor of $R$ as
  having been `Montgomerized'.  While this is an ugly word, it does describe
  a useful concept.} %
This is the reason that $R^2 \bmod m$ is precomputed.  The value $R \bmod m$
is the Montgomery multiplicative identity and is often useful as an initial
value when forming products.

\subsubsection{Overview of functions provided}

All of the types and declarations needed for Montgomery reduction are in the
header file \hdr{<catacomb/mpmont.h>}.

Catacomb maintains precomputed values in a \emph{Montgomery reduction
context}, which has type \code{mpmont}.  A reduction context is initialized
by calling \code{mpmont_create} (page~\pageref{fn:mpmont-create}) and
disposed of by \code{mpmont_destroy} (page~\pageref{fn:mpmont-destroy}).

The Montgomery reduction context is a structure containing the following
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 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 integer $R^2 \bmod m$.
\end{description}
All of these are set up by \code{mpmont_create} and should not be changed.

Montgomery reduction can be performed by passing a Montgomery reduction
context and a number to be reduced to \code{mpmont_reduce}
(page~\pageref{fn:mpmont-reduce}).  Simultaneous multiplication and reduction
are performed by \code{mpmont_mul} (page~\pageref{fn:mpmont-mul}).

Catacomb can perform modular exponentiation using Montgomery reduction.  The
function \code{mpmont_exp} (page~\pageref{fn:mpmont-exp}) performs a standard
modular exponentiation; \code{mpmont_expr} (page~\pageref{fn:mpmont-expr})
does the same but doesn't remove the factor of $R$ from the result, which can
be useful if you want to perform further computations.

Catacomb can also perform multiple simultaneous modular exponentations,
multiplying the results together, without much of a performance penalty over
a single exponentation.  The function \code{mpmont_mexp}
(page~\pageref{fn:mpmont-mexp}) calculates the product $x_0^{e_0} x_1^{e_1}
\ldots x_{k - 1}^{e_{k - 1}} \bmod m$ given an array of $x_i$ and $e_i$;
\code{mpmont_mexpr} (page~\pageref{fn:mpmont-mexpr}) does the same but leaves
a factor of $R$ in the result.  These functions are particularly useful for
verifying discrete-logarith-based digital signatures.

\subsubsection{Calculation using Montgomery reduction}

Before any reduction work can be done, a Montgomery context must be
initialized.  Suppose $m$ is the modulus we need to work with:
\begin{listing}
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:
\begin{listing}
xr = mpmont_mul(&mm, MP_NEW, x, mm.r2);
\end{listing}
Because multiplication is distributive over addition, addition and
subtraction work normally, although it can be worth doing a little repeated
subtraction in case the result ends up larger than the modulus.
\begin{listing}
a = mp_add(a, a, b);
a = mp_add(a, a, c);
a = mp_add(a, a, d);
while (MP_CMP(a, >=, mm.m))
  a = mp_sub(a, a, mm.m);
\end{listing}
Multiplication of different numbers is best done with the supplied
multiply-and-reduce operation.
\begin{listing}
ab = mpmont_mul(&mm, MP_NEW, a, b);
\end{listing}
However, squaring is best done using separate square and reduce steps.
\begin{listing}
a2 = mp_sqr(MP_NEW, a);
a2 = mpmont_reduce(&mm, a2, a2);
\end{listing}
When the computation has finished, the results can be read out using
straightforward Montgomery reduction.  Don't forget to dispose of the
reduction context.
\begin{listing}
x = mpmont_reduce(&mm, x, xr);
mpmont_destroy(&mm);
\end{listing}

For particularly simple calculations, it can save a multiplication if you
reduce first and correct afterwards.
\begin{listing}
ab = mpmont_mul(&mm, MP_NEW, a, b);
ab = mpmont_mul(&mm, ab, ab, mm.r2);
\end{listing}
The first stage computes $ab R^{-1} \bmod m$.  The second computes $abR^{-1}
\cdot R^2 \cdot R^{-1} \bmod m = ab \bmod m$.  Doing this the usual way is
nowhere near as efficient:
\begin{listing}
ar = mpmont_mul(&mm, MP_NEW, a, mm.r2);
br = mpmont_mul(&mm, MP_NEW, b, mm.r2);
ab = mpmont_mul(&mm, ar, ar, br);
ab = mpmont_reduce(&mm, ab, ab);
mp_drop(br);
\end{listing}

Remember that for all of this magic to work, $m$ and $b$ must have no common
factors.  Since in Catacomb $b$ is a power of 2, this means simply that
\emph{$m$ must be odd}.  If you want to work with an even modulus, you'll
have to use Barrett reduction instead (section~\ref{sec:mp-barrett},
page~\pageref{sec:mp-barrett}).

\begin{note}[Important]
  Montgomery reduction only works when the modulus is an odd number.
\end{note}

\subsubsection{How it works}

Let $x$ be the integer we wish to reduce.  Compute $u = x m' \bmod R$, and
$v = x + u 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$.

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}
\label{fn:mpmont-create}

\fsec{Synopsis}

\begin{listinglist}
|#include <catacomb/mpmont.h>| \\
|void mpmont_create(mpmont *|$\mm$|, mp *|$m$|);|
\end{listinglist}

\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.

If Catacomb is compiled with debugging support enabled, it will abort if $m$
is negative or even.

\subsubsection{The \code{mpmont_destroy} function}
\label{fn:mpmont-destroy}

\fsec{Synopsis}

\begin{listinglist}
|#include <catacomb/mpmont.h>| \\
|void mpmont_destroy(mpmont *|$\mm$|);|
\end{listinglist}

\fsec{Description}

Destroys the Montgomery reduction context $\mm$, releasing any resources it
had acquired.

\subsubsection{The \code{mpmont_reduce} function}
\label{fn:mpmont-reduce}

\fsec{Synopsis}

\begin{listinglist}
|#include <catacomb/mpmont.h>| \\
|mp *mpmont_reduce(mpmont *|$\mm$|, mp *|$d$|, mp *|$x$|);|
\end{listinglist}

\fsec{Description}

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~$mR$.  In practice,
it's probably safest to ensure that $x < m^2$.

This function can be particularly efficient if $d = x$, i.e., if you
overwrite $x$ with its Montgomery reduction.

\fsec{Returns}

An integer $d$ such that $0 \le d < m$ and $d R \equiv x \pmod m$.

\subsubsection{The \code{mpmont_mul} function}
\label{fn:mpmont-mul}

\fsec{Synopsis}

\begin{listinglist}
|#include <catacomb/mpmont.h>| \\
|mp *mpmont_mul(mpmont *|$\mm$|, mp *|$d$|, mp *|$x$|, mp *|$y$|);|
\end{listinglist}

\fsec{Description}

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
if you ensure that $x, y < m$.

\fsec{Returns}

An integer $d$ such that $0 \le d < m$ and $d R \equiv x y \bmod m$.

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "catacomb"
%%% End: