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
|
\documentclass[12pt]{amsart}
\usepackage{fullpage}
\pagestyle{empty}
\title{The MPFR Library: Algorithms and Proofs}
\author{The MPFR team}
\date{\tt www.mpfr.org}
\begin{document}
\maketitle
\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}.
\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}
\bibliographystyle{acm}
\bibliography{algo}
\end{document}
|
