1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
\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}}
\begin{document}
\maketitle
\section{Low level functions}
\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 = \frac{\exp(-n)}{n}
\sum_{k=0}^{\infty} \frac{k!}{(-n)^k}. \]
This identity is attributed to Sweeney by Brent \cite{Brent78}.\footnote{
Brent only states that both expression for $R(n)$ are equivalent, but
they are in fact equal.}
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 = \lfloor \alpha n \rfloor$.
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}$.
\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{\exp(-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.
The remainder term is
\[ |R(n) - R'(n)| = \frac{\exp(-n)}{n}
\sum_{k=n-1}^{\infty} \frac{k!}{(-n)^k}. \]
{\bf Problem: the series for $R(n)$ diverges!}
\section{High level functions}
\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}.
\bibliographystyle{acm}
\bibliography{algo}
\end{document}
|
