path: root/algorithms.tex

\documentclass[12pt]{amsart}
\usepackage{fullpage}
\pagestyle{empty}
\title{The MPFR Library: Algorithms and Proofs}
\author{The MPFR team}
\date{\tt www.mpfr.org}
\def\O{{\mathcal O}}
\def\q{\hspace*{5mm}}
\def\ulp{{\rm ulp}}
\begin{document}
\maketitle

\section{Useful formulas}

With a working precision of $n$ bits, and $x = m \cdot 2^e$
with $\frac{1}{2} \le m < 1$, $e := {\rm EXP}(x)$, we have:
\[ x \cdot 2^{-1-n} \le \ulp(x) = 2^{{\rm EXP}(x) - n} < x \cdot 2^{-n}. \]

We also have:
\[ \frac{1}{2} a \ulp(b) \le \ulp(a b) \le a \ulp(b). \]

\section{Low level functions}

\subsection{The {\tt mpfr\_add} function}

\begin{verbatim}
   mpfr_add (A, B, C, rnd)
   /* on suppose B et C de me^me signe, et EXP(B) >= EXP(C) */

   0. d = EXP(B) - EXP(C) /* d >= 0 par hypothe`se */
   1. Soient B1 les prec(A) premiers bits de B, et B0 le reste
	     C1 les bits de C correspondant a` B1, C0 le reste
   /* B0, C1, C0 peuvent e^tre vides, mais pas B1 */

	  <----------- A ---------->
	  <----------- B1 ---------><------ B0 ----->
	     <---------- C1 -------><------------ C0 ----------->

   2. A <- B1 + (C1 >> d)
   3. q <- compute_carry (B0, C0, rnd)
   4. A <- A + q
\end{verbatim}

\subsection{The {\tt mpfr\_cmp2} function}

This function computes the exponent shift when subtracting $c > 0$ from
$b \ge c$. In other terms, if ${\rm EXP}(x) := 
\lfloor \frac{\log b}{\log 2} \rfloor$, it returns:
it returns ${\rm EXP}(b) - {\rm EXP}(b-c)$.

This function admits the following specification in terms of the binary
representation of the mantissa of $b$ and $c$: if $b = u 1 0^n r$ and
$c = u 0 1^n s$, where $u$ is the longest common prefix to $b$ and $c$,
and $(r,s)$ do not start with $(0, 1)$, then ${\tt mpfr\_cmp2}(b,c)$ returns
$|u| + n$ if $r \ge s$, and $|u| + n + 1$ otherwise, where $|u|$ is the number
of bits of $u$.

As it is not very efficient to compare $b$ and $c$ bit-per-bit, we propose
the following algorithm, which compares $b$ and $c$ word-per-word.
Here $b[n]$ represents the $n$th word from the mantissa of $b$, starting from
the most significant word $b[0]$, which has its most significant bit set.
The values $c[n]$ represent the words of $c$, after a possible shift if the
exponent of $c$ is smaller than that of $b$.
\begin{verbatim}
   n = 0; res = 0;
   while (b[n] == c[n])
      n++;
      res += BITS_PER_MP_LIMB;

   /* now b[n] > c[n] and the first res bits coincide */

   dif = b[n] - c[n];
   while (dif == 1)
      n++;
      dif = (dif << BITS_PER_MP_LIMB) + b[n] - c[n];
      res += BITS_PER_MP_LIMB;

   /* now dif > 1 */

   res += equal_leading_bits(b[n], c[n]);

   if (!is_power_of_two(dif))
      return res;

   /* otherwise result is res + (low(b) < low(c)) */
   do
      n++;
   while (b[n] == c[n]);
   return res + (b[n] < c[n]);
\end{verbatim}

\section{Mathematical constants}

\subsection{Euler's constant}

Euler's constant is computed using the formula $\gamma = S(n) - R(n) - \log n$
where:
\[ S(n) = \sum_{k=1}^{\infty} \frac{n^k (-1)^{k-1}}{k! k}, \quad
   R(n) = \int_n^{\infty} \frac{\exp(-u)}{u} du \sim \frac{\exp(-n)}{n}
	\sum_{k=0}^{\infty} \frac{k!}{(-n)^k}. \]
This identity is attributed to Sweeney by Brent \cite{Brent78}.
We have $S(n) = _2 F_2(1,1;2,2;-n)$ and $R(n) = {\rm Ei}(1, n)$.

\paragraph{Evaluation of $S(n)$.}
As in \cite{Brent78}, let $\alpha \sim 4.319136566$ the positive root
of $\alpha + 2 = \alpha \log \alpha$, and $N = \lceil \alpha n \rceil$.
We approximate $S(n)$ by
\[ S'(n) = \sum_{k=1}^{N} \frac{n^k (-1)^{k-1}}{k! k}. \]
%         = \frac{1}{N!} \sum_{k=1}^N \frac{a_k}{k}, 
% where $a_k = n^k (-1)^{k-1} N!/k!$ is an integer.
% Therefore $a_k$ is exactly computed, and when dividing it by $k$
% (integer division)
% the error is at most $1$, and thus the absolute error on $S'(n)$ is
% at most $N/N! = 1/(N-1)!$.
The remainder term $S(n) - S'(n)$ is bounded by:
\[ |S(n) - S'(n)| \le \sum_{k=N+1}^{\infty} \frac{n^k}{k!}. \]
Since $k! \ge (k/e)^k$, and $k \ge N+1 \ge \alpha n$, we have:
\[ |S(n) - S'(n)| \le \sum_{k=N+1}^{\infty} \left( \frac{n e}{k} \right)^k
                  \le \sum_{k=N+1}^{\infty} \left( \frac{e}{\alpha} \right)^k
                  \le 2 \left( \frac{e}{\alpha} \right)^N
                  \le 2 e^{-2n} \]
since $(e/\alpha)^{\alpha} = e^{-2}$.

To approximate $S'(n)$, we use the following algorithm, where $m$ is the
working precision, and $a, s, t$ are integer variables:
\begin{quote}
$a \leftarrow 2^m$ \\
$s \leftarrow 0$ \\
{\bf for} $k$ {\bf from} $1$ {\bf to} $N$ {\bf do} \\
\q $a \leftarrow \lfloor \frac{n a}{k} \rfloor$ \\
\q $t \leftarrow \lfloor \frac{a}{k} \rfloor$ \\
\q $s \leftarrow s + (-1)^{k-1} t$ \\
return $x = s/2^m$
\end{quote}
The absolute error $\epsilon_k$ on $a$ at step $k$ satisfies
$\epsilon_k \le 1 + n/k \epsilon_{k-1}$ with $\epsilon_0=0$.
The maximum error is $\epsilon_n \le \frac{n^n}{n!} \le e^n$.
Thus the error on $t$ at step $k$ is less than $1 + e^n/k$,
and the total error on $s$ is bounded by $N (e^n + 1)$.
Hence to get a precision of $n$ bits, we need to use
$m ge n (1 + \frac{1}{\log 2})$.
In such a case, the value $s/2^m$ converted to a floating-point number
will have an error of at most $\ulp(x)$.

\paragraph{Evaluation of $R(n)$.}
We estimate $R(n)$ using the terms up to $k=n-2$, again 
as in \cite{Brent78}:
\[ R'(n) = \frac{e^{-n}}{n} \sum_{k=0}^{n-2} \frac{k!}{(-n)^k}. \]
% = \frac{\exp(-n)}{n^{n-1}} \sum_{k=0}^{n-2} (-1)^k \frac{k!} {n^{n-2-k}}.
% Here again, the integer sum is computed exactly, converting it to a 
% floating-point number introduces at most one ulp of error,
% $\exp(-n)$ is computed within one ulp,
% and $n^{n-1}$ within at most $n-2$ ulps.
% The division by $n^{n-1}$ and the multiplication introduce one more ulp of
% error, thus the total error on $R'(n)$ is at most $n+2$ ulps.
Let us introduce $I_k = \int_n^{\infty} \frac{e^{-u}}{u^k} du$.
We have $R(n) = I_1$ and the recurrence $I_k = \frac{e^{-n}}{n^k} - k I_{k+1}$,
which gives
\[ R(n) = \frac{e^{-n}}{n} \sum_{k=0}^{n-2} \frac{k!}{(-n)^k}
	+ (-1)^{n-1} (n-1)! I_n. \]
Bounding $n!$ by $(n/e)^n \sqrt{2 \pi (n+1)}$ which holds\footnote{
Formula 6.1.38 from \cite{AbSt73} gives
$x! = \sqrt{2\pi} x^{x+1/2} e^{-x+\frac{\theta}{12x}}$ for $x>0$ and
$0 < \theta < 1$.
Using it for $x \ge 1$, we have $0 < \frac{\theta}{6x} < 1$, and
$e^t < 1+2t$ for $0 < t < 1$, thus 
$(x!)^2 \le 2\pi x^{2x} e^{-2x} (x+\frac{1}{3})$.}
for $n \ge 1$, we have:
\[ |R(n) - R'(n)| = (n-1)! I_n 
	\le \frac{n!}{n} \int_n^{\infty} \frac{e^{-n}}{u^n} du
	\le \frac{n^{n-1}}{e^n} \sqrt{2 \pi (n+1)} \frac{e^{-n}}{(n-1) 
	n^{n-1}} \]
and since $\sqrt{2 \pi (n+1)}/(n-1) \le 1$ for $n \ge 9$:
\[ |R(n) - R'(n)| \le e^{-2n} \quad \mbox{for $n \ge 9$}. \]
Thus we have:
\[ |\gamma - S'(n) - R'(n) - \log n| \le 3 e^{-2n} \quad \mbox{for $n\ge 9$}.\]
% If the working precision is $p$, then choose $n \ge \frac{\log 2}{2} (p+2)$
% such that $3 e^{-2n}$ represents at most one ulp.

To approximate $R'(n)$, we use the following:
\begin{quote}
$m \leftarrow {\rm prec}(x) - \lfloor \frac{n}{\log 2} \rfloor$ \\
$a \leftarrow 2^m$ \\
$s \leftarrow 1$ \\
{\bf for} $k$ {\bf from} $1$ {\bf to} $n$ {\bf do} \\
\q $a \leftarrow \lfloor \frac{k a}{n} \rfloor$ \\
\q $s \leftarrow s + (-1)^{k} a$ \\
$t \leftarrow \lfloor s/n \rfloor$ \\
$x \leftarrow t/2^m$ \\
return $r = e^{-n} x$
\end{quote}
The absolute error $\epsilon_k$ on $a$ at step $k$ satisfies
$\epsilon_k \le 1 + k/n \epsilon_{k-1}$ with $\epsilon_0=0$.
As $k/n \le 1$, we have $\epsilon_k \le k$, whence the error
on $s$ is bounded by $n(n+1)/2$, and that on $t$ by 
$1 + (n+1)/2 \le n+1$ since $n \ge 1$.
The operation $x \leftarrow t/2^m$ is exact as soon as ${\rm prec}(x)$ is large
enough, thus the error on $x$ is at most $(n+1) \frac{e^n}{2^{{\rm prec}(x)}}$.
If $e^{-n}$ is computed with $m$ bits, then
the error on it is at most $e^{-n} 2^{1-m}$.
The error on $r$ is then $(n + 1 + 2/n) 2^{-{\rm prec}(x)} +
\ulp(r)$.
Since $x \ge \frac{2}{3} n$ for $n \ge 2$, and $x 2^{-{\rm prec}(x)}
< \ulp(x)$, this gives an error bounded by 
$\ulp(r) + (n + 1 + 2/n) \frac{3}{2n} \ulp(r)
\le 4 \ulp(r)$ for $n \ge 2$ --- if ${\rm prec}(x) = {\rm prec}(r)$.
Now since $r \le \frac{e^{-n}}{n} \le \frac{1}{8}$, that error
is less than $\ulp(1/2)$.

\paragraph{Final computation.} We use the formula
$\gamma = S(n) - R(n) - \log n$ with $n$ such that $e^{-2n} \le 
\ulp(1/2) = 2^{-{\rm prec}}$, i.e.~$n \ge {\rm prec} \frac{\log 2}{2}$:
\begin{quote}
$s \leftarrow S'(n)$ \\
$l \leftarrow \log(n)$ \\
$r \leftarrow R'(n)$ \\
return $(s - l) - r$
\end{quote}
Since the final result is in $[\frac{1}{2}, 1]$, and $R'(n) \le 
\frac{e^{-n}}{n}$, then $S'(n)$ approximates $\log n$.
If we use $m + \lceil \log_2({\rm prec}) \rceil$ bits to evaluate $s$ and $l$,
then the error on $s-l$ will be at most $3 \ulp(1/2)$,
so the error on $(s - l) - r$ is at most $5 \ulp(1/2)$,
and adding the $3 e^{-2n}$ truncation error, we get a bound of
$8 \ulp(1/2)$.

\subsubsection{A faster formula}

Brent and McMillan give in \cite{BrMc80} a faster algorithm (B2) using the
modified Bessel functions, which was
used by Gourdon and Demichel to compute 108,000,000 digits of $\gamma$ in
October 1999:
\[ \gamma = \frac{S_0 - K_0}{I_0} - \log n, \]
where $S_0 = \sum_{k=1}^{\infty} \frac{n^{2k}}{(k!)^2} H_k$,
$H_k = 1 + \frac{1}{2} + \cdots + \frac{1}{k}$ being the $k$-th harmonic
number,
$K_0 = \sqrt{\frac{\pi}{4n}} e^{-2n} \sum_{k=0}^{\infty}
	(-1)^k \frac{[(2k)!]^2}{(k!)^3 (64n)^k}$,
and $I_0 = \sum_{k=0}^{\infty} \frac{n^{2k}}{(k!)^2}$.

We have $I_0 \ge \frac{e^{2n}}{\sqrt{4 \pi n}}$ (see \cite{BrMc80} page 306).
From the remark following formula 9.7.2 of \cite{AbSt73},
the remainder in the truncated expansion for $K_0$ up to $k$ does not
exceed the $(k+1)$-th term, whence
$K_0 \le \sqrt{\frac{\pi}{4n}} e^{-2n}$ and
$\frac{K_0}{I_0} \le \pi e^{-4n}$ as in formula (5) of \cite{BrMc80}.
Let $I'_0 = \sum_{k=0}^{K-1} \frac{n^{2k}}{(k!)^2}$ with
$K = \lceil \beta n \rceil$, and $\beta$ is the root of
$\beta (\log \beta - 1) = 3$
($\beta \sim 4.971...$) then
\[ |I_0 - I'_0| \le \frac{\beta}{2 \pi (\beta^2-1)} \frac{e^{-6n}}{n}. \]
Let $K'_0 = \sqrt{\frac{\pi}{4n}} e^{-2n} \sum_{k=0}^{4n-1} (-1)^k
	\frac{[(2k)!]^2}{(k!)^3 (64n)^k}$, then bounding by the next term:
\[ |K_0 - K'_0| \le \frac{(8n+1)}{16 \sqrt{2} n} \frac{e^{-6n}}{n}
	\le \frac{1}{2} \frac{e^{-6n}}{n}. \]
We get from this
\[ \left| \frac{K_0}{I_0} - \frac{K'_0}{I'_0} \right| 
	\le \frac{1}{2 I_0} \frac{e^{-6n}}{n} \le \sqrt{\frac{\pi}{n}}
	e^{-8n}. \]
Let $S'_0 = \sum_{k=1}^{K-1} \frac{n^{2k}}{(k!)^2} H_k$,
then using $\frac{H_{k+1}}{H_k} \le \frac{k+1}{k}$ and the same bound $K$
than for $I'_0$ ($4n \le K \le 5n$), we get:
\[ |S_0 - S'_0| \le \frac{\beta}{2 \pi (\beta^2-1)} H_k \frac{e^{-6n}}{n}. \]
We deduce:
\[ \left| \frac{S_0}{I_0} - \frac{S'_0}{I'_0} \right|
	\le e^{-8n} H_K \frac{\sqrt{4 \pi n}}{\pi (\beta^2-1)}
	\frac{\beta}{n} \le e^{-8n}. \]
Hence we have
\[ \left| \gamma - \left( \frac{S'_0 - K'_0}{I'_0} - \log n \right) \right|
	\le (1 + \sqrt{\frac{\pi}{n}}) e^{-8n} 
	\le 3 e^{-8n}. \]

\section{High level functions}

\subsection{The cosine function}

To evaluate $\cos x$ with a target precision of $n$ bits, we use the following
algorithm with working precision $m$:
\begin{quote}
$k \leftarrow \lfloor \sqrt{n/2} \rfloor$ \\
$r \leftarrow x^2$ rounded up \\ % err <= ulp(r)
$r \leftarrow r/2^{2k}$ \\ % err <= ulp(r)
$s \leftarrow 1, t \leftarrow 1$ \\ % err = 0
{\bf for} $l$ {\bf from} $1$ {\bf while} ${\rm EXP}(t) \ge -m$ \\
\q $t \leftarrow t \cdot r$ rounded up \\ % err <= (3*l-1)*ulp(t)
\q $t \leftarrow \frac{t}{(2l-1)(2l)}$ rounded up \\ % err <= 3*l*ulp(t)
\q $s \leftarrow s + (-1)^l t$ rounded down\\ % err <= l/2^m
{\bf do} $k$ times \\
\q $s \leftarrow 2 s^2$ rounded up \\
\q $s \leftarrow s - 1$ \\
return $s$ \\
\end{quote}
The error on $r$ after $r \leftarrow x^2$
is at most $1 \ulp(r)$ and remains $1 \ulp(r)$ after
$r \leftarrow r/2^{2k}$ since that division is just an exponent shift.
By induction, the error on $t$ after step $l$ of the for-loop is at most
$3 l \ulp(t)$.
Hence as long as $3 l \ulp(t)$ remains less than $\le 2^{-m}$
during that loop
(this is possible as soon as $r < 1/\sqrt{2}$)
and the loop goes to $l_0$, the error on $s$ after the for-loop is at most
$2 l_0 2^{-m}$ (for $|r| < 1$, it is easy to check that $s$ will remain
in the interval $[\frac{1}{2}, 1[$, thus $\ulp(s) = 2^{-m}$).
(An additional $2^{-m}$ term represents the truncation error,
but for $l=1$ the value of $t$ is exact, giving $(2 l_0 - 1) + 1 = 2 l_0$.)

Denoting by $\epsilon_i$ the maximal error on $s$ after the $i$th step
in the do-loop, we have $\epsilon_0 = 2 l_0 2^{-m}$ and
$\epsilon_{k+1} \le 4 \epsilon_k + 2^{-m}$,
giving $\epsilon_k \le (2 l_0+1/3) 2^{2k-m}$.

\subsection{The sine function}

The sine function is computed from the cosine, with a working precision of
$m$ bits:
\begin{quote}
$c \leftarrow \cos x$ rounded to zero \\
$t \leftarrow c^2$ rounded up \\
$u \leftarrow 1 - t$ rounded to nearest \\
$s \leftarrow {\rm sign}(x) \sqrt{u}$ rounded to nearest \\
\end{quote}
The absolute error on $c$ is at most $\ulp(c)
\le 2^{-m}$ since $|c| < 1$,
then that on
$t$ is at most $3 \cdot 2^{-m}$, that on $u$ is at most $3 \cdot 2^{-m} + 
\frac{1}{2}\ulp(u) \le \frac{7}{2} \cdot 2^{-m}$
since $|t|, |u| < 1$, so that on $s$ is at most ${\rm error}(\sqrt{u})
 + \frac{1}{2} \ulp(s) \le 2^{-m} (1/2 + \frac{7}{4 \sqrt{u}})$.
(By Rolle's theorem,
$|\sqrt{u} - \sqrt{u'} \le \frac{1}{2 \sqrt{v}} |u-u'|$ for
$v \in [u, u']$.)
% |s'-s| <= ulp(s) + |sqrt(u)-sqrt(u')| <= ulp(s) + (u-u')/(2*sqrt(v))
% for v in [u,u']
Let $u = m 2^{-2e}$ with $1 \le m < 4$ and $e \ge 1$, then 
the error on $s$ is at most $2^{-m} (1/2 + \frac{7 \cdot 2^{e}}{4})
	\le 2^{e + 1 - m}$.

\subsection{The tangent function}

The tangent function is computed from the cosine, 
using $\tan x = {\rm sign}(x) \sqrt{\frac{1}{\cos^2 x} - 1}$,
with a working precision of $m$ bits:
\begin{quote}
$c \leftarrow \cos x$ rounded down \\ % c <= cos(x) <= 1
$t \leftarrow c^2$ rounded down \\    % t <= cos(x)^2 <= 1
$v \leftarrow 1/t$ rounded up \\      % v >= 1/cos(x)^2 >= 1
$u \leftarrow v - 1$ rounded up \\    % u >= 1/cos(x)^2 - 1
$s \leftarrow {\rm sign}(x) \sqrt{u}$ rounded away from $0$ \\
\end{quote}
The absolute error on $c$ is at most $\ulp(c)$.
Hence the error on
$t$ is at most $\ulp(t) + 2 c \ulp(c) \le 5 \ulp(t)$,
% err(t) <= ulp(t) + |c^2-c'^2| <= ulp(t) + |c-c'|*|c+c'|
%        <= ulp(t) + 2*ulp(c)*c <= 5*ulp(t)
that on $v$ is at most $\ulp(v) + 5 \ulp(t)/t^2
	\le \ulp(v) + 10 \ulp(1/t) \le 11 \ulp(v)$,
% err(v) <= ulp(v) + |1/t-1/t'| <= ulp(v) + |t-t'|/t/t'
%        <= ulp(v) + 5*ulp(t)/t^2 <= ulp(v) + 10*ulp(1/t) <= 11*ulp(v)
that on $u$ is at most $\ulp(u) + 11 \ulp(v) \le (1 + 11 \cdot 2^e) \ulp(u)$
where $e$ is the exponent difference between $v$ and $u$.
% err(u) <= ulp(u) + err(v) <= ulp(v) + 11*ulp(v) <= (1+11*2^e)*ulp(u)
The final error on $s$ is $\le \ulp(s) + (1+11 \cdot 2^e) \ulp(u)/2/\sqrt{u}
	\le \ulp(s) + (1+11 \cdot 2^e) \ulp(u/\sqrt{u})
	\le (2 + 11 \cdot 2^e) \ulp(s)$.
% err(s) <= ulp(s) + |u-u'|/2/sqrt(u) <= ulp(s) + (1+11*2^e)*ulp(u)/2/sqrt(u)
%        <= ulp(s) + (1+11*2^e)*ulp(u/sqrt(u))
%        <= (2+11*2^e)*ulp(s) 

\subsection{The exponential function}

The {\tt mpfr\_exp} function implements three different algorithms.
For very large precision, it uses a $\O(M(n) \log^2 n)$ algorithm
based on binary splitting, based on the generic implementation for
hypergeometric functions in the file {\tt generic.c} (see \cite{Jeandel00}).
Currently this third algorithm is used only for precision greater
than $13000$ bits.

For smaller precisions, it uses Brent's method~;
if $r = (x - n \log 2)/2^k$ where $0 \le r < \log 2$, then 
\[ \exp(x) = 2^n \cdot \exp(r)^{2^k} \]
and $\exp(r)$ is computed using the Taylor expansion:
\[ \exp(r) =  1 + r + \frac{r^2}{2!} + \frac{r^3}{3!} + \cdots \]
As $r < 2^{-k}$, if the target precision is $n$ bits, then only
about $l = n/k$ terms of the Taylor expansion are needed.
This method thus requires the evaluation of the Taylor series to
order $n/k$, and $k$ squares to compute $\exp(r)^{2^k}$.
If the Taylor series is evaluated using a na\"{\i}ve way, the optimal
value of $k$ is about $n^{1/2}$, giving a complexity of $\O(n^{1/2} M(n))$.
This is what is implemented in {\tt mpfr\_exp2\_aux}.

If we use a baby step/giant step approach, the Taylor series
can be evaluated in $\O(l^{1/2})$ operations, 
thus the evaluation requires $(n/k)^{1/2} + k$ multiplications,
and the optimal $k$ is now about $n^{1/3}$,
giving a total complexity of $\O(n^{1/3} M(n))$.
This is implemented in the function {\tt mpfr\_exp2\_aux2}.

\subsection{The error function}

The error function admits the following expansion at zero:
% \cite[formula 7.1.5]{AbSt73}:
% \[ {\rm erf} \, z = \frac{2}{\sqrt{\pi}} \sum_{k=0}^{\infty} \frac{(-1)^k}
%	{k! (2k+1)} z^{2k+1}, \]
\cite[formula 7.1.6]{AbSt73}:
\[ {\rm erf} \, z = \frac{2}{\sqrt{\pi}} e^{-z^2}
	\sum_{k=0}^{\infty} \frac{2^k}{1 \cdot 3 \cdots (2k+1)} z^{2k+1}, \]
and the following asymptotic expansion for ${\rm erfc} z = 1 - {\rm erf} z$
\cite[formula 7.1.23]{AbSt73}:
\[ \sqrt{\pi} z e^{z^2} {\rm erfc} \, z \sim 1 +
	\sum_{k=1}^{\infty} (-1)^k \frac{1 \cdot 3 \cdots (2k-1)}{(2z^2)^k}. \]
The former formula requires $m \sim n \frac{\log 2}{\log(m/(ez^2))}$ terms
% same number of terms for 7.1.5 and 7.1.6
to get $n$ correct bits, and the latter requires
$m \sim n \frac{\log 2}{\log(ez^2/m)}$ terms.
Thus, we use the expansion at $z=0$ for $n \ge e z^2$, and the asymptotic
expansion for $n < e z^2$.

\medskip

If we use the series at $z=0$, the maximum term is obtained for
$k \sim z^2$, and is of the order of $e^{z^2}$; this means
we need $z^2/(\log 2)$ additional bits. As $z^2 \le n/e$ in that case,
this is bounded by $n/(e \log 2) \sim 0.531 n$.


The series at $z=0$ is implemented as follows, 
$m$ representing the working precision,
$x, y, s, t, u$ being integer variables, and $p, r$ floating-point
variables:
\newpage
\begin{quote}
\verb|erf_0|$(z, n)$, assumes $z^2 \le n/e$ \\
$m \leftarrow n + z^2/(\log 2)$ \\
$x \leftarrow \lceil {\rm msb}(z, m) \rceil$ \\
$y \leftarrow \lceil {\rm msb}(x^2,m) \rceil$ such that $y \sim z^2 2^{e_y}$ \\
$s \leftarrow 2^n, t \leftarrow 2^n$ \\
{\bf for} $k$ {\bf from} $1$ {\bf do} \\
\q $t \leftarrow \lceil y t/k \rceil$ \\
\q $t \leftarrow \lceil t/2^{e_y} \rceil$ \\
\q $u \leftarrow \lceil t/(2k+1) \rceil$ \\
\q $s \leftarrow {\mathcal N}(s + (-1)^k u)$ \\
\q {\bf if} $u \le 1$ {\bf then} break \\
$r \leftarrow 2 \lceil z s/2^n \rceil$ \\
$p \leftarrow \pi, p \leftarrow \sqrt{p}$ \\
return $r/p$
\end{quote}
The variable $u$ contains the current term $\frac{z^{2k}}{k! (2k+1)}$,
scaled by $2^m$. Suppose $u \le 1$ for index $k_0$:
as $u \ge 2$ for index $k_0-1$, and the 
ratio between two consecutive terms decreases, then $u \le 1/2$ for index 
$k_0+1$ and the alternating series $\sum_{k_0+1}^{\infty} \frac{(-1)^k
z^{2j}}{k! (2k+1)}$
is bounded by its first term, i.e.~$2^{-m-1}$ after rescaling.

Now the relative error on $x$ is at most $2^{1-n}$,
that on $y$ is at most $2x/2^m + 1$,
that on $s$ and $t$ is zero initially.
Let $\varepsilon_k$ and $\tau_k$ the errors on $y$ and $t$
at the beginning of loop $k$,
then that for $t$ after $t \leftarrow \lceil y t/2^m \rceil$
is at most $(\varepsilon_k t + \tau_k y)/2^m + 1$.


\bibliographystyle{acm}
\bibliography{algorithms}

\end{document}