\documentclass[12pt]{amsart} \usepackage{fullpage,amssymb} \pagestyle{empty} \title{The MPFR Library: Algorithms and Proofs} \author{The MPFR team} \date{\tt www.mpfr.org} \def\O{{\mathcal O}} \def\n{\textnormal} \def\pinf{\bigtriangleup} \def\minf{\bigtriangledown} \def\q{\hspace*{5mm}} \def\ulp{{\rm ulp}} \def\Exp{{\rm EXP}} \def\prec{{\rm prec}} \def\sign{{\rm sign}} \def\Paragraph#1{\noindent {\sc #1}} \def\Z{{\mathcal Z}} \def\N{{\mathcal N}} \def\If{{\bf if}} \def\then{{\bf then}} \def\Else{{\bf else}} \def\elif{{\bf elif}} \def\for{{\bf for}} \def\to{{\bf to}} \def\while{{\bf while}} \def\err{{\rm err}} \newcommand{\U}[1]{\quad \mbox{[Rule~\ref{#1}]}} \newtheorem{Rule}{Rule} \begin{document} \maketitle \tableofcontents \section{Basic rules} In the following, $n$ is the precision (number of bits of the mantissa), each floating-point number is written $x = m \cdot 2^e$ with $\frac{1}{2} \le |m| < 1$ and $e := {\rm EXP}(x)$, and $\ulp(x) := 2^{{\rm EXP}(x) - n}$. \begin{Rule} \label{R1} $2^{-n} |x| < \ulp(x) \le 2^{-n+1} |x|$. \end{Rule} \begin{proof} Obvious from $x = m \cdot 2^e$ with $\frac{1}{2} \le |m| < 1$. \end{proof} \begin{Rule} \label{R2} If $a$ and $b$ have same precision $n$, and $|a| \le |b|$, then $\ulp(a) \le \ulp(b)$. \end{Rule} \begin{proof} Write $a = m_a \cdot 2^{e_a}$ and $b = m_b \cdot 2^{e_b}$. Then $|a| \le |b|$ implies $e_a \le e_b$, thus $\ulp(a) = 2^{e_a-n} \le 2^{e_b-n} = \ulp(b)$. \end{proof} \begin{Rule} \label{R3} Let $\ulp(x) := 2^{{\rm EXP}(x) - n}$, then $2^{n-1} \ulp(x) \le |x| < 2^{n} \ulp(x)$. \end{Rule} \begin{proof} Obvious from $x = m \cdot 2^e$ with $\frac{1}{2} \le |m| < 1$. \end{proof} \begin{Rule} \label{R4} $\frac{1}{2} |a| \cdot \ulp(b) < \ulp(a b) < 2 |a| \cdot \ulp(b)$. \end{Rule} \begin{proof} Write $a = m_a 2^{e_a}$, $b = m_b \cdot 2^{e_b}$, and $a b = m 2^e$ with $\frac{1}{2} \le m_a, m_b, m < 1$, then $\frac{1}{4} 2^{e_a+e_b} \le |a b| < 2^{e_a+e_b}$, thus $e = e_a + e_b$ or $e = e_a + e_b - 1$, which implies $\frac{1}{2} 2^{e_a} \ulp(b) \le \ulp(a b) \le 2^{e_a} \ulp(b)$ using $2^{e_b-n} = \ulp(b)$, and the rule follows from the fact that $|a| < 2^{e_a} \le 2|a|$ (equality on the right side can occur only if $e = e_a + e_b$ and $m_a = \frac{1}{2}$, which are incompatible). \end{proof} \begin{Rule} \label{R5} $\ulp(2^k a) = 2^k \ulp(a)$. \end{Rule} \begin{proof} Easy: if $a = m_a \cdot 2^{e_a}$, then $2^k a = m_a \cdot 2^{e_a+k}$. \end{proof} \begin{Rule} \label{R6} Let $x > 0$, $o(\cdot)$ be any rounding, and $u := o(x)$, then $\frac{1}{2} u < x < 2 u$. \end{Rule} \begin{proof} Assume $x \geq 2 u$, then $2u$ is another representable number which is closer from $x$ than $u$, which leads to a contradiction. The same argument proves $\frac{1}{2} u < x$. \end{proof} \begin{Rule} \label{R7} $\frac{1}{2} |a| \cdot \ulp(1) \leq \ulp(a) \leq |a| \cdot \ulp(1)$. \end{Rule} \begin{proof} The left inegality comes from Rule~\ref{R4} with $b=1$, and the right one from $|a| \ulp(1) \geq \frac{1}{2} 2^{e_a} 2^{1-n} =\ulp(a)$. \end{proof} \begin{Rule} \label{R8} For any $x \neq 0$ and any rounding mode $o(\cdot)$, we have $\ulp(x) \leq \ulp(o(x))$, and equality holds when rounding towards zero, towards $-\infty$ for $x>0$, or towards $+\infty$ for $x<0$. \end{Rule} \begin{proof} Without loss of generality, assume $x > 0$. Let $x = m \cdot 2^e$ with $\frac{1}{2} \leq m < 1$. As $\frac{1}{2} 2^{e_x}$ is a machine number, necessarily $o(x) \geq \frac{1}{2} 2^{e_x}$, thus by Rule~\ref{R2}, then $\ulp(o(x)) \geq 2^{e_x - n} = \ulp(x)$. If we round towards zero, then $o(x) \leq x$ and by Rule~\ref{R2} again, $\ulp(o(x)) \leq \ulp(x)$. \end{proof} \begin{Rule} \label{R9} \begin{eqnarray}\nonumber &&\n{For}\;\; error(u) \leq k_u \ulp(u),\;\; u.c_u^- \leq x \leq u.c_u^+)\\\nonumber &&\n{with}\;\; c_u^{-}=1- k_u 2^{1-p} \n{ and } c_u^{+}=1+ k_u 2^{1-p} \end{eqnarray} \begin{eqnarray}\nonumber &&\n{For}\;\; u=o(x),\;\; u.c_u^- \leq x \leq u.c_u^+\\\nonumber &&\n{if}\;\; u=\pinf(x),\n{ then } c_u^+=1\\\nonumber &&\n{if}\;\; u=\minf(x),\n{ then } c_u^-=1\\\nonumber &&\n{if}\;\; \n{for $x<0$ and } u=Z(x),\n{ then } c_u^+=1 \\\nonumber &&\n{if}\;\; \n{for $x>0$ and } u=Z(x),\n{ then } c_u^-=1 \\\nonumber &&\n{else}\;\; c_u^{-}=1-2^{1-p} \n{ and } c_u^{+}=1+2^{1-p} \end{eqnarray} \end{Rule} \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 += BITS_PER_MP_LIMB - number_of_bits(dif); 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} \subsection{The {\tt mpfr\_sub} function} The algorithm used is as follows, where $w$ denotes the number of bits per word. We assume that $a$, $b$ and $c$ denote different variables (if $a:=b$ or $a:=c$, we have first to copy $b$ or $c$), and that the rounding mode is either $\N$ (nearest), $\Z$ (towards zero), or $\infty$ (away from zero). \begin{quote} Algorithm {\tt mpfr\_sub}. \\ Input: $b$, $c$ of opposite sign with $b > c > 0$, a rounding mode $\circ \in \{ \N, \Z, \infty \}$ \\ Side effect: store in $a$ the value of $\circ(b - c)$ \\ Output: $0$ if $\circ(b - c) = b-c$, $1$ if $\circ(b - c) > b-c$, and $-1$ if $\circ(b - c) < b-c$ \\ ${\tt an} \leftarrow \lceil \frac{\prec(a)}{w} \rceil$, ${\tt bn} \leftarrow \lceil \frac{\prec(b)}{w} \rceil$, ${\tt cn} \leftarrow \lceil \frac{\prec(c)}{w} \rceil$ \\ ${\tt cancel} \leftarrow {\tt mpfr\_cmp2}(b, c)$; \quad ${\tt diff\_exp} \leftarrow \Exp(b)-\Exp(c)$ \\ ${\tt shift_b} \leftarrow (-{\tt cancel}) \bmod w$; \quad ${\tt cancel_b} \leftarrow ({\tt cancel} + {\tt shift_b})/w$ \\ \If\ ${\tt shift_b} > 0$ \then\ ${\tt b}[0 \dots \mbox{\tt bn}] \leftarrow {\tt mpn\_rshift}({\tt b}[0 \dots {\tt bn}-1], {\tt shift_b})$; ${\tt bn} \leftarrow {\tt bn} + 1$ \\ ${\tt shift_c} \leftarrow ({\tt diff\_exp}-{\tt cancel}) \bmod w$; \quad ${\tt cancel_c} \leftarrow ({\tt cancel} + {\tt shift_c}-{\tt diff\_exp})/w$ \\ \If\ ${\tt shift_c} > 0$ \then\ ${\tt c}[0 \dots \mbox{\tt cn}] \leftarrow {\tt mpn\_rshift}({\tt c}[0 \dots {\tt cn}-1], {\tt shift_c})$; ${\tt cn} \leftarrow {\tt cn} + 1$ \\ $\Exp(a) \leftarrow \Exp(b) - {\tt cancel}$; \quad $\sign(a) \leftarrow \sign(b)$ \\ $a[0 \dots {\tt an}-1] \leftarrow b[{\tt bn} - {\tt cancel_b} - {\tt an} \dots {\tt bn} - {\tt cancel_b} - 1]$ \\ $a[0 \dots {\tt an}-1] \leftarrow a[0 \dots {\tt an}-1] - c[{\tt cn} - {\tt cancel_c} - {\tt an} \dots {\tt cn} - {\tt cancel_c} - 1]$ \\ ${\tt sh} \leftarrow {\tt an} \cdot w - \prec(a)$; \quad $r \leftarrow a[0] \bmod 2^{\tt sh}$; \quad $a[0] \leftarrow a[0] - r$ \\ \If\ $\circ = \N$ and ${\tt sh} > 0$ \then \\ \q \If\ $r > 2^{{\tt sh}-1}$ \then\ $a \leftarrow a + \ulp(a)$; return $1$ \elif\ $0 < r < 2^{{\tt sh}-1}$ \then\ return $-1$ \\ \elif\ $\circ \in \{ \Z, \infty \}$ and $r > 0$ \then \\ \q \If\ $\circ = \Z$ return $-1$ \Else\ $a \leftarrow a + \ulp(a)$; return $1$ \\ ${\tt bl} \leftarrow {\tt bn} - {\tt an} - {\tt cancel_b}$ \\ ${\tt cl} \leftarrow {\tt cn} - {\tt an} - {\tt cancel_c}$ \\ \for\ $k=0$ \while\ ${\tt bl} > 0$ and ${\tt cl} > 0$ {\bf do} \\ \q ${\tt bl} \leftarrow {\tt bl} - 1$; ${\tt bp} \leftarrow b[{\tt bl}]$ \\ \q ${\tt cl} \leftarrow {\tt cl} - 1$; ${\tt cp} \leftarrow c[{\tt cl}]$ \\ \q \If\ $\circ = \N$ and $k=0$ and ${\tt sh}=0$ \then \\ \q \q \If\ ${\tt cp} \ge 2^{w-1}$ \then\ return $-1$ \\ \q \q $r \leftarrow {\tt bp} - {\tt cp}$; \quad ${\tt cp} \leftarrow {\tt cp} + 2^{w-1}$ \\ \q \If\ ${\tt bp} < {\tt cp}$ \then \\ \q \q \If\ $\circ = \Z$ \then\ $a \leftarrow a - \ulp(a)$; \quad \If\ $\circ = \infty$ \then\ return $1$ \Else\ return $-1$ \q \If\ ${\tt bp} > {\tt cp}$ \then \\ \q \q \If\ $\circ = \Z$ \then\ return $-1$ \Else\ $a \leftarrow a + \ulp(a)$; return $1$ \\ \If\ $\circ = \N$ and $r > 0$ \then \\ \q \If\ $a[0] \, {\rm div} \, 2^{\tt sh}$ is odd \then\ $a \leftarrow a + \ulp(a)$; return $1$ \Else\ return $-1$ \\ Return $0$. \end{quote} where $b[i]$ and $c[i]$ is meant as $0$ for negative $i$, and $c[i]$ is meant as $0$ for $i \ge {\tt cn}$ (${\tt cancel_b} \ge 0$, but ${\tt cancel_c}$ may be negative). \subsection{The {\tt mpfr\_div} function} The goals of the code of the {\tt mpfr\_div} function include the fact that the complexity should, as much as possible while preserving exact rounding, depend on the precision required on the result rather than on the precision given on the operands. Let $u = u_n 2^{u_e}$, $v = v_n 2^{v_e}$, where $u_n$ and $v_n$ are in $[1/2, 1[$. Let $q_p$ be the precision required on $q$. Put $b_p = \min(v_p, q_p + \varepsilon_p)$, $a_p = b_p + q_p + \varepsilon_p$, where $\varepsilon_p$ is a small value to be chosen. First, a integer division of $u_{hi} = \lfloor u_n 2^{a_p} \rfloor$ by $v_{hi} = \lfloor v_n 2^{b_p} \rfloor$ is performed. Write $u_{hi} = \tilde{q} v_{hi} + \tilde{r}$. If this division is not a full one, to obtain the real value of the quotient, if $\delta = max(u_p, v_p)$, we have to divide $u_n 2^{q_p + \varepsilon_p + \delta}$ by $v_n 2^{\delta}$. In that case, $2^{q_p + \varepsilon_p + \delta} u_n = \tilde{q}v_n 2^{\delta} + \tilde{r} 2^{\delta - q_p - \varepsilon_p} + u_{lo} - \tilde{q}v_{lo}$, with obvious notations. A positive correction on $q$ has to come from the contribution of $\tilde{r} 2^{\delta - q_p - \varepsilon_p} + u_{lo}$. The first term is at most $v_{hi} 2^{\delta - q_p - \varepsilon_p}$. As for $u_{lo}$, we have $u_{lo} < 2^{\delta-q_p-\varepsilon_p}$. Hence, the sum $u_{lo} + \tilde{r} 2^{\delta - q_p - \varepsilon_p} < 2v$, and the positive correction is at most 1. We now have to estimate $\tilde{q}v_{lo}$. It is easily seen that $\tilde{q} < 2^{q_p + \varepsilon_p + 1}$. As for $v_{lo}$, we have $v_{lo} < 2^{\delta - q_p - \varepsilon_p}$, so that $\tilde{q} v_{lo} < 2^{\delta + 1}$, to be compared with $v_n 2^{\delta}$, so that a negative correction is at most 3. As a consequence, to be able to decide rounding after the first stage, one should choose $\varepsilon_p \geq 3$ (to include the round-to-nearest case). \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: \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$. \subsection{Generic error of addition/soustraction}\label{generic:sous} We want to compute the generic error of the soustraction, this following rules can be apply on addition too. \begin{eqnarray}\nonumber \textnormal{Note:}&& error(u) \leq k_u \, \ulp(u), \;\; error(v) \leq k_v \, \ulp(v) \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&& \ulp(w)=2^{e_w-p}, \;\; \ulp(u)=2^{e_u-p},\;\; \ulp(v)=2^{e_v-p}\;\;\textnormal{with} \; p \; \textnormal{the accuracy} \\\nonumber && \ulp(u)=2^{d+e_w-p}, \;\; \ulp(u)=2^{d+e_w-p},\;\;\textnormal{with} \;\;d=e_u-e_w \;\; d^{'}=e_v-e_w \end{eqnarray} \begin{eqnarray}\nonumber error(w)& \leq &c_w \ulp(w) + k_u \ulp(u) + k_v \ulp(v) \\\nonumber &\leq&(c_w+k_u 2^d+ k_v 2^{d^{'}}) \ulp(w) \end{eqnarray} \begin{eqnarray}\nonumber &&\textnormal{If} \;\; ( u \geq 0 \;\;\textnormal{and}\;\; v \geq 0) \;\;\textnormal{or}\;\; (u \leq 0 \;\;\textnormal{and}\;\; v \leq 0) \end{eqnarray} \begin{eqnarray}\nonumber error(w)& \leq&(c_w + k_u + k_v) \, \ulp(w) \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&&\textnormal{If}\;\; w=N(u+v) \;\;\textnormal{Then}\;\; c_w =\frac{1}{2} \;\;\textnormal{else}\;\; c_w =1\\\nonumber \end{eqnarray} \subsection{Generic error of multiplication}\label{generic:mul} We want to compute the generic error of the multiplication. \begin{eqnarray}\nonumber w&=&o(u.v) \\\nonumber \textnormal{Note:}&& error(u) \leq k_u \, \ulp(u), \;\; error(v) \leq k_v \, \ulp(v) \end{eqnarray} \begin{eqnarray}\nonumber error(w)& = &|w - x.y| \\\nonumber & \leq &|w - u.b| +|u.y - x.y| \\\nonumber & \leq & c_w \ulp(w) + \frac{1}{2} [|u.v-u.y|+|u.y-x.y|+|u.v-x.v|+|x.v-x.y|]\\\nonumber & \leq & c_w \ulp(w) + \frac{u+x}{2} k_v \ulp(v) + \frac{v+y}{2} k_u \ulp(u)\\\nonumber & \leq & c_w \ulp(w) + \frac{u(1+c_u^+)}{2} k_v \ulp(v) + \frac{v(1+c_v^+)}{2} k_u \ulp(u) \U{R9}\\\nonumber & \leq & c_w \ulp(w) + (1+c_u^+) k_v \ulp(u.v) + (1+c_v^+) k_u \ulp(u.v) \U{R4}\\\nonumber & \leq & [ c_w + (1+c_u^+) k_v + (1+c_v^+) k_u ] \ulp(w)\U{R8}\\\nonumber \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&&\textnormal{If}\;\; w=N(u+v) \;\;\textnormal{Then}\;\; c_w =\frac{1}{2} \;\;\textnormal{else}\;\; c_w =1 \end{eqnarray} \subsection{Generic error of inverse}\label{generic:inv} We want to compute the generic error of the inverse. \begin{eqnarray}\nonumber w&=&o(\frac{1}{v}) \\\nonumber \textnormal{Note:}&& error(u) \leq k_u \, \ulp(u) \end{eqnarray} \begin{eqnarray}\nonumber error(w)& = &|w - \frac{1}{x}| \\\nonumber & \leq &|w - \frac{1}{u}| +|\frac{1}{u} - \frac{1}{x}| \\\nonumber & \leq & c_w \ulp(w) + \frac{1}{ux}|u-x| \\\nonumber & \leq & c_w \ulp(w) + \frac{k_u}{ux} \ulp(u) \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&& \frac{u}{c_u} \leq x\;\; \U{R6}\\\nonumber &&\n{for } u=\minf(x),\;c_u=1 \n{ else } c_u=2\\\nonumber && \n{then: } \frac{1}{x} \leq c_u \frac{1}{u} \end{eqnarray} \begin{eqnarray}\nonumber error(w)& \leq & c_w \ulp(w) + c_u\frac{k_u}{u^2} \ulp(u)\\\nonumber & \leq & c_w \ulp(w) + 2.c_u.k_u \ulp(\frac{u}{u^2}) \U{R4}\\\nonumber & \leq & [c_w + 2.c_u.k_u].\ulp(w) \U{R8} \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&&\textnormal{If}\;\; w=N(\frac{1}{u}) \;\;\textnormal{Then}\;\; c_w =\frac{1}{2} \;\;\textnormal{else}\;\; c_w =1\\\nonumber\end{eqnarray} \subsection{Generic error of division}\label{generic:div} We want to compute the generic error of the division. \begin{eqnarray}\nonumber w&=&o(\frac{u}{v}) \\\nonumber \textnormal{Note:}&& error(u) \leq k_u \, \ulp(u), \;\; error(v) \leq k_v \, \ulp(v) \end{eqnarray} \begin{eqnarray}\nonumber error(w)& = &|w - \frac{x}{y}| \\\nonumber & \leq &|w - \frac{u}{v}| +|\frac{u}{v} - \frac{x}{y}| \\\nonumber & \leq & c_w \ulp(w) + \frac{1}{vy}|uy-vx| \\\nonumber & \leq & c_w \ulp(w) + \frac{1}{vy}[|uy-xy|+|xy-vx| ]\\\nonumber & \leq & c_w \ulp(w) + \frac{1}{vy}[y k_u \ulp(u)+ x k_v \ulp(v)]\\\nonumber & \leq & c_w \ulp(w) + \frac{k_u}{v} \ulp(u)+ \frac{k_v x}{vy} \ulp(v) \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&& \frac{\ulp(u)}{v} \leq 2 \ulp(\frac{u}{v}) \;\; \U{R4}\\\nonumber && 2 \ulp(\frac{u}{v}) \leq 2 \ulp(w) \;\; \U{R8} \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&& x \leq c_u u \textnormal{ and } \frac{v}{c_v} \leq y\;\; \U{R6}\\\nonumber &&\n{with } \n{for } u=\pinf(x),\;c_u=1 \n{ else } c_u=2\\\nonumber &&\n{ and }\n{for } v=\minf(y),\;c_v=1 \n{ else } c_v=2\\\nonumber && \n{then: } \frac{x}{y} \leq c_u c_v \frac{u}{v} \end{eqnarray} \begin{eqnarray}\nonumber error(w)& \leq & c_w \ulp(w) + 2.k_u \ulp(w)+ c_u.c_v.\frac{k_v u}{vv} \ulp(v)\\\nonumber & \leq & c_w \ulp(w) + 2.k_u \ulp(w)+ 2.c_u.c_v.k_v \ulp(\frac{u.v}{v.v}) \U{R4}\\\nonumber & \leq & [c_w + 2.k_u+ 2.c_u.c_v.k_v].\ulp(w) \U{R8} \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&&\textnormal{If}\;\; w=N(\frac{u}{v}) \;\;\textnormal{Then}\;\; c_w =\frac{1}{2} \;\;\textnormal{else}\;\; c_w =1 \end{eqnarray} \subsection{Generic error of square root}\label{generic:sqrt} We want to compute the generic error of the square root. \begin{eqnarray}\nonumber v&=&o(\sqrt{u}) \\\nonumber \textnormal{Note:}&& error(u) \leq k_u \, \ulp(u)\\\nonumber \end{eqnarray} \begin{eqnarray}\nonumber error(v)& = &|v - \sqrt{x}| \\\nonumber & \leq &|v - \sqrt{u}| +|\sqrt{u} - \sqrt{x}| \\\nonumber & \leq & c_v \ulp(v) + \frac{1}{\sqrt{u} + \sqrt{x}}|u-x| \\\nonumber & \leq & c_v \ulp(v) + \frac{1}{\sqrt{u} + \sqrt{x}} k_u \ulp(u) \\\nonumber \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&& u.c_u^- \leq x \;\; \U{R9}\\\nonumber && \sqrt{u.c_u^-} \leq \sqrt{x} \\\nonumber && \sqrt{u}.(1+\sqrt{c_u^-}) \leq \sqrt{x}+\sqrt{u} \\\nonumber && \frac{1}{\sqrt{x}+\sqrt{u}} \leq \frac{1}{\sqrt{u}.(1+\sqrt{c_u^-})} \end{eqnarray} \begin{eqnarray}\nonumber error(v)& \leq & c_v \ulp(v) + \frac{1}{\sqrt{u}.(1+\sqrt{c_u^-})} k_u \ulp(u) \\\nonumber & \leq & c_v \ulp(v) + \frac{2}{1+\sqrt{c_u^-}} k_u \ulp(\sqrt{u}) \;\; \U{R4}\\\nonumber & \leq & (c_v + \frac{2.k_u}{1+\sqrt{c_u^-}}) \ulp(v) \;\; \U{R8}\\\nonumber \end{eqnarray} \subsection{Generic error of the exponential }\label{generic:exp} We want to compute the generic error of the exponential. \begin{eqnarray}\nonumber v&=&o(e^{u}) \\\nonumber \textnormal{Note:}&& error(u) \leq k_u \, \ulp(u)\\\nonumber \end{eqnarray} \begin{eqnarray}\nonumber error(v)& = &|v - e^{x}| \\\nonumber & \leq &|v - e^{u}| +|e^{u} - e^{x}| \\\nonumber & \leq & c_v \ulp(v) + e^t |u-x| \n{ with Rolle's theorem, for } t\in[x,u]\n{ or }t\in[u,x] \end{eqnarray} \begin{eqnarray}\nonumber error(v)& \leq & c_v \ulp(v) + c_u^* e^u k_u \ulp(u) \\\nonumber & \leq & c_v \ulp(v) + 2 c_u^* u k_u \ulp(e^u) \;\U{R4}\\\nonumber & \leq & (c_v + 2 c_u^* u k_u )\ulp(v) \;\U{R8}\\\nonumber & \leq & (c_v + c_u^* 2^{\Exp(u)+1} k_u )\ulp(v) \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&& u= m_u 2^{e_u} \n{ and }\ulp(u)=2^{e_u-p} \n{ with } p \n{ the accuracy} \\\nonumber \n{ Case }&x \leq u& c_u^* =1 \\\nonumber \n{ Case }&u \leq x& \\\nonumber && x \leq u + k_u \ulp(u)\\\nonumber && e^x \leq e^u e^{k_u \ulp(u)}\\\nonumber &&e^x \leq e^u e^{k_u 2^{e_u-p}}\\\nonumber &\n{then}& c_u^* = e^{k_u 2^{\Exp(u)-p}}\\\nonumber \end{eqnarray} \subsection{Generic error of the logarithm }\label{generic:log} We want to compute the generic error of the logarithm. \begin{eqnarray}\nonumber v&=&o(\log{u}) \\\nonumber \textnormal{Note:}&& error(u) \leq k_u \, \ulp(u)\\\nonumber \end{eqnarray} \begin{eqnarray}\nonumber error(v)& = &|v - \log{x}| \\\nonumber & \leq &|v - \log{u}| +|\log{u} - \log{x}| \\\nonumber & \leq & c_v \ulp(v) + \log{|\frac{x}{u}|} \\\nonumber & \leq & c_v \ulp(v) + \frac{|x-u|}{u} \\\nonumber & \leq & c_v \ulp(v) + \frac{k_u \, \ulp(u)}{u}\\\nonumber & \leq & c_v \ulp(v) + k_u \, \ulp(1)\;\; \U{R7}\\\nonumber \end{eqnarray} \begin{eqnarray}\nonumber \textnormal{Note:}&& \ulp(1)=2^{1-p}, \n{ and } \ulp(v)=2^{e_v-p} \n{ with } p \n{ the accuracy} \\\nonumber && \ulp(1)= 2^{1-e_v+e_v-p}=2^{1-e_v} \ulp(v) \end{eqnarray} \begin{eqnarray}\nonumber error(v)& \leq & c_v \ulp(v) + k_u 2^{1-e_v} \ulp(v)\\\nonumber & \leq & (c_v + k_u 2^{1-e_v} )\ulp(v)\\\nonumber \end{eqnarray} \subsection{The hyperbolic cosine function} The {\tt mpfr\_cosh} ($\cosh{x}$) function implements the hyperbolic cosine as : $$ \cosh x = \frac{1}{2} \left( e^{x} + \frac{1}{e^x} \right) $$ The algorithm used for the calculation of the hyperbolic cosine is as follows\footnote{$o()$ represent the arrondi error and $error(u)$ the error associate with the calcualtion of $u$}: \begin{eqnarray}\nonumber u&\leftarrow&o(e^x)\\\label{coshalgo1} v&\leftarrow&o({u}^{-1})\\\label{coshalgo2} w&\leftarrow&o(u+v)\\\label{coshalgo3} s&\leftarrow&\frac{1}{2} w\\\label{coshalgo4} \end{eqnarray} Now, we have to bound the rounding error for each step of this algorithm. First, let consider the parity of hyperbolic cosine ($\cosh(-x)=\cosh(x)$) : the problem is reduced to calculate $\cosh x$ with $x \geq 0$. We can deduce $e^x \geq 1$ and $0 \leq e^{-x} \leq 1$. \begin{center} \begin{tabular}{l l l} \begin{minipage}{2.5cm} ${\textnormal{error}}(u)$ $u \leftarrow o(e^x)$\\ $-\infty \;\; (\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber |u-e^x| &\leq& \ulp(u)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} {\hspace{7cm}} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(v)$ $v \leftarrow o({u}^{-1}) $\\ $+\infty \;\; (\bullet\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|v-e^{-x}| \\\nonumber & \leq& |v - u^{-1}| + |u^{-1} - e^{-x}|\\\nonumber & \leq& \ulp(v) + \frac{1}{u \cdot e^x} |u-e^{x}|\\\nonumber & \leq& \ulp(v) + \frac{1}{u^2} \ulp(u)\;\;(\star)\\\nonumber & \leq& \ulp(v) + 2 \ulp(\frac{1}{u})\;\;(\star\star)\\\nonumber & \leq& 3 \, \ulp(v)\;\;(\star\star\star) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} $(\star)$ With $\frac{1}{e^x} \leq \frac{1}{u}$,\\ for that we must have $u \leq e^x$,\\ it is possible with a rounding of\\ $u$ to $-\infty \;\; (\bullet)$ $(\star\star)$ From inequation \U{R4}, \[ a \cdot \ulp(b) \leq 2 \cdot \ulp(a \cdot b)\] if $a =\frac{1}{u^2},\;b = u$ then \[ \frac{1}{u^2} \ulp(u) \leq 2 \ulp(\frac{1}{u})\] $(\star\star\star)$ If $\ulp(\frac{1}{u}) \leq ulp(v)$,\\ it is possible with a rounding of \\ $v$ to $+\infty \;\; (\bullet)$\\ \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(w)$ $w \leftarrow o(u+v) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|w-(e^{x}+e^{-x})| \\\nonumber & \leq& |w - (u+v)|+|u - e^x|+|v - e^{-x}|\\\nonumber & \leq& \ulp(w) + \ulp(u) + 3\ulp(v)\\\nonumber & \leq& \ulp(w) + 4\ulp(u)\;\;(\star)\\\nonumber & \leq& 5\ulp(w)\;\;(\star\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} $(\star)$ With $v \leq 1\leq u$ then $\ulp(v) \leq \ulp(u)$ $(\star\star)$ With $u \leq w$ then $\ulp(u) \leq \ulp(w)$ \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(s)$ $s \leftarrow o(\frac{w}{2}) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{center} \begin{eqnarray}\nonumber {\textnormal{error}}(s) & = & {\textnormal{error}}(w)\\\nonumber & \leq & 5\ulp(s) \end{eqnarray} \end{center} \end{minipage} & \begin{minipage}{6cm} \end{minipage} \end{tabular} \end{center} That shows the rounding error on the calculation of $\cosh x$ can be bound by $5 \;\; \ulp$ on the result. So, to calculate the size of intermediary variables, we have to add, at least, $\lceil \log_2 5 \rceil=3$ bits the wanted precision. \subsection{The inverse hyperbolic cosine function} The {\tt mpfr\_acosh} ($\n{acosh}{x}$) function implements the inverse hyperbolic cosine as : $$ \n{acosh} = \log \left( \sqrt{x+1} \sqrt{x-1} + x \right) $$ The algorithm used for the calculation of the inverse hyperbolic cosine is as follows \begin{eqnarray}\nonumber q&\leftarrow&o(x+1)\\\nonumber r&\leftarrow&o(x-1)\\\nonumber s&\leftarrow&o(\sqrt{q})\\\nonumber t&\leftarrow&o(\sqrt{r})\\\nonumber u&\leftarrow&o(s \times t)\\\nonumber v&\leftarrow&o(u+x)\\\nonumber w&\leftarrow&o(\log v) \end{eqnarray} Now, we have to bound the rounding error for each step of this algorithm. First, let consider the function field : {\tt mpfr\_acosh} is define for $x \geq 1$. \begin{center} \begin{tabular}{l l l} \begin{minipage}{2.5cm} ${\textnormal{error}}(q)$ $q \leftarrow \minf(x+1) $ $(\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|q-(x+1)| \\\nonumber & \leq& 2 \ulp(q)\;\;(\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:sous} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(r)$ $r \leftarrow \minf(x-1) $ $(\bullet\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|r-(x-1)| \\\nonumber & \leq& (1+2^{\Exp(x)-\Exp(r)}) \ulp(r)\;\;(\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:sous} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(s)$ $s \leftarrow \pinf(\sqrt{q}) $ ($\bullet\bullet\bullet$) \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|s-\sqrt{x+1}| \\\nonumber & \leq& \ulp(s) + |\sqrt{q}-\sqrt{x+1}|\\\nonumber & \leq& \ulp(s) + \frac{1}{\sqrt{q} + \sqrt{x+1}}|q-(x+1)|\\\nonumber & \leq& \ulp(s) + \frac{1}{\sqrt{q} + \sqrt{x+1}}.2.\ulp(q) \\\nonumber & \leq& \ulp(s) + \frac{1}{2\sqrt{q}}.2.\ulp(q) \;\;(\star) \\\nonumber & \leq& \ulp(s) + 2.\ulp(\sqrt{q}) \;\;(\star\star)\\\nonumber & \leq& (1 + 2) \ulp(s) \;\;(\star\star\star) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) If $q \leftarrow \minf(x+1) \;(\bullet)$ Then $q \leq x+1$ or $\frac{1}{\sqrt{x+1}+\sqrt{v}} \leq \frac{1}{2.\sqrt{q}}$ ($\star\star$) $\U{R4}$ ($\star\star\star$) $\U{R8}$ \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(t)$ $t \leftarrow \pinf(\sqrt{r}) $ ($\bullet\bullet\bullet\bullet$) \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|t-\sqrt{x-1}| \\\nonumber & \leq& \ulp(t) + |\sqrt{r}-\sqrt{x-1}|\\\nonumber & \leq& \ulp(t) + \frac{1}{\sqrt{r} + \sqrt{x-1}}|r-(x-1)|\\\nonumber & \leq& \ulp(t) + \frac{1}{\sqrt{r} + \sqrt{x+1}} \\\nonumber & \cdots & (1+2^{\Exp(x)-\Exp(r)}) \ulp(r) \\\nonumber & \leq& \ulp(t) + \frac{1}{2\sqrt{r}} \cdots \;\;(\star) \\\nonumber &\cdots& (1+2^{\Exp(x)-\Exp(r)}) \ulp(r) \;\;(\star) \\\nonumber & \leq& \ulp(t) + (1+2^{\Exp(x)-\Exp(r)})\ulp(\sqrt{r}) \;\;(\star\star)\\\nonumber & \leq& (2+2^{\Exp(x)-\Exp(r)}) \ulp(t) \;\;(\star\star\star) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) If $r \leftarrow \minf(x-1) \;(\bullet\bullet)$ Then $q \leq x+1$ or $\frac{1}{\sqrt{x-1}+\sqrt{r}} \leq \frac{1}{2.\sqrt{r}}$ ($\star\star$) $\U{R4}$ ($\star\star\star$) $\U{R8}$ \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(u)$ $u \leftarrow o(t \times s) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|u-(\sqrt{x+1} \sqrt{x-1})| \\\nonumber & \leq& (1+2 \times 3 +2 \times (2+2^{\Exp{x}-\Exp{r}}))\\\nonumber & \cdots& \ulp(w) \;(\star)\\\nonumber & \leq& (13+2^{\Exp{x}-\Exp{r}+1}) \ulp(u) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:mul} with $(\bullet\bullet\bullet)$ and $(\bullet\bullet\bullet\bullet)$ \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(v)$ $v \leftarrow o(u+x) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|v-(\sqrt{x+1} \sqrt{x-1} +x)| \\\nonumber & \leq& (15+2^{\Exp(x)-\Exp(r)+1}) \ulp(v) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:sous} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(w)$ $w \leftarrow o(\log{v}) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|w-\log(\sqrt{x+1} \sqrt{x-1} +x)| \\\nonumber & \leq& (1+(15+2^{\Exp(x)-\Exp(r)+1}).2^{1-\Exp(w)} \ulp(w) \\\nonumber & \leq& (1+15.2^{1-\Exp(w)}+2^{\Exp(x)-\Exp(r)-\Exp(w)+2}) \ulp(w) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:log} \end{minipage} \end{tabular} \end{center} That shows the rounding error on the calculation of $\n{acosh} x$ can be bound by $ (1+15.2^{1-\Exp(w)}+2^{\Exp(x)-\Exp(r)-\Exp(w)+2})\;\; \ulp$ on the result. So, to calculate the size of intermediary variables, we have to add, at least, $\lceil \log_2 (1+15.2^{1-\Exp(w)}+2^{\Exp(x)-\Exp(r)-\Exp(w)+2}) \rceil$ bits the wanted precision. \subsection{The hyperbolic sine function} The {\tt mpfr\_sinh} ($\sinh{x}$) function implements the hyperbolic sine as : $$ \sinh x = \frac{1}{2} \left( e^{x} - \frac{1}{e^x} \right) $$ The algorithm used for the calculation of the hyperbolic sine is as follows\footnote{$o()$ represent the arrondi error and $error(u)$ the error associate with the calcualtion of $u$}: \begin{eqnarray}\nonumber u&\leftarrow&o(e^x)\\\nonumber v&\leftarrow&o({u}^{-1})\\\nonumber w&\leftarrow&o(u-v)\\\nonumber s&\leftarrow&\frac{1}{2} w \end{eqnarray} Now, we have to bound the rounding error for each step of this algorithm. First, let consider the parity of hyperbolic sine ($\sinh(-x)=-\sinh(x)$) : the problem is reduced to calculate $\sinh x$ with $x \geq 0$. We can deduce $e^x \geq 1$ and $0 \leq e^{-x} \leq 1$. \begin{center} \begin{tabular}{l l l} \begin{minipage}{2.5cm} ${\textnormal{error}}(u)$ $u \leftarrow \minf(e^x)$\\ $(\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber |u-e^x| &\leq& \ulp(u)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} {\hspace{7cm}} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(v)$ $v \leftarrow \pinf({u}^{-1}) $\\ $(\bullet\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|v-e^{-x}| \\\nonumber & \leq& |v - u^{-1}| + |u^{-1} - e^{-x}|\\\nonumber & \leq& \ulp(v) + \frac{1}{u \cdot e^x} |u-e^{x}|\\\nonumber & \leq& \ulp(v) + \frac{1}{u^2} \ulp(u)\;\;(\star)\\\nonumber & \leq& \ulp(v) + 2 \ulp(\frac{1}{u})\;\;(\star\star)\\\nonumber & \leq& 3 \, \ulp(v)\;\;(\star\star\star) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} $(\star)$ With $\frac{1}{u} \leq \frac{1}{e^x}$,\\ for that we must have $e^x \leq u$,\\ it is possible with $u=\minf(e^x)$ $(\bullet)$ $(\star\star)$ From inequation \U{R4}, \[ a \cdot \ulp(b) \leq 2 \cdot \ulp(a \cdot b)\] if $a =\frac{1}{u^2},\;b = u$ then \[ \frac{1}{u^2} \ulp(u) \leq 2 \ulp(\frac{1}{u})\] $(\star\star\star)$ If $\ulp(\frac{1}{u}) \leq \ulp(v)$,\\ it is possible with $v=\pinf(u^{-1})$ $(\bullet\bullet)$ \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(w)$ $w \leftarrow o(u+v) $ \end{minipage} & \begin{minipage}{7.8cm} \begin{eqnarray}\nonumber &&|w-(e^{x}-e^{-x})| \\\nonumber & \leq& |w - (u-v)|+|u - e^x|+|-v + e^{-x}|\\\nonumber & \leq& \ulp(w) + \ulp(u) + 3\ulp(v)\\\nonumber & \leq& \ulp(w) + 4\ulp(u)\;\;(\star)\\\nonumber & \leq& (1+ 4 \cdot 2^{\Exp(u)-\Exp(w)}) \ulp(w)\;\;(\star\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} $(\star)$ With $v \leq 1\leq u$ then $\ulp(v) \leq \ulp(u)$ $(\star\star)$ see subsection \ref{generic:sous} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(s)$ $s \leftarrow o(\frac{w}{2}) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{center} \begin{eqnarray}\nonumber {\textnormal{error}}(s) & = & {\textnormal{error}}(w)\\\nonumber & \leq & (1+ 4 \cdot 2^{\Exp(u)-\Exp(w)}) \ulp(w) \end{eqnarray} \end{center} \end{minipage} & \begin{minipage}{6cm} \end{minipage} \end{tabular} \end{center} That show the rounding error on the calculation of $\sinh x$ can be bound by $(1+ 4 \cdot 2^{\Exp(u)-\Exp(w)}) \ulp(w)$, then the number of bits need to add to the want accuracy to define intermediary variable is : \[ N_t=\lceil \log_2(1+ 4 \cdot 2^{\Exp(u)-\Exp(w)}) \rceil \] \subsection{The inverse hyperbolic sine function} The {\tt mpfr\_asinh} ($\n{acosh}{x}$) function implements the inverse hyperbolic sine as : $$ \n{asinh} = \log \left( \sqrt{x^2+1} + x \right) $$ The algorithm used for the calculation of the inverse hyperbolic sine is as follows \begin{eqnarray}\nonumber s&\leftarrow&o(x^2)\\\nonumber t&\leftarrow&o(s+1)\\\nonumber u&\leftarrow&o(\sqrt{t})\\\nonumber v&\leftarrow&o(u+x)\\\nonumber w&\leftarrow&o(\log v) \end{eqnarray} Now, we have to bound the rounding error for each step of this algorithm. First, let consider the parity of hyperbolic arcsine ($\n{asinh}(-x)=-\n{asinh}(x)$) : the problem is reduced to calculate $\n{asinh} x$ with $x \geq 0$. \begin{center} \begin{tabular}{l l l} \begin{minipage}{2.5cm} ${\textnormal{error}}(s)$ $s \leftarrow o(x^2) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|s-x^2| \\\nonumber & \leq& \ulp(s)\;\;(\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(t)$ $t \leftarrow \minf(s+1) $ $(\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|t-(x^2+1)| \\\nonumber & \leq& 2 \ulp(t) \;\;(\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:sous} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(u)$ $u \leftarrow o(\sqrt{t}) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|u-\sqrt{x^2+1}| \\\nonumber & \leq& 3 \ulp(u) \;(\star) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:sqrt} with ($\bullet$) \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(v)$ $v \leftarrow o(u+x) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|v-(\sqrt{x^2+1}+x)| \\\nonumber & \leq& 5 \ulp(v) \;(\star) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:sous} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(w)$ $w \leftarrow o(\log v) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|w-\log(\sqrt{x^2+1}+x)| \\\nonumber & \leq& (1+5.2^{1-\Exp(w)}) \ulp(w) \;\star \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:log} \end{minipage} \end{tabular} \end{center} That shows the rounding error on the calculation of $\n{asinh} x$ can be bound by $ (1+5.2^{1-\Exp(w)})\;\; \ulp$ on the result. So, to calculate the size of intermediary variables, we have to add, at least, $\lceil \log_2 (1+5.2^{1-\Exp(w)}) \rceil$ bits the wanted precision. \subsection{The hyperbolic tangent function} The {\tt mpfr\_tanh} ($\tanh{x}$) function implements the hyperbolic tangent as : $$ \tanh x = \frac{ e^{2 \cdot x} -1 }{ e^{2 \cdot x} +1} $$ The algorithm used for the calculation of the hyperbolic tangent is as follows\footnote{$o()$ represent the arrondi error and $error(u)$ the error associate with the calcualtion of $u$}: \begin{eqnarray}\nonumber u&\leftarrow&o(2 \cdot x)\\\nonumber v&\leftarrow&o(e^u)\\\nonumber w&\leftarrow&o(v+1)\\\nonumber r&\leftarrow&o(v-1)\\\nonumber s&\leftarrow&o(\frac{r}{w}) \end{eqnarray} Now, we have to bound the rounding error for each step of this algorithm. First, let consider the parity of hyperbolic tangent ($\tanh(-x)=-\tanh(x)$) : the problem is reduced to calculate $\tanh x$ with $x \geq 0$. We can deduce $e^x \geq 1$. \begin{center} \begin{tabular}{l l l} \begin{minipage}{2.5cm} ${\textnormal{error}}(u)$ $u \leftarrow o(2 \cdot x)$ \end{minipage} & \begin{minipage}{7.5cm} exact \end{minipage} & \begin{minipage}{6cm} {\hspace{7cm}} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(v)$ $v \leftarrow o(e^{u}) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|v-e^{2 \cdot x}| \\\nonumber & \leq& |v - e^{u}| + |e^{u} - e^{2 \cdot x}|\\\nonumber & \leq& \ulp(v) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(w)$ $w \leftarrow \minf(v+1) $ $(\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|w-(e^{2 \cdot x}+1)| \\\nonumber &\leq& |w - (v+1)|+|(v+1) - (e^{2 \cdot x}+1)|\\\nonumber &\leq& \ulp(w) + \ulp(v)\\\nonumber &\leq& 2 \cdot \ulp(w)\;(\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} $(\star)$ With $v \leq w$ then $\ulp(v) \leq \ulp(w)$ \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(r)$ $r \leftarrow \pinf(v-1) $ $(\bullet\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|r-(e^{2 \cdot x}-1)| \\\nonumber &\leq& |r - (v-1)|+|(v-1) - (e^{2 \cdot x}-1)|\\\nonumber &\leq& \ulp(r) + \ulp(v)\\\nonumber &\leq& (1+2^{\Exp(v)-\Exp(r)}) \ulp(r)\;(\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} $(\star)$ see subsection \ref{generic:sous} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(s)$ $s \leftarrow o(\frac{r}{w}) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|s-\frac{e^{2x}-1}{e^{2x}+1}| \\\nonumber &\leq& (1+2 \times 2+ \hdots\\\nonumber &\hdots& 2(1+2^{\Exp(v)-\Exp(r)})) \ulp(s) \;\;(\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} $(\star)$ see subsection \ref{generic:div} with $(\bullet)$ and $(\bullet\bullet)$ \end{minipage} \end{tabular} \end{center} That show the rounding error on the calculation of $\tanh x$ can be bound by $(1+2 \times 2+2(1+2^{\Exp(v)-\Exp(r)})) \ulp(s)$, then the number of bits need to add to the want accuracy to define intermediary variable is : \[ N_t=\lceil \log_2(7+2^{\Exp(v)-\Exp(r)+1}) \rceil \] \subsection{The inverse hyperbolic tangent function} The {\tt mpfr\_atanh} ($\n{acosh}{x}$) function implements the inverse hyperbolic tangent as : $$ \n{atanh} = \frac{1}{2} \log \frac{1+x}{1-x} $$ The algorithm used for the calculation of the inverse hyperbolic tangent is as follows \begin{eqnarray}\nonumber s&\leftarrow&o(1+x)\\\nonumber t&\leftarrow&o(1-x)\\\nonumber u&\leftarrow&o(\frac{s}{t})\\\nonumber v&\leftarrow&o(\log u)\\\nonumber w&\leftarrow&o(\frac{1}{2} v) \end{eqnarray} Now, we have to bound the rounding error for each step of this algorithm. First, let consider the parity of hyperbolic arctan ($\n{atanh}(-x)=-\n{atanh}(x)$) : the problem is reduced to calculate $\n{atanh} x$ with $x \geq 0$. \begin{center} \begin{tabular}{l l l} \begin{minipage}{2.5cm} ${\textnormal{error}}(s)$ $s \leftarrow \pinf(1+x) $ $(\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|s-(1+x)| \\\nonumber & \leq& 2 \ulp(s)\;\;(\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} see subsection \ref{generic:sous} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(t)$ $t \leftarrow \minf(1-x) $ $(\bullet\bullet)$ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|t-(1-x)| \\\nonumber & \leq& (1+2^{\Exp(x)-\Exp(t)}) \ulp(t) \;\;(\star)\\\nonumber \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:sous} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(u)$ $u \leftarrow o(\frac{s}{t}) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|u-\frac{1+x}{1-x}| \\\nonumber & \leq& (1+ 2 \times 2 + \\\nonumber & \cdots& 2 \times (1+2^{\Exp(x)-\Exp(t)}))\ulp{u} \;(\star)\\\nonumber & \leq& (7+2^{\Exp(x)-\Exp(t)+1})\ulp(u) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:inv} with ($\bullet$) and ($\bullet\bullet$) \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(v)$ $v \leftarrow o(\log(u)) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|v-(\log{\frac{1+x}{1-x}})| \\\nonumber & \leq& (1+(7+2^{\Exp(x)-\Exp(t)+1}) \\\nonumber & \cdots& \times 2^{1-\Exp(v)}) \ulp(v)\;(\star)\\\nonumber & \leq& (1+7 \times 2^{1-\Exp(v)} +\\\nonumber & \cdots& 2^{\Exp(x)-\Exp(t)-\Exp(v)+2}) \ulp(v) \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) see subsection \ref{generic:log} \end{minipage}\\%\hline \begin{minipage}{2.5cm} ${\textnormal{error}}(w)$ $w \leftarrow o(\frac{1}{2} v) $ \end{minipage} & \begin{minipage}{7.5cm} \begin{eqnarray}\nonumber &&|w-\frac{1}{2}\log{\frac{1+x}{1-x}}| \\\nonumber & \leq& (1+7 \times 2^{1-\Exp(v)} + \\\nonumber & \cdots& 2^{\Exp(x)-\Exp(t)-\Exp(v)+2}) \ulp(w) \;\star \end{eqnarray} \end{minipage} & \begin{minipage}{6cm} ($\star$) exact \end{minipage} \end{tabular} \end{center} That shows the rounding error on the calculation of $\n{atanh} x$ can be bound by $ (1+7 \times 2^{1-\Exp(v)} + 2^{\Exp(x)-\Exp(t)-\Exp(v)+2})\;\; \ulp$ on the result. So, to calculate the size of intermediary variables, we have to add, at least, $\lceil \log_2 (1+7 \times 2^{1-\Exp(v)} + 2^{\Exp(x)-\Exp(t)-\Exp(v)+2}) \rceil$ bits the wanted precision. \subsection{The arc-sine function} \begin{enumerate} \item We use the formula $arcsin\,x=\arctan\,\frac{x}{\sqrt{1-x^2}}$ \item We will have the when $x$ is near $1$ we will experience uncertainty problems: \item If $x=a(1+\delta)$ with $\delta$ being the relative error then we will have \begin{equation*} 1-x=1-a-a\delta=(1-a)[1-\frac{a}{1-a}\delta] \end{equation*} Ans so when using the arctangent programs we need to take into account that decrease in precision. \item We will have \end{enumerate} \subsection{The arc-cosine function} % from Mathieu Dutour \begin{enumerate} \item Obviously, we used the formula \begin{equation*} \arccos\,x=\frac{\pi}{2}-\arcsin\,x \end{equation*} \item The problem of $\arccos$ is that it is $0$ at $1$, so, we have a cancellation problem to treat at $1$. \item (Suppose $x\geq 0$, this is where the problem happens) The derivative of $\arccos$ is $\frac{-1}{\sqrt{1-x^2}}$ and we will have \begin{equation*} \frac{1}{2\sqrt{1-x}} \leq |\frac{-1}{\sqrt{1-x^2}}|=\frac{1}{\sqrt{(1-x)(1+x)}} \leq \frac{1}{\sqrt{1-x}} \end{equation*} So, integrating the above inequality on $[x,1]$ we get \begin{equation*} \sqrt{1-x}\leq \arccos\,x\leq 2\sqrt{1-x} \end{equation*} \item The important part is the lower bound that we get which tell us a upper bound on the cancellation that will occur:\\ The terms that are canceled are $\pi/2$ and $\arcsin\,x$, their order is $2$. The number of canceled terms is so \begin{verbatim} 1-1/2*MPFR_EXP(1-x) \end{verbatim} \end{enumerate} \subsection{The arc-tangent function} % from Mathieu Dutour \subsubsection{Binary splitting} \noindent The Taylor serie for $\arctan$ is suitable for analysis using Binary splitting. \par This method is detailed for example in ``Pi and The AGM'' p 334. It is efficient for rational numbers and is non efficient for non rational numbers. \par The efficiency of this method is then quite limited. One can then wonder how to use it for non rational numbers. \par Using the formulas \begin{equation*} \arctan\,(-x)=-\arctan\,x\;\;\mbox{and}\;\;\arctan\,x+\arctan\,\frac{1}{x}=\frac{\pi}{2}{\rm sign}(x) \end{equation*} we can restrict ourself to $0\leq x\leq 1$. \par Writing \begin{equation*} x=\sum_{i=1}^{\infty} \frac{u_i}{2^i}\;\;\mbox{with}\;\;u_i\in\{0,1\} \end{equation*} or \begin{equation*} x=\sum_{i=1}^{\infty} \frac{u_i}{2^{2^i}}\;\;\mbox{with}\;\;u_i\in\{0,1,\dots,2^{2^{i-1}}_1\}\mbox{~if~}i>1\mbox{~and~}u_1\in \{0,1\} \end{equation*} we can compute $\cos$, $\sin$ or $\exp$ using the formulas \begin{equation*} \begin{array}{c} \cos\,(a+b)=\cos\,a\cos\,b-\sin\,a\sin\,b\\ \sin\,(a+b)=\sin\,a\cos\,b+\cos\,a\sin\,b\\ \exp(a+b)=(\exp\,a)(\exp\,b) \end{array} \end{equation*} Unfortunately for $\arctan$ there is no similar formulas. The only formula known is \begin{equation*} \arctan\,x+\arctan\,y=\arctan\,\frac{x+y}{1-xy}+k\pi\;\;\mbox{with}\;\;k\in\Z \end{equation*} we will use \begin{equation*} \arctan\,x=\arctan\,y+\arctan\,\frac{x-y}{1+xy} \end{equation*} with $x,y>0$ and $y0 $, with $\frac{1}{8}\leq N <1$, $\frac{2x_k+N^{\frac{1}{3}}}{3x^2_k} >0 $ and then $\epsilon_{k+1} >0$ in othe hand $ N^{\frac{1}{3}}\leq x_{k+1}$, we can obtain this inegality for $1 \leq k$: \begin{equation} N^{\frac{1}{3}}\leq x_{k} \end{equation} We can now bound $\epsilon_{k+1}$ : \begin{eqnarray}\nonumber \epsilon_{k+1} = \epsilon^2_{k} \frac{2x_k+N^{\frac{1}{3}}}{3x^2_k} \leq \epsilon^2_{k}\frac{2x_k+x_k}{3x^2_k}\leq 2 \epsilon^2_{k}\\\nonumber \end{eqnarray} Because $\frac{1}{2}\leq n + \log_2 n_k + \log_2 n \] \begin{table} \begin{center} \begin{tabular}{|c|c|c|} \hline $w$ & $\err(w)/\ulp(w) \le c_w + \ldots$ &special case\\ \hline\hline $o(u+v)$ & $k_u 2^{e_u-e_w} + k_v 2^{e_v-e_w}$ & $k_u + k_v$ if $u v \ge 0$\\ $o(u \cdot v)$ & $(1+c^{+}_u)k_u + (1+c^{+}_v)k_v$ & $2k_u + 2k_v$ if $u \ge x$, $v \ge y$\\ $o(1/v)$ & $4 k_v$ & $2 k_v$ if $v \le y$ \\ $o(u/v)$ & $4 k_u + 4 k_v$ & $2 k_u + 2 k_v$ if $v \le y$ \\ $o(\sqrt{u})$ & $2 k_u/(1+\sqrt{c^{-}_u})$ & $k_u$ if $u \le x$ \\ $o(e^u)$ & $e^{k_u 2^{e_u-p}} 2^{e_u+1} k_u$ & $2^{e_u+1} k_u$ if $u \ge x$ \\ $o(\log u)$ & $k_u 2^{1-e_w}$ & \\ \hline \end{tabular} \end{center} \caption{Generic error for several operations, assuming all variables have a mantissa of $p$ bits, and no overflow/underflow occurs. The inputs $u$ and $v$ are approximations of $x$ and $y$ with $|u-x| \le k_u \ulp(u)$ and $|v-y| \le k_v \ulp(v)$. The additional rounding error $c_w$ is $1/2$ for rounding to nearest, and $1$ otherwise. The value $c^{\pm}_u$ equals $1 \pm k_u 2^{1-p}$. } \end{table} Remark 1: in the addition case, when $u v > 0$, necessarily one of $2^{e_u-e_w}$ and $2^{e_v-e_w}$ is less than $1/2$, thus $\err(w)/\ulp(w) \le c_w + {\rm max}(k_u + k_v/2, k_u/2 + k_v) \le c_w + \frac{3}{2} {\rm max}(k_u, k_v)$. \nocite{BoBo98} \bibliographystyle{acm} \bibliography{algorithms} \end{document}