added description of algorithm for mpfr_remainder

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@4445 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 47 insertions, 0 deletions

diff --git a/algorithms.tex b/algorithms.tex
index 8049c9da5..b522554e3 100644
--- a/algorithms.tex
+++ b/algorithms.tex
@@ -10,6 +10,7 @@
 \title{The MPFR\footnote{\lowercase{\url{http://www.mpfr.org/}}} Library:
 Algorithms and Proofs}
 \author{The MPFR team}
+
 \def\O{{\mathcal O}}
 \def\n{\textnormal}
 \def\pinf{\bigtriangleup}
@@ -31,6 +32,7 @@ Algorithms and Proofs}
 \def\to{{\bf to}}
 \def\while{{\bf while}}
 \def\err{{\rm err}}
+\def\cmod{\,\mathrm{cmod}\,}
 
 \newcommand{\U}[1]{\quad \mbox{[Rule~\ref{#1}]}}
 
@@ -811,6 +813,51 @@ Note that equality can occur --- i.e.\ the ``nearest even rounding rule'' ---
 only when the input has at least $2p+1$ bits; in particular it can not
 happen in the common case when input and output have the same precision.
 
+\subsection{The \texttt{mpfr\_remainder} and \texttt{mpfr\_remquo} functions}
+
+The \texttt{mpfr\_remainder} and \texttt{mpfr\_remquo} are two useful
+functions for argument reduction. Given two floating-point numbers $x$ and
+$y$, \texttt{mpfr\_remainder} computes the correct rounding of $x - q y$, where
+$q = \lfloor x/y \rceil$, with ties rounded to the nearest even integer,
+as in the rounding to nearest mode.
+Additionally, \texttt{mpfr\_remquo} returns $q \bmod 2^n$, where $n$
+is one less the machine word size ($n=31$ on a 32-bit machine, $n=63$ on a
+$64$-bit machine).
+
+Whatever the input $x$ and $y$, it should be noted that if $\ulp(x) \geq
+\ulp(y)$, then $x - q y$ is always
+exactly representable in the precision of $y$. To see this, 
+let $\ulp(y) = 2^{-k}$; 
+multiplying $x$ and $y$ by $2^k$ we get $X = 2^k x$ and 
+$Y = 2^k y$ such that $\ulp(Y)=1$,
+and $\ulp(X) \geq \ulp(Y)$, thus both $X$ and $Y$ are integers.
+Now perform the division of $X$ by $Y$, with quotient rounded to nearest:
+$X = q Y + R$, with $|R| \leq Y/2$. Since $R$ is an integer, it is
+necessarily representable with the precision of $Y$, and thus of $y$.
+The quotient $q$ of $x/y$ is the same as that of $X/Y$, the remainder
+$x - q y$ is $2^{-k} R$.
+
+We assume without loss of generality that $x, y > 0$, and that 
+$\ulp(y) = 1$, i.e., $y$ is an integer. 
+If $a$ and $b$ are two integers, we denote $a \cmod b$ the centered
+remainder of $a$ divided by $b$, with ties rounded to even quotient.
+\begin{quote}
+Algorithm Remainder. \\
+Input: $x$, $y$ with $\ulp(y)=1$, a rounding mode $\circ$ \\
+Output: $x \cmod y$, rounded according to $\circ$ \\
+1. If $\ulp(x) < 1$, decompose $x$ into $x_h + x_l$ with $\ulp(x_h) \geq 1$
+   and $0 \leq x_l < 1$. \\
+1a. $r \leftarrow \mathrm{Remainder}(x_h, y)$ [exact, $-y/2 \leq r \leq y/2$]\\
+1b. if $r < y/2$ or $x_l = 0$ then return $\circ(r + x_l)$ \\
+1c. else return $\circ(r + x_l - y) = \circ(x_l - r)$ \\
+2. Write $x = m \cdot 2^k$ with $k \geq 0$ \\
+3. $z \leftarrow 2^k \bmod y$ [binary exponentiation] \\
+4. Return $\circ(mz \cmod y)$.
+\end{quote}
+Note: at step (1a) the auxiliary variable $r$ has the precision of $y$;
+since $x_h$ and $y$ are integers, so it $r$ and the result is exact by the
+above reasoning.
+
 \section{High level functions}
 
 \subsection{The cosine function}