xref: /dports/math/reduce/Reduce-svn5758-src/libraries/crlibm/docs/latex/log.tex

There are two versions of the logarithm.
\begin{itemize}
\item The first relies on 80-bit double-extended arithmetic, and is
  well suited to IA32 and IA64 architectures which have hardware
  support for such arithmetic. It computes the quick step in
  double-extended arithmetic, and the accurate step in
  double-double-extended arithmetic.
\item The second relies only on double-precision arithmetic, and is
  portable. It uses double-double for the quick step, and
  triple-double for the accurate step.
\end{itemize}

Both implementations use the same algorithm, which is detailed in
\ref{sec:logoutline}. Sections \ref{sec:logdeproof} and
\ref{sec:logtdproof} detail the proof of both implementations, and
\ref{sec:logperf} give some performance results.


\section{General outline of the algorithm\label{sec:logoutline}}
The algorithm used is mainly due to
Wong and Goto\cite{WG94} and has been discussed further in \cite{Muller97}. In the case we are given here,
both quick and accurate phase use principally the same algorithm however optimized for different accuracies.

The function's argument $x \in \F$ is first checked for special cases, such as $x\leq0$, $+\infty$, $\nan$ etc.
These checks are mainly implemented using integer arithmetics and will be further explained in section
\ref{subsec:reduction}. Then, the argument is reduced using integer arithmetics as follows:
$$x = 2^{E^\prime} \cdot m$$
where $E^\prime$ is the exponent of $x$ and $m$ a double corresponding
to the mantissa of $x$. This decomposition is done such that in any
case, i.e. even if $x$ is subnormal, $1 \leq m < 2$. In the subnormal
case, the exponent of $x$ is adjusted accordingly.  This first
argument reduction corresponds to the equality
$$\log\left( x \right) = E^\prime \cdot \log\left(2\right) + \log\left(m \right)$$
Using this term directly would lead to catastrophic cancellation in the case where $E^\prime = -1$ and
$m \approx 2$. To overcome this difficulty, a second adjustment is done as follows:
\begin{center}
  \begin{tabular}{cc}
    \begin{minipage}{50mm}
      $$E = \left \lbrace \begin{array}{lcl} E^\prime & \mbox{ if } & m \leq \sqrt{2} \\
          E^\prime +1 & \mbox{ if } & m > \sqrt{2} \end{array} \right.$$
    \end{minipage}
    &
    \begin{minipage}{50mm}
      $$y = \left \lbrace \begin{array}{lcl} m & \mbox{ if } & m \leq \sqrt{2} \\
          \frac{m}{2} & \mbox{ if } & m > \sqrt{2} \end{array} \right.$$
    \end{minipage}
  \end{tabular}
\end{center}

The decision whether $m \leq \sqrt{2}$ or not is performed using integer arithmetics on
the high order bits of the mantissa $m$. The test is therefore not completely exact which is no
disadvantage since, in any case, the bound $\sqrt{2}$ is somewhat arbitrary.\par
All the previous reduction steps can be implemented exactly as they consist mainly in decompositions
of a floating point number, multiplications by powers of $2$ and integer additions on the corresponding exponent value.
All this leads to the following equation
$$\log\left( x \right) = E \cdot \log\left( 2 \right) + \log\left( y \right)$$
where
$$-\frac{1}{2} \cdot \log\left( 2 \right) \leq \log\left( y \right) \leq \frac{1}{2} \cdot \log\left( 2 \right)$$
The magnitude of $y$ is thus still too great for allowing for a direct polynomial approximation of $\log\left(y\right)$.
Therefore, a second argument reduction step is performed using a table of $128$ entries as follows:
using the high order bits of $y$ as an index $i$, a tabulated value $r_i$ is looked up which approximated very well
$\frac{1}{y}$. Setting $z = y \cdot r_i - 1$, one obtains
$$\log\left( y \right) = \log\left( 1 + z \right) - \log\left( r_i \right)$$
Since $y = \frac{1}{r_i} + \delta$ the magnitude of $z$ is finally
small enough (typically $\left \vert z \right \vert < 2^{-8}$) for
approximating $\log\left(1+z\right)$ by a Remez polynomial
$p\left(z\right)$. The values for $\log\left(r_i\right)$ are of course also tabulated.

It is important to notice that the reduction step
$$z = y \cdot r_i - 1$$
can be implemented exactly which eases the correctness proof of the algorithm. This property will be proven in
section \ref{subsec:reduction}. The reduced argument $z$ will
be represented as a double-double number $z_\hi + z_\lo$ that will be fed into the polynomial approximation
algorithms of both quick and accurate phase. Each of these phases will take into account the lower significant value
$z_\lo$ for more or less
higher monomial degrees.

Both phases will finally reconstruct the function's value as follows:
$$\log\left( x \right) \approx E \cdot \log\left( 2 \right) + p\left( z \right) - \log\left( r_i \right)$$
using a double (respectively a triple for the accurate phase) double value for each
$\log\left( 2 \right)$ and $-\log\left( r_i \right)$. The computations necessary for performing this reconstruction
are carried out in double-double arithmetics for the quick phase and triple-double for the accurate phase.

The quick phase uses a modified Remez polynomial of degree $7$ of the form
$$p\left( z \right) = z - \frac{1}{2} \cdot z^2 + z^3 \cdot
\left( c_3 + z \cdot \left( c_4 + z \cdot \left( c_5 + z \cdot \left( c_6 + z \cdot c_7 \right) \right) \right) \right)$$
with $c_i \in \F$.
This polynomial is evaluated as indicated by the parenthesis in the following term:
$$p\left( z_\hi + z_\lo \right) \approx \left( \left(z_\hi + z_\lo \right) - \frac{1}{2} \cdot  z_\hi^2\right) +
\left( \left( - z_\hi \cdot z_\lo \right) +
\left(z_\hi^2 \cdot z_\hi \right) \cdot
\left( c_3 + z_\hi \cdot \left( c_4 + z_\hi \cdot \left( c_5 + z_\hi \cdot \left( c_6 + z_\hi \cdot c_7 \right) \right) \right) \right) \right)$$
The mathematical relative approximation error of the polynomial $p\left( z \right)$ defined as
$$\epsilon_{\mbox{\tiny meth}} = \frac{p\left( z \right) - \log\left( 1 + z \right)}{\log\left(1 + z \right)}$$ is bounded by
$$\left \vert \epsilon_{\mbox{\tiny meth}} \right \vert \leq 2^{-62.99}$$
This methodical error is joined by the arithmetical error induced by the evaluation of $p\left( z \right)$
and by the rounding of the constants $\log\left( 2 \right)$ and $\log\left( r_i \right)$.
As will be shown in section \ref{subsec:quickphase}, the overall error of the quick phase defined as
$$\epsilon_{\mbox{\tiny quick}} = \frac{\left(log_\hi + log_\lo\right) - \log\left(x\right)}{\log\left(x\right)}$$
is bounded by
$$\left \vert \epsilon_{\mbox{\tiny quick}} \right \vert \leq 5 \cdot 2^{-65} \leq 2^{-62.6}$$ ~ \par
After the computation of the quick phase double-double value $\left( log_\hi + log_\lo \right)$ a rounding test is performed
using the rounding constants according to \ref{th:roundingRN1}. If the rounding cannot be decided, the accurate
phase is launched. \par
The accurate phase performs all its computations on the same reduced argument $z = z_\hi + z_\lo$ which will be shown to be
exact. An approximation polynomial of degree $14$ is used. It is once again a modified Remez polynomial and has the
following form:
$$p\left( z \right) = z + \frac{1}{2} \cdot z + z^3 \cdot q\left( z \right)$$
where
$$q\left( z \right) =
c^\prime_3 + z \cdot \left(
c^\prime_4 + z \cdot \left(
c^\prime_5 + z \cdot \left(
c^\prime_6 + z \cdot \left(
c^\prime_7 + z \cdot \left(
c^\prime_8 + z \cdot \left(
c^\prime_9 + z \cdot r\left( z \right) \right) \right) \right) \right) \right) \right)$$
with
$c^\prime_i = c_{i\hi} + c_{i\lo} \in \F + \F$ and
$$r\left( z \right) =
c_{10} + z \cdot \left(
c_{11} + z \cdot \left(
c_{12} + z \cdot \left(
c_{13} + z \cdot c_{14} \right) \right) \right)$$
with $c_i \in \F$.
The mathematical relative error
$$\epsilon_{\mbox{\tiny meth}} = \frac{p\left( z \right) - \log\left( 1 + z \right)}{\log\left( 1 + z \right)}$$
is bounded by
$$\left \vert \epsilon_{\mbox{\tiny meth}} \right \vert \leq  2^{-125}$$
The polynomial is evaluated using double precision for $r\left( z \right)$, double-double arithmetic for
$q\left( z \right)$ and a triple-double representation for $p\left( z \right)$ and the final reconstruction.

The overall error
$$\epsilon_{\mbox{\tiny accurate}} = \frac{\left( log_\hi + log_\mi + log_\lo \right) - \log\left( x \right)}{\log\left( x \right)}$$
is bounded by
$$\left \vert \epsilon_{\mbox{\tiny accurate}} \right \vert \leq 5735 \cdot 2^{-132} \leq 2^{-119.5}$$
as will be shown in section \ref{subsec:accuratephase}. Here $\left( log_\hi + log_\mi + log_\lo \right)$
are obtained by reconstructing the logarithm as indicated by the parenthesis in the following term:
$$log_\hi + log_\mi + log_\lo = \left(E \cdot \left( log2_\hi + log2_\mi + log2_\lo \right) \right) +
\left( \left( p_\hi + p_\mi + p_\lo \right) + \left(logi_\hi + logi_\mi + logi_\lo \right) \right)$$
where
$log2_\hi + log2_\mi + log2_\lo \approx \log\left( 2 \right)$ and $logi_\hi + logi_\mi + logi_\lo \approx -\log\left(r_i \right)$.

Since the critical accuracy of the double precision $\log$ function is $118$ bits according to
\cite{DinDefLau2004LIP}, rounding $\log_\hi + log_\mi + log_\lo \approx \log\left( x \right)$ to double precision is equivalent
to rounding the infinite precision value $\log\left( x \right)$ to double precision.
Using the final rounding sequences presented in \cite{Lauter2005LIP:tripledouble}, which are supposed to be correct,
the double precision value returned by the function is the correctly rounded double precision value of
$\log\left( x \right)$.



\section{Proof of correctness of the triple-double implementation \label{sec:logtdproof}}
Proving that an implementation of an elementary function is correctly rounded means mainly proving two
bounds on the relative error $\epsilon_{\mbox{\tiny quick}}$ and $\epsilon_{\mbox{\tiny accurate}}$, using the appropriate lemma for proving the
correctness of the rounding test and conluding by means of the theorem stating the critical accuracy of the
function considered. The computation of the error bounds will be done mainly using the Gappa tool\cite{Melqu05} but
some parts of the proof are still based on paper or Maple computations. These parts will be shown in sections
\ref{subsec:reduction}, \ref{subsec:quickphase} and \ref{subsec:accuratephase} and mainly comprise the following:
\begin{itemize}
\item the demonstration that all special cases are handled correctly,
\item a proof that $z_\hi + z_\lo = r_i \cdot y - 1$ exactly,
\item the bounds for the mathematical approximation errors for the polynoms,
\item a proof of the exactness of some multiplications in the code,
\item the proof for the accuracy of all basic addition and multiplication code sequences on
double-double and triple-double numbers,
\item the correctness proof of the final rounding sequences for rounding triple-double numbers to double precision and
\item the mathematical equality of the term rewriting hints in the Gappa code.
\end{itemize}
The proofs for the accuracy of the basic operation bricks and the correctness proof of the final rounding sequences
are somewhat lengthy and are not given here; they can be found in \cite{Lauter2005LIP:tripledouble}.
\subsection{Exactness of the argument reduction\label{subsec:reduction}}
In this section, we will show that all special cases are handled correctly and that the
reduced argument consisting in $E$ and $z_\hi + z_\lo$ is exact, which means that we have the mathematically exact
equation
$$\log\left( x \right) = E \cdot \log\left( 2 \right) + \log\left( 1 + \left( z_\hi + z_\lo \right) \right) - \log\left( r_i \right)$$
This part of the algorithm is performed by the following code sequences which we will analyse line by line:
\begin{lstlisting}[caption={Handling of special cases and table access},firstnumber=1]
E=0;
xdb.d=x;

/* Filter cases */
if (xdb.i[HI] < 0x00100000){                     /* x < 2^(-1022)    */
  if (((xdb.i[HI] & 0x7fffffff)|xdb.i[LO])==0){
    return -1.0/0.0;
  }                    		                 /* log(+/-0) = -Inf */
  if (xdb.i[HI] < 0){
    return (x-x)/0;                              /* log(-x) = Nan    */
  }
  /* Subnormal number */
  E = -52;
  xdb.d *= ((db_number) ((double) two52)).d; 	 /* make x a normal number    */
}

if (xdb.i[HI] >= 0x7ff00000){
  return  x+x;				         /* Inf or Nan       */
}


/* Do argument reduction */
E += (xdb.i[HI]>>20)-1023;                       /* extract the exponent */
index = (xdb.i[HI] & 0x000fffff);
xdb.i[HI] =  index | 0x3ff00000;	         /* do exponent = 0 */
index = (index + (1<<(20-L-1))) >> (20-L);

/* reduce  such that sqrt(2)/2 < xdb.d < sqrt(2) */
if (index >= MAXINDEX){                          /* corresponds to xdb>sqrt(2)*/
  xdb.i[HI] -= 0x00100000;
  E++;
}
y = xdb.d;
index = index & INDEXMASK;

ed = (double) E;

ri = argredtable[index].ri;

logih = argredtable[index].logih;
logim = argredtable[index].logim;
\end{lstlisting}
Analysis of the code:
{\renewcommand{\labelenumi}{}
\begin{enumerate}
\item line 1 and 2: Initialization of integer $E$ and {\tt db\_number} {\tt xdb} which is now equal to $x$.
\item line 5: As the integer ordering and the ordering on floating point numbers are compatible,
$x < +2^{-1022}$, i.e. negative, negative infinite, equal to zero or a subnormal.
\item line 6: {\tt xdb.i[HI] \& 0x7fffffff} is the high order word of $x$ without the sign bit. If the test is true,
$\left \vert x \right \vert = 0$. As the logarithm of $0$ is not defined but as the limit $-\infty$ is known, returning
$-1.0 / 0.0$ is correct.
\item line 9: Since the integer ordering and the ordering on floating point numbers are compatible,
{\tt xdb.i[HI] < 0} implies $x < 0$. The logarithm is not defined for negative numbers, so the result must be $\nan$.
$0.0 / 0.0$ leads to a $\nan$; one uses $\left(x - x\right) / 0.0$ in order to overcome the static tests of the compiler.
\item line 13 and 14: if this code lines are reached, $x$ is a subnormal. Since $E$ equals $0$ at this point,
setting it to $-52$ and multipliying {\tt xdb} by $2^{-52}$ means bringing {\tt xdb} to the normal number range and
rescaling the internal representation $x = 2^E \cdot m = 2^E \cdot${\tt xdb} in consequence.
\item line 17: As the integer ordering and the ordering on floating point numbers are compatible and as
{\tt 0x7fefffff ffffffff} is the greatest normal, the test being true implies that $x$ is equal to $+\infty$ or $\nan$.
In the case of $x=+\infty$, $+\infty$ must be returned which is done. In the other case, $\nan$ must be returned which is
still assured by $x + x$.
\item line 23: At this point of the code, the most significant bit of the high order word of {\tt xdb} must be $0$ as
the case where $x < 0$ is already filtered out. So {\tt xdb.i[HI] > > 20} is equal to the biased exponent of {\tt xdb}
because a double number consists in $1$ sign bit, $11$ exponent bits and the word bit length is supposed to be $32$.
Subtracting $1023$ yields to the unbiased exponent which is written to $E$.
\item line 24 and 25: Since a double number consists in $1$ sign bit and $11$ exponent bits, the operation
{\tt xdb.i[HI] \& 0x000fffff} masks out the mantissa bits in the higher order word of {\tt xdb}.
Rewriting {\tt xdb.i[HI] = index | 0x3ff00000} means setting the exponent of {\tt xdb} to $0$ because
{\tt 0x3ff}$ - 1023 = 0$.
\item line 26: Before execution of this line of code, {\tt index} contains the high order bits of the normalized mantissa
of $x$ stored as a double in {\tt xdb.d} and verifying thus $1 \leq m < 2$. The second argument reduction step
will slice this intervall in $128$ intervalls for each of which we dispose of a table entry. For reasons of possible
cancellation in the reconstruction step on the operation $p\left( z \right) - \log\left( r_i \right)$, we want the
small intervalls to be centered around $1$. That means e.g. for the intervall around $1$ and a table indexed by $7$ bits
that mantissas (as doubles) with the high order word {\tt 0x3fefffff} through {\tt 0x3ff00fff} must be mapped to $0$.
The decision is therefore made at the $7+1$th bit of the mantissa part of the double depending on whether this bit is $0$
-- in which case the value falls in the lower intervall -- or $1$ -- in which case the value goes to the next higher
intervall. So adding $1$ to the $\left(20 - 7 - 1\right)$ rightmost bit ($L = 7$) increases the index value by $1$ iff this bit is $1$.
So after execution of the line, {\tt index} contains the number of the intervall for the second argument reduction step
centered in $1$.
\item line 29 through 31: The second adjustment to be made on $E^\prime$ and $m$ is the decision whether $m > \sqrt{2}$ as
indicated in section \ref{sec:logoutline}. The high order word of $\sqrt{2}$ rounded to a double is {\tt 0x3ff6a09e}.
As one can simply verify, the value for {\tt index} calculated for this value is $53$. As the integer ordering and
the ordering of floating point numbers are compatible and as the computations for computing {\tt index} are monotone,
{\tt index} being greater or equal than $53$ implies that the (normalized) mantissa of $x$ is greater than
$\sqrt{2} + \delta$ with a neglectable error $\delta$.
As {\tt MAXINDEX} is equal to $53$, the test will be true iff the adjustment on $E^\prime$ leading
to $E$ and $m$ yielding $y$ is to be made. It is trivial to see that the code in the {\tt if}'s body implements the
adjustment correctly.
\item lines 33 and 34: the final value of the reduced argument $y$ -- still stored in {\tt xdb.d} -- is copied to
a {\tt double} variable (or register) named {\tt y}. The final index value is masked out by means of an {\tt INDEXMASK}
which is equal to $127 = 2^7-1$.
\item lines 36: The integer value of the exponent $E$ stored in {\tt E} is cast to a {\tt double ed}.
\item lines 38 through 41: The table is indexed by {\tt index} and values {\tt ri}$=r_i$ and
{\tt logih}$=logi_\hi$ and {\tt logim}$=logi_\mi$ are read.
Since the latter form a double-double precision value, we know that
$logi_\hi + logi_\mi = \log\left( r_i \right) \cdot \left( 1 + \epsilon \right)$ with $\left \vert \epsilon \right \vert \leq 2^{-106}$.
The value {\tt ri} is stored as a single precision variable and a Maple procedure assures that for each value
$y$ the following inegality is verified:
$$\left \vert z \right \vert = \left \vert y \cdot r_i - 1 \right \vert \leq 2^{-8}$$
\end{enumerate}
}
Let us show now that the following line calculate $z_\hi$ and $z_\lo$ such that for each $y$ and corresponding $r_i$,
we obtain exactly
$$z_\hi + z_\lo = y \cdot r_i - 1$$
\begin{lstlisting}[caption={Argument reduction},firstnumber=42]
Mul12(&yrih, &yril, y, ri);
th = yrih - 1.0;
Add12Cond(zh, zl, th, yril);
\end{lstlisting}
We know that we can suppose that the multiplication and addition sequences \Mul~ and \Add~ used at lines
42 and 44 are exact. Thus, it suffices to show that
$$yri_\hi - 1.0 = yri_\hi \ominus 1.0$$
because in that case, we can note
$$z_\hi + z_\lo = th + yri_\lo = yri_\hi \ominus 1.0 + yri_\lo = y \cdot r_i - 1.0$$
We will show this property using Sterbenz' lemma. It suffices thus to prove that
$$\frac{1}{2} \leq yri_\hi \leq 2$$
We know that
\begin{eqnarray*}
yri_\hi & = & \circ\left( y \cdot r_i \right) \\
& \leq & \circ \left( 1 + 2^{-8} \right) \\
& = & 1 + 2^{-8} \\
& < & 2
\end{eqnarray*}
since the rounding function $\circ$ is monotonic and the accuracy of the format is greater than $9$ bits.

The other way round, we get
\begin{eqnarray*}
yri_\hi & = & \circ \left( y \cdot r_i \right) \\
& \geq & \circ \left( 1 - 2^{-8} \right) \\
& = & 1 - 2^{-8} \\
& > & \frac{1}{2}
\end{eqnarray*}
for the same reasons.

Thus $z_\hi + z_\lo = y \cdot r_i$ exactly. Since the previous phases of the argument reduction were all exact, the reduced argument
verifies $x = 2^{E} \cdot y$ exactly.

Still in this section, let us show that neither the reduced argument of the logarithm function nor its result may be
a sub-normal double number. The first property has already been assured by special case handling as shown above. The
latter can be proven as follows: the $\log\left( x \right)$ function has one zero for $x = 1$ and only one.
As it is monotone, for $x = 1 \pm 1 \mUlp = 1 \pm 2^{-52}$ we will obtain $\log\left( 1 \pm 2^{-52} \right) = 0 \pm 2^{-52} + \delta$
with
$\left \vert \delta \right \vert \leq 2^{-103}$. As $0 \pm 2^{-1022}$ is the least normal, the result of the logarithm function will
always be a normal. Further, in both double-double and triple-double representations for the final intermediate result
for the function, as its critical accuracy is $118$, the least significant double in the representation will still be
a normal as $52 + 106 = 158 < 1022$.
\subsection{Accuracy proof of the quick phase\label{subsec:quickphase}}
As already mentionned, the accuracy proof of the quick phase is mainly based on the Gappa tool. To prove the desired
accuracy bound defined as
$$\epsilon_{\mbox{\tiny quick}} = \frac{\left(log_\hi + log_\lo\right) - \log\left(x\right)}{\log\left(x\right)}$$
and given by
$$\left \vert \epsilon_{\mbox{\tiny quick}} \right \vert \leq 5 \cdot 2^{-65} \leq 2^{-62.6}$$
three different Gappa proof files are necessary depending on the following cases:
\begin{itemize}
\item for $E \geq 1$ and all indexes to the table $0 \leq i \leq 127$, a general proof file named {\tt log-td.gappa} is used
\item for $E = 0$ and all indexes to the table except $0$, i.e. $1 \leq i \leq 127$, a proof file named {\tt log-td-E0.gappa}
comes to hand and
\item for $E = 0$ and the table index $i = 0$, a proof file called {\tt log-td-E0-logir0.gappa} is employed.
This latter file
uses relative error computations in opposition to the other two cases where absolute error estimates suffice. This
is necessary because in this case and in this one only, the logarithm function has a zero in the intervall considered.
\end{itemize}
In each of the three proof files, we will ask the Gappa tool to verify the accuracy bound expressed in its syntax as
follows:
\begin{lstlisting}[caption={Accuracy bound to prove},firstnumber=109]
->
((logh + logm) - Log) / Log in [-5b-65,5b-65]
\end{lstlisting}
Still in any proof file, some hypothesis are made on the correctness of one multiplication sequence and the
accuracy of the constants and resting operations carried out in double-double arithmetic.
These hypothesis are the following:
\begin{itemize}
\item The operations in the following code sequence are exact since the constants are stored with enough trailing zeros:
\begin{lstlisting}[caption={Multiplication by $E$},firstnumber=50]
Add12(log2edh, log2edl, log2h * ed, log2m * ed);
\end{lstlisting}
This means that $log2ed_\hi + log2ed_\lo = E \cdot \left( log2_\hi + log2_\lo \right)$ exactly.
\item The operations in the following code sequence are exact since multiplications with a power of $2$ are exact
as long as the result is not underflowed:
\begin{lstlisting}[caption={Multiplication by $-0.5$},firstnumber=60]
zhSquareHalfh = zhSquareh * -0.5;
zhSquareHalfl = zhSquarel * -0.5;
\end{lstlisting}
i.e. $zhSquareHalf_\hi + zhSquareHalf_\lo = -0.5 \cdot \left( zhSquare_\hi + zhSquare_\lo \right)$.
\item The following hypothesis on the accuracy bounds, expressed here in Gappa syntax, are verified:
\begin{lstlisting}[caption={Gappa hypothesis},firstnumber=100]
(T2hl - T2) / T2 in [-1b-103,1b-103]
/\ (Phl - PE) / PE in [-1b-103,1b-103]
/\ (LogTabPolyhl - LogTabPoly) / LogTabPoly in [-1b-103,1b-103]
/\ (Loghm - LogE) / LogE in [-1b-103,1b-103]
/\ (Log2hm - Log2) / Log2 in [-1b-84,1b-84]
/\ (Logihm - Logir) / Logir in [-1b-106,1b-106]
/\ Z in [_zmin,_zmax]
/\ (P - Log1pZ) / Log1pZ in [-_epsilonApproxQuick,_epsilonApproxQuick]
/\ ((logh + logm) - Loghm) / Loghm in [-1b-106,1b-106]
\end{lstlisting}
Here, {\tt \_zmin}, {\tt \_zmax} and {\tt \_epsilonApproxQuick} are replaced by Maple calculated values, typically
$-zmin = zmax = 2^{-8}$ and $epsilonApproxQuick = 2^{-62.99}$.
\end{itemize}
Let us now show each of this hypothesises.
\begin{enumerate}
\item The operations yielding {\tt log2edh} and {\tt log2edl} are all exact because the \Add~ sequence is supposed
to be exact in any case and because the constants {\tt log2h} and {\tt log2m} are calculated by the following Maple
code and have in consequence at least 11 trailing zeros and {\tt ed} $=E$ is less than $1024$ in magnitude since $1024$ is
the maximum exponent value for double precision.
\begin{lstlisting}[caption={Maple code for computing {\tt log2h} and {\tt log2m}},firstnumber=21,label={list:maplelog2}]
log2acc := log(2):
log2h := round(log2acc * 2**(floor(-log[2](abs(log2acc))) + (53 - 11))) /
         2**(floor(-log[2](abs(log2acc))) + (53 - 11)):
log2m := round((log2acc - log2h) * 2**(floor(-log[2](abs((log2acc - log2h)))) +
         (53 - 11))) / 2**(floor(-log[2](abs((log2acc - log2h)))) + (53 - 11)):
\end{lstlisting}
\item To show that $zhSquareHalf_\hi + zhSquareHalf_\lo = -0.5 \cdot \left( zhSquare_\hi + zhSquare_\lo \right)$ we just have to show
that both values $zhSquare_\hi$ and $zhSquare_\lo$ are either equal to $0$ or greater than $2$ times the smallest
normal. Let us first give the definitions of both values:
\begin{eqnarray*}
zhSquare_\hi & = & \circ \left( z_\hi \cdot z_\hi \right) \\
zhSquare_\lo & = & z_\hi \cdot z_\hi - zhSquare_\hi
\end{eqnarray*}
where $z_\hi = \circ \left( z \right)$.
Let us suppose that $z \not = 0$. Otherwise all values are equal to $0$ and we can conclude.

Let us first show that $\left \vert zhSquare_\hi \right \vert$ is greater than $2^{54}$ times the smallest normal.
Let us therefore suppose that this
is not the case, i.e. $\left \vert zhSquare_\hi \right \vert < 2^{-948}$. Since the rounding function is monotonic,
this implies that $\left \vert z_\hi \right \vert \leq 2^{-424}$. For the same reason, we can note that
$\left \vert z \right \vert \leq 2^{-424}$. As we have $z = y \cdot r_i - 1$, clearly neither $y$ nor $r_i$ can be exactly $1$.
If this were the case for both, we would obtain $z=0$ which we excluded; if there were one of them only that was
exactly $1$, the other being a floating point number in the intervall $\left[ 0.5; 1.5 \right]$,
the resulting inegality $\left \vert z \right \vert \geq 2^{-53}$ which would be contradictory.

Otherwise, since we know that $1 - 2^{-8} \leq y \cdot r_i \leq 1 + 2^{-8}$ and since the precision of all formats used is greater than
$9$, the hypothesis that $1 - 2^{-424} \leq y \cdot r_i \leq 1 + 2^{-424}$ and $y \cdot r_i \not = 0$
would imply that the infinite precision mantissa of
$y \cdot r_i$ contains a $1$ weighted with $2^0$ and a $1$ weighted with less than $2^{-424}$. So its length would be greater than
$423$ bits. As it is the product of two floating point numbers which have $52$ and $23$ significant bits, there
cannot be a $1$ weighted with less than $76$ if there is a $1$ weighted with $2^0$ which is the case. Contradiction.

So $-0.5 \cdot zhSquare_\hi$ is not underflowed. Additionally, with a similar argument, since {\tt zh} is a double precision
number, $zhSquare_\lo$ is either $0$ or greater in magnitude than $2^{-53} \cdot \left \vert zhSquare_\hi \right \vert$ which is
$2^{52}$ times greater in magnitude than the smallest normal. So $zhSquare_\lo$ is either $0$ or $2$ times greater in
magnitude than the smallest normal.

So, the floating point multiplication of $zhSquare_\hi$ and $zhSquare_\lo$ with $-0.5$ can be considered to be exact.
\item {\tt (T2hl - T2) / T2 in [-1b-103,1b-103]} which means that
$$\left \vert \frac{T2hl - T2}{T2} \right \vert \leq 2^{-103}$$ is verified as $T2hl$ and $T2$ are defined as follows:
$$T2hl = t2_\hi + t2_\lo \gets \mAddDD \left( z_\hi, z_\lo, zhSquareHalf_\hi, zhSquareHalf_\lo \right)$$
$$T2 = \left( z_\hi + z_\lo \right) + \left( zhSquareHalf_\hi + zhSquareHalf_\lo \right)$$
The given bound is thus just the accuracy bound of the \AddDD~ sequence for which a proof can be found in
\cite{Lauter2005LIP:tripledouble}.
\item {\tt (Phl - PE) / PE in [-1b-103,1b-103]} is verified for the same reason; let us just recall the definitions
$$Phl = p_\hi + p_\lo \gets \mAddDD \left( t2_\hi, t2_\lo, t1_\hi, t1_\lo \right)$$
$$PE = \left( t2_\hi + t2_\lo\right) + \left( t1_\hi + t1_\lo \right)$$
\item {\tt (LogTabPolyhl - LogTabPoly) / LogTabPoly in [-1b-103,1b-103]} falls still into the same case with
$$LogTabPolyhl = logTabPoly_\hi + logTabPoly_\lo \gets \mAddDD \left( logi_\hi, logi_\mi, p_\hi, p_\lo \right)$$
$$LogTabPoly = \left( logi_\hi, + logi_\mi \right) + \left( p_\hi +  p_\lo \right)$$
\item And finally, {\tt (Loghm - LogE) / LogE in [-1b-103,1b-103]}
which is also just the accuracy bound of the \AddDD~ sequence for
$$Loghm = log_\hi + log_\mi \gets \mAddDD \left( log2ed_\hi, log2ed_\lo, logTabPoly_\hi, logTabPoly_\lo \right)$$
$$LogE = \left( log2ed_\hi + log2ed_\lo \right) + \left( logTabPoly_\hi + logTabPoly_\lo \right)$$
\item {\tt (Log2hm - Log2) / Log2 in [-1b-84,1b-84]} is verified
since $log2_\hi$ and $log2_\mi$ are computed as already indicated in listing \ref{list:maplelog2}.
This means that at least $11$ trailing zeros are stored in each in the doubles in this (pseudo-)double-double number,
so it is exact to $2^{-106-2 \cdot 11} = 2^{-84}$.
\item {\tt (Logihm - Logir) / Logir in [-1b-106,1b-106]} which means
$$\left \vert \frac{\left( logi_\hi + logi_\mi \right) - \log\left( r_i \right)}{\log\left( r_i \right)} \right \vert \leq 2^{-106}$$
is verified by construction as $logi_\hi$ and $logi_\mi$ are computed by the following Maple code:
\begin{lstlisting}[caption={Maple code for computing $logi_\hi$ and $logi_\mi$},firstnumber=35]
(logih[i], logim[i], logil[i]) := hi_mi_lo(evalf(-log(r[i]))):
\end{lstlisting}
where {\tt hi\_mi\_lo} is the procedure for rounding an arbitrary precision number to a triple-double number the higher
significant numbers of which form a double-double number.
\item The hypothesis {\tt Z in [\_zmin,\_zmax]} simply recalls the bounds for $z$ as calculated by Maple.
\item The same can be said on the hypothesis \\
{\tt (P - Log1pZ) / Log1pZ in [-\_epsilonApproxQuick,\_epsilonApproxQuick]} \\
which gives the mathematical approximation error of the polynomial. This bound is computed by Maple using the following
instructions:
\begin{lstlisting}[caption={Maple code for computing the relative error of the polynomial},firstnumber=129]
epsilonApproxQuick := numapprox[infnorm]( 1-polyQuick/log(1+x), x=zminmin..zmaxmax)
\end{lstlisting}
\item Finally, Gappa's hypothesis {\tt ((logh + logm) - Loghm) / Loghm in [-1b-106,1b-106]}
simply restates the fact that a double-double precision number is exact to
at least $2^{-106}$ in terms of its relative error.
\end{enumerate}
The Gappa tool itself is not capable of proving the final accuracy bound it is asked for a complex algorithm as the
one given here. Its user must provide hints to help it to rewrite the interval arithmetics terms it encounters
in the program. These hints are generally given in the form {\tt $\alpha$ -> $\beta$}
where $\beta$ is an expression we want the
tool to rewrite the expression $\alpha$ by. Generally speaking, the idea behind each hint is one of the following:
\begin{itemize}
\item For computing intervall bounds on differences like $\alpha = a - A$ where both $a$ and $A$ are sums of terms like
$a = c + C$ and $B = d + D$, it is often useful to rewrite $\alpha$ by $\beta = \left(c - d \right) + \left( C - D \right)$.
\item An intervall bound can often be easier found for a term $A$ representing an exact mathematical value that for $a$
which is its arithmetical equivalent. So it is useful to rewrite $a$ by $A \cdot \left( 1 + \frac{a - A}{A} \right)$ when
an intervall for $\frac{a - A}{A}$ is known.
\item Fractional left hand sides like $\frac{a}{b}$ where both expressions $a$ and $b$ are functions in a common argument
$x$ that can be written like $a = a\left( x \right) = x^n \cdot a^\prime\left( x \right)$ and
$b = b\left( x \right) = x^m \cdot b^\prime\left( x \right)$ should usually be rewritten as follows:
$$\frac{a\left(x\right)}{b\left( x \right)} = \frac{x^n \cdot a^\prime\left( x \right)}{x^m \cdot b^\prime\left( x \right)} =
x^{n - m} \cdot \frac{a^\prime\left( x \right)}{b^\prime\left( x \right)}$$ In particular, this kind of hint is needed when an
intervall for the denominator of a fractional left-hand-side comprises $0$.
\item Fractional left-hand-sides of the form $\frac{a - A}{A}$ with an unknown $A$ can easily be written like
$$\frac{a - A}{A} = \frac{a - B}{B} + \frac{B - A}{A} + \frac{a - B}{B} \cdot \frac{B - A}{A}$$
We can show this equivalence like this
\begin{eqnarray*}
\frac{a - A}{A} & = & \frac{a - B + B - A}{A} \\
& = & \frac{a - B}{A} + \frac{B - A}{A} \\
& = & \frac{a - B}{B} \cdot \frac{B}{A} + \frac{B - A}{A} \\
& = & \frac{a - B}{B} \cdot \left( 1 + \frac{B - A}{A} \right) + \frac{B - A}{A} \\
& = & \frac{a - B}{B} + \frac{B - A}{A} + \frac{a - B}{B} \cdot \frac{B - A}{A}
\end{eqnarray*}
This is particularly useful when a bound on the relative error of some term $a$ with regard to $B$ should be
extended to the next approximation level.
\end{itemize}
Clearly, the left-hand-side $A$ and right-hand-side $B$ of an hint must be mathematically equivalent to provide a
correct result. The Gappa tool checks for this equivalence and sometimes is able to prove it. If not, it emits a
warning indicating that the formal proof it is generating for the accuracy bound computations is valid only under
the hypothesis that both sides of the rewriting hint are mathematically equivalent. Further, it prints out the
difference $A - B$ of both sides $A$ and $B$ which it has already reduced using the equivalences given in the
Gappa code. It is relatively simple to verify that all this differences are equal to $0$ modulo the definitions
given in the Gappa code by means of Maple-scripts. This work can even been done automatically. Thus, we refrain
from giving a paper proof of each hint in the Gappa files used for proving the logarithm function but just
give the exhaustive list of the hints in files {\tt log-td.gappa} and {\tt log-td-E0-logir0.gappa}:
\begin{lstlisting}[caption={Gappa term rewriting hints in file {\tt log-td.gappa}},firstnumber=115]
T2hl - T2 -> ((T2hl - T2) / T2) * T2;
T2hl -> (T2hl - T2) + T2;

Phl - PE -> ((Phl - PE) / PE) * PE;
Phl -> (Phl - PE) + PE;


LogTabPolyhl -> (LogTabPolyhl - LogTabPoly) + LogTabPoly;

Loghm -> (Loghm - LogE) + LogE;

Log2 -> Log2hm * (1 / (((Log2hm - Log2) / Log2) + 1));

Logir -> Logihm * (1 / (((Logihm - Logir) / Logir) + 1));


LogTabPolyhl - LogTabPoly -> ((LogTabPolyhl - LogTabPoly) / LogTabPoly) * LogTabPoly;

HZZsimp -> (-0.5 * zh * zh) - (0.5 * zl * zl);

T2hl - ZpHZZsimp -> (0.5 * zl * zl) + delta1;

zhCube - ZZZ -> (Z * (zhSquareh - Z * Z)) - (zl * zhSquareh);

polyUpper - ZZZPhigher -> ZZZ * (polyHorner - Phigher) + polyHorner * delta3 + delta2;

ZpHZZ + ZZZPhigher -> ZpHZZsimp + ZZZPhigherPzhzl;

Phl - P -> (T2hl - ZpHZZsimp) + (T1hl - ZZZPhigherPzhzl) + delta4;

Log1pZ -> P * (1 / (((P - Log1pZ) / Log1pZ) + 1));
P - Log1pZ -> ((P - Log1pZ) / Log1pZ) * Log1pZ;

Phl - Log1pZ -> (Phl - P) + delta6;

LogTabPolyhl - Log1pZpTab -> (Logihm - Logir) + (Phl - Log1pZ) + delta7;

Loghm - Log -> (Log2edhm - Log2E) + (LogTabPolyhl - Log1pZpTab) + delta5;

(logh + logm) - Loghm -> (((logh + logm) - Loghm) / Loghm) * Loghm;

(logh + logm) - Log -> ((logh + logm) - Loghm) + (Loghm - Log);
\end{lstlisting}
\begin{lstlisting}[caption={Gappa term rewriting hints in file {\tt log-td-E0-logir0.gappa}},firstnumber=81]
T2hl - T2 -> ((T2hl - T2) / T2) * T2;
T2hl -> (T2hl - T2) + T2;

Phl - PE -> ((Phl - PE) / PE) * PE;
Phl -> (Phl - PE) + PE;


(ZhSquarehl - ZZ) / ZZ -> 2 * ((zh - Z) / Z) + ((zh - Z) / Z) * ((zh - Z) / Z);

(zhSquareh - ZZ) / ZZ -> ((ZhSquarehl - ZZ) / ZZ) + ((zhSquareh - ZhSquarehl) / ZZ);

(zhSquareh - ZhSquarehl) / ZZ -> ((zhSquareh - ZhSquarehl) / ZhSquarehl) * (ZhSquarehl / ZZ);

ZhSquarehl / ZZ -> ((ZhSquarehl - ZZ) / ZZ) + 1;

(ZhCube - ZZZ) / ZZZ -> (((zh * zhSquareh) - ZZZ) / ZZZ) + ((ZhCube - (zh * zhSquareh)) / ZZZ);

((zh * zhSquareh) - ZZZ) / ZZZ -> (1 + ((zh - Z) / Z)) * (1 + ((zhSquareh - ZZ) / ZZ)) - 1;

((ZhCube - (zh * zhSquareh)) / ZZZ) -> ((ZhCube - (zh * zhSquareh)) / (zh * zhSquareh)) * (((zh - Z) / Z) + 1) * (((zhSquareh - ZZ) / ZZ) + 1);

polyHorner / Phigher -> ((polyHorner - Phigher) / Phigher) + 1;

(polyUpper - ZZZPhigher) / ZZZPhigher -> ((polyHorner - Phigher) / Phigher) + ((ZhCube - ZZZ) / ZZZ) * (polyHorner / Phigher) +
					  + ((polyUpper - (polyHorner * ZhCube)) / (polyHorner * ZhCube)) * (polyHorner / Phigher) +
					  + ((ZhCube - ZZZ) / ZZZ) * ((polyUpper - (polyHorner * ZhCube)) / (polyHorner * ZhCube)) *
					    (polyHorner / Phigher);


((ZhSquareHalfhl - (zh * zl)) - HZZ) / HZZ -> - ((zh - Z) / Z) * ((zh - Z) / Z);

(ZhSquareHalfhl - HZZ) / HZZ -> (ZhSquarehl - ZZ) / ZZ;

((T2hl - (zh * zl)) - ZpHZZ) / ZpHZZ -> ((HZ * (((ZhSquareHalfhl - (zh * zl)) - HZZ) / HZZ)) + ((T2hl - T2) / T2)
                                        + (HZ * ((T2hl - T2) / T2))
					+ (HZ * ((ZhSquareHalfhl - HZZ) / HZZ) * ((T2hl - T2) / T2))) / (1 + HZ);

(PE - P) / P -> (((1 + HZ) * (((T2hl - (zh * zl)) - ZpHZZ) / ZpHZZ)) +
		((1 + ((zh - Z) / Z)) * (Z * ((zh - Z) / Z)) * ((Flzhzl - (zh * zl)) / (zh * zl)))
	        + (ZZ * Phigher * ((polyUpper - ZZZPhigher) / ZZZPhigher))) / (1 + HZ + ZZ * Phigher);

(Phl - P) / P -> ((PE - P) / P) + ((((PE - P) / P) + 1) * ((Phl - PE) / PE));

(Loghm - Log) / Log -> ((Loghm - P) / P) + ((P - Log) / Log) + ((Loghm - P) / P) * ((P - Log) / Log);

(((logh + logm) - Log) / Log) -> (((logh + logm) - Loghm) / Loghm) + ((Loghm - Log) / Log) + (((logh + logm) - Loghm) / Loghm) * ((Loghm - Log) / Log);
\end{lstlisting}
For the reasons mentionned, we can consider the accuracy proof of the quick phase to be correct.
\subsection{Accuracy proof of the accurate phase\label{subsec:accuratephase}}
The accuracy proof of the accurate phase is also based mainly on the use of the Gappa tool.
Nevertheless, since the tool is currently not directly supporting triple-double representations, some additional
hand-proven accuracy bound results for the main addition and multiplication operators are needed. They can be
found in \cite{Lauter2005LIP:tripledouble}. Since all these accuracy bounds are parameterized by the maximal overlap bound
for the triple-double numbers along the computations, before being able to give a numerical value for
these error bounds understood by the Gappa tool, it is necessary to do a maximal overlap bound analysis using
the theorems given in \cite{Lauter2005LIP:tripledouble}.\par
Eventually, since not an overlapped triple-double intermediate result is to be returned by the logarithm function but a
double precision number that is the correct rounding according to the rounding mode chosen, the algorithm effectuates
a renormalizing operation on the final result and rounds this non-overlapped result down to a double using an
appropriate rounding sequence. All this renormalization and rounding sequences are exact and have been shown to be
correct in \cite{Lauter2005LIP:tripledouble}. The same way, all properties shown in section \ref{subsec:reduction}
concerning the special case handling and exactness argument reduction can be reused because the algorithm implemented in
the accurate phase uses the same reduced argument and is substantially the same as for the quick phase. \par
We will thus rely on all these properties and simply show the following accuracy bound
$$\epsilon_{\mbox{\tiny accurate}} = \frac{\left( log_\hi + log_\mi + log_\lo \right) - \log\left( x \right)}{\log\left( x \right)}$$
is bounded by
$$\left \vert \epsilon_{\mbox{\tiny accurate}} \right \vert \leq 5735 \cdot 2^{-132} \leq 2^{-119.5}$$
which will be expressed in Gappa syntax as follows:
\begin{lstlisting}[caption={Accuracy bound to prove for the accurate phase},firstnumber=165]
->
((logh + logm + logl) - MLog) / MLog in [-5735b-132,5735b-132]
\end{lstlisting}
The Gappa proof files still make the hypothesis that two of the multiplications in the accurate phase code can be
considered to be
exact. This property must therefore be shown in a paper proof in the following.

The first of these multiplications is the following sequence:
\begin{lstlisting}[caption={Multiplication of triple-double $\circ\left( Z \cdot Z \right)$ by $-\frac{1}{2}$},firstnumber=99]
  zSquareHalfh = zSquareh * -0.5;
  zSquareHalfm = zSquarem * -0.5;
  zSquareHalfl = zSquarel * -0.5;
\end{lstlisting}
As it will be shown below, the relative error $\epsilon_{ZSquare}$ defined as
$$\epsilon_{ZSquare} = \frac{\left( zSquare_\hi + zSquare_\mi + zSquare_\lo \right) - Z^2}{Z^2}$$
is bounded by $\left \vert \epsilon_{ZSquare} \right \vert \leq 2^{-149}$.
Using the same argument as the one given in section \ref{subsec:quickphase}, one can show that $Z$ is either $0$
or greater in magnitude than at least $2^{-77}$. So the following is true
$$Z^2 = 0 \lor \left \vert Z^2 \right \vert \geq 2^{-154}$$
If $Z^2=0$, $ZSquarehml = zSquare_\hi + zSquare_\mi + zSquare_\lo$ trivially is $0$, too, and the multiplication is
with $-\frac{1}{2}$ is therefore exact.
Since we can note $ZSquarehml = Z^2 \cdot \left( 1 + \epsilon_{ZSquare} \right)$, we know that in the other case,
$$\left \vert ZSquarehml \right \vert \geq 2^{-155}$$
We can suppose that in the triple-double number $zSquare_\hi + zSquare_\mi + zSquare_\lo$, $zSquare_\mi$ and
$zSquare_\lo$ are not overlapped at all (since $zSquare_\mi = \circ \left( zSquare_\mi + zSquare_\lo \right)$)
and that $zSquare_\hi$ and $zSquare_\mi$ are not fully overlapped.
So we can note $\left \vert zSquare_\mi \right \vert \leq 2^{-\beta_o} \cdot \left \vert zSquare_\hi \right \vert$ and
$\left \vert zSquare_\lo \right \vert \leq 2^{-\beta_u} \cdot \left \vert zSquare_\mi \right \vert$ with $\beta_o \geq 1$ and
$\beta_u \geq 53$.
We will show this property below we are just supposing here.
So we can verify the following
\begin{eqnarray*}
\left \vert ZSquarehml \right \vert & = & \left \vert zSquare_\hi + zSquare_\mi + zSquare_\lo \right \vert \\
& \leq & \left \vert zSquare_\hi \right \vert + \left \vert zSquare_\mi \right \vert + \left \vert zSquare_\lo \right \vert \\
& \leq & \left \vert zSquare_\hi \right \vert +
2^{-\beta_o} \cdot \left \vert zSquare_\hi \right \vert +
2^{-\beta_o} \cdot 2^{-\beta_u} \cdot \left \vert zSquare_\hi \right \vert \\
& \leq & 2 \cdot \left \vert zSquare_\hi \right \vert
\end{eqnarray*}
In consequence, we obtain
$$\left \vert zSquare_\hi \right \vert \geq \frac{1}{2} \cdot \left \vert ZSquarehml \right \vert$$
and thus
$$\left \vert zSquare_\hi \right \vert \geq 2^{-156}$$ under the hypothesis that it is not exactly zero.
So $zSquareHalf_\hi = -\frac{1}{2} \cdot zSquare_\hi$ will never be underflowed.

Let us now show first that the operations for computing $zSquareHalf_\mi$ and $zSquareHalf_\lo$ cannot both be
inexact. We will use the fact that $\left \vert zSquare_\lo \right \vert \leq 2^{-53} \cdot \left \vert zSquare_\mi \right \vert$.
Suppose first that
$$zSquareHalf_\mi \gets - \frac{1}{2} \otimes zSquare_\mi$$ is inexact. So $\left \vert zSquare_\mi \right \vert < 2^{-1022}$ and
in consequence $\left \vert zSquare_\lo \right \vert < 2^{-1022-53}$. Note that the inegality is strict.
Since the least (in magnitude) representable denormalized double precision floating point number is $2^{-52} \cdot 2^{-1023}$,
$zSquare_\lo = 0$ in this case. So $$zSquareHalf_\lo \gets - \frac{1}{2} \otimes zSquare_\lo$$ is exact because trivially, a
multiplication with $0$ is exact.

Suppose now that $$zSquareHalf_\lo \gets - \frac{1}{2} \otimes zSquare_\lo$$ is inexact.
So $\left \vert zSquare_\lo \right \vert < 2^{-1022}$. Further, the least significant bit of the mantissa of $zSquare_\lo$ is
$1$ because otherwise, a bit-shift in its mantissa by 1 would be an exact operation.
Thus $\left \vert zSquare_\lo \right \vert \geq 2^{-52} \cdot 2^{-1023}$ and $\left \vert zSquare_\mi \right \vert \geq 2^{-1022}$. So
$$zSquareHalf_\mi \gets - \frac{1}{2} \otimes zSquare_\mi$$ cannot be inexact because in this case we would have
$\left \vert zSquare_\mi \right \vert < 2^{-1022}$.

So, in any case, if ever $zSquareHalf_\mi + zSquareHalf_\lo$ are not exactly
$-\frac{1}{2} \cdot \left( zSquare_\mi + zSquare_\lo \right)$, the error made will be $\frac{1}{2} \cdot d$
in magnitude, where $d = 0^+$ is the smallest representable denormalized non-zero double. So we can note down in this case
$$zSquareHalf_\hi + zSquareHalf_\mi + zSquareHalf_\lo = - \frac{1}{2} \cdot \left( zSquare_\hi + zSquare_\mi + zSquare_\lo \right) + \delta$$
with $\left \vert \delta \right \vert \leq 2^{-1075}$. Since we know that
$\left \vert - \frac{1}{2} \cdot \left( zSquare_\hi + zSquare_\mi + zSquare_\lo \right) \right \vert \geq 2^{-156}$, we can give the
following bound
$$\left \vert \frac{\delta}{-\frac{1}{2} \cdot \left( zSquare_\hi + zSquare_\mi + zSquare_\lo \right)} \right \vert \leq
\frac{2^{-1075}}{2^{-156}} = 2^{-919}$$
So we get
$$ZSquareHalfhml = - \frac{1}{2} \cdot ZSquarehml \cdot \left(1 + \epsilon\right)$$
with
$\left \vert \epsilon \right \vert \leq 2^{-919}$

In contrast, since we know that $\left \vert Z \right \vert \leq 2^{-8}$
thus that $\left \vert Z^2 \right \vert \leq 2^{-16}$ but that $\left \vert Z^2 \right \vert \geq 2^{-154}$, we can
assume that the infinite precision mantissa of $Z^2$ can always be written exactly with at most $154 - 16 = 138 < 149$
bits. As we can show that
$\frac{1}{2} \cdot \left \vert ZSquarehml \right \vert \leq \left \vert zSquare_\hi \right \vert \leq
2 \cdot \left \vert ZSquarehml \right \vert$ we know that if ever one of $zSquare_\mi$ or $zSquare_\lo$ is such that
the multiplication with $-\frac{1}{2}$ is not exact, the error made has already been accounted for in the error bound
for $ZSquarehml$ with regard to $Z^2$.
So the operation computing $ZSquareHalfhml$ out of $ZSquarehml$ can be considered to be exact. \par
Let us now analyse the following sequence
\begin{lstlisting}[caption={Multiplication of triple-double $Log2hml$ by $E$},firstnumber=126]
  log2edhover = log2h * ed;
  log2edmover = log2m * ed;
  log2edlover = log2l * ed;
\end{lstlisting}
Similar to the argumentation that has been given in section \ref{subsec:quickphase}, since $E=ed$ is bound
in magnitude by $1024=2^{10}$ and since $log2_\hi$, $log2_\mi$ are stored with at least $11$ trailing bits at zero,
the multiplications in these components are exact. The constant $log2_\lo$ is not stored with $11$ trailing bits at zero
but it could be because we will be just supposing the bound $\left \vert \epsilon_{Log2hml} \right \vert \leq 2^{-3 \cdot 53 + 33} = 2^{-126}$ for
$$\epsilon_{Log2hml} = \frac{log2_\hi + log2_\mi + log2_\lo - \log\left( 2 \right)}{\log\left( 2 \right)}$$
So the multiplication is not exact in itself but the final result is exacter than the bound we are using for it.

Let us finally just recall the Maple code for computing the constants:
\begin{lstlisting}[caption={Maple code for computing $Log2hml$},firstnumber=21]
log2acc := log(2):
log2h := round(log2acc * 2**(floor(-log[2](abs(log2acc))) + (53 - 11))) / 2**(floor(-log[2](abs(log2acc))) + (53 - 11)):
log2m := round((log2acc - log2h) * 2**(floor(-log[2](abs((log2acc - log2h)))) + (53 - 11))) / 2**(floor(-log[2](abs((log2acc - log2h)))) + (53 - 11)):
log2l := log2acc - (log2h + log2m):
\end{lstlisting}
So the multiplication can be considered to be exact as long the less accurate bound for $\epsilon_{Log2hml}$ is used.

Let us know analyse the bounds that we can give for the maximal overlap of the components of the triple-double numbers
in the logarithm implementation. For doing this, we will assign each triple-double number in the code an overlap bound as
follows. Call the number in consideration e.g. $a_\hi + a_\mi + a_\lo$. So we will give the bounds expressed like this:
\begin{eqnarray*}
\left \vert a_\mi \right \vert & \leq & 2^{-\alpha_o} \cdot \left \vert a_\hi \right \vert \\
\left \vert a_\lo \right \vert & \leq & 2^{-\alpha_u} \cdot \left \vert a_\mi \right \vert
\end{eqnarray*}
where $\alpha_o, \alpha_u \geq 2$.
We will then propagate this information following the flow of control in the implementation and using the overlap
bound theorems given in \cite{Lauter2005LIP:tripledouble}. Here, we understand by ``propagating'' checking a system of constraints
of the bounds under the limitations provided by the theorems. As the control-flow-graph of our implementation
is completely linear, this check is linear, too. The theorems mentionned can be summarized as follows:
\begin{center}
\begin{tabular}{|l|ll|ll|ll|}
\hline
Operation & 1st arg. & 2nd arg. & result high & result low \\
\hline
\AddTT & $\alpha_o \geq 4$, $\alpha_u \geq 1$ & $\beta_o \geq 4$, $\beta_u \geq 1$ & $\gamma_o \geq \min\left( \alpha_o, \beta_o \right) - 5$ & $\gamma_u \geq 53$ \\
\hline
\AddDTT & - & $\beta_o \geq 2$, $\beta_u \geq 1$ & $\gamma_o \geq \min\left( 45, \beta_o - 4, \beta_o + \beta_u - 2 \right)$ & $\gamma_u \geq 53$ \\
\hline
\MulDT & - & - & $\gamma_o \geq 48$ & $\gamma_u \geq 53$ \\
\hline
\MulDTT & - & $\beta_o \geq 2$, $\beta_u \geq 1$ & $\gamma_o \geq \min\left( 48, \beta_o - 4, \beta_o + \beta_u - 4 \right)$ & $\gamma_u \geq 53$ \\
\hline
\end{tabular}
\end{center}
So let us analyse the following code:
\begin{lstlisting}[caption={Triple-double computations},firstnumber=90,label={list:tripledouble}]
Mul23(&zSquareh, &zSquarem, &zSquarel, zh, zl, zh, zl);
Mul233(&zCubeh, &zCubem, &zCubel, zh, zl, zSquareh, zSquarem, zSquarel);
Mul233(&higherPolyMultZh, &higherPolyMultZm, &higherPolyMultZl, t14h, t14l, zCubeh, zCubem, zCubel);
zSquareHalfh = zSquareh * -0.5;
zSquareHalfm = zSquarem * -0.5;
zSquareHalfl = zSquarel * -0.5;
Add33(&polyWithSquareh, &polyWithSquarem, &polyWithSquarel,
      zSquareHalfh, zSquareHalfm, zSquareHalfl,
      higherPolyMultZh, higherPolyMultZm, higherPolyMultZl);
Add233(&polyh, &polym, &polyl, zh, zl, polyWithSquareh, polyWithSquarem, polyWithSquarel);
Add33(&logyh, &logym, &logyl, logih, logim, logil, polyh, polym, polyl);
log2edhover = log2h * ed;
log2edmover = log2m * ed;
log2edlover = log2l * ed;
log2edh = log2edhover;
log2edm = log2edmover;
log2edl = log2edlover;
Add33(&loghover, &logmover, &loglover, log2edh, log2edm, log2edl, logyh, logym, logyl);
\end{lstlisting}
This code will finally generate triple-double numbers respecting the following overlap bounds as will be
shown below:
\begin{center}
\begin{tabular}{|l|l|l|l|}
\hline
Variable & Line(s) & $\alpha_o \geq$ & $\alpha_u \geq$ \\
\hline
$ZSquarehml$ & 90 & $48$ & $53$ \\
\hline
$ZCubehml$ & 91 & $44$ & $53$ \\
\hline
$HigherPolyMultZhml$ & 92 & $40$ & $53$ \\
\hline
$ZSquareHalfhml$ & 93-95 & $48$ & $53$ \\
\hline
$PolyWithSquarehml$ & 96-98 & $35$ & $53$ \\
\hline
$Polyhml$ & 99 & $31$ & $53$ \\
\hline
$Logyhml$ & 100 & $26$ & $53$ \\
\hline
$Log2edhml$ & 101-106 & $40$ & $40$ \\
\hline
$Logoverhml$ & 107 & $21$ & $53$ \\
\hline
\end{tabular}
\end{center}
So let us verify exemplarily some of these bounds:
\begin{itemize}
\item At line 90, $ZSquarehml$ is computed out of the double-double number $z_\hi + z_\lo$ by use of the \MulDT~ sequence.
Since the inputs of this function are not triple-double, the overlap bound is just the bound provided by the sequence
itself, i.e. $\alpha_o \geq 48$, $\alpha_u \geq 53$.
\item $ZCubehml$ is the result of a \MulDTT~ sequence at line 91. Its overlap bound depends therefore on the one
for $ZSquarehml$, which is the second argument of the function. Since we know the bound for this variable, we easily
verify the one for $ZCubehml$ which is $\alpha_o \geq 44$ and $\alpha_u \geq 53$.
\item $Log2edhml$ is the exact pairwise product of the triple-double constant $Log2hml$ and double $E$. Since $E$ may be
as small as $0$ in magnitude and further, since the multiplication is pairwise, the overlap bound we dispose of for
$Log2edhml$ is the same as for $Log2hml$ which is stored with at least $11$ bit trailing zeros.
So an appropriate bound is $\alpha_o \geq 52 - 11 \geq 40$ and $\alpha_u \geq 40$.
\end{itemize}
All other bounds can be verified the same way using the theorems given in \cite{Lauter2005LIP:tripledouble} and indicated above. \par
Since we have computed the overlap bounds for the different triple-double operands in the code, we can now
calculate the accuracy bounds for the operations. Doing this is only possible with the knowledge of the
overlap of the operations because all accuracy bound theorems given in \cite{Lauter2005LIP:tripledouble} are parameterized with this
overlap expressions.

Let us first give a list of the accuracy of the different basic operations which is not exhaustive with regard to
its lack of listing almost all preconditions on the sequences required for theorems to hold. We refrain from explicitely
verifying each of this preconditions in this document as this is only fastidious work but not of special interest.
\begin{center}
\begin{tabular}{|l|l|l|l|}
\hline
Operation & Overlap 1st arg. & Overlap 2nd arg. & Relative error $\epsilon$ \\
\hline
\AddDD & - & - & $\left \vert \epsilon \right \vert \leq 2^{-103.5} \leq 2^{-103}$ \\
\hline
\MulDD & - & - & $\left \vert \epsilon \right \vert \leq 2^{-102}$\\
\hline
\AddTT & $\alpha_o \geq 4$, $\alpha_u \geq 1$ & $\beta_o \geq 4$, $\beta_u \geq 1$ &
$\left \vert \epsilon \right \vert \leq 2^{-\min\left( \alpha_o + \alpha_u, \beta_o + \beta_u \right) -47} + 2^{-\min\left( \alpha_o, \beta_o \right) - 98}$
\\
\hline
\AddDTT & - & $\beta_o \geq 2$, $\beta_u \geq 1$ &
$\left \vert \epsilon \right \vert \leq 2^{-\beta_o - \beta_u - 52} + 2^{-\beta_o-104} + 2^{-153}$
\\
\hline
\MulDT & - & - &
$\left \vert \epsilon \right \vert \leq 2^{-149}$
\\
\hline
\MulDTT & - & $\beta_o \geq 2$, $\beta_u \geq 1$ &
$\left \vert \epsilon \right \vert \leq 2^{-97-\beta_o} + 2^{-97-\beta_o-\beta_u} + 2^{-150}$
\\
\hline
\end{tabular}
\end{center}
Still analyzing the following double-double computations code and the code
given at listing \ref{list:tripledouble}, one can now easily check the bounds for
the relative error of the different operations listed in the table below.
We define here the relative error of an operation $\ast$ and its arithmetical equivalent $\circledast$ as follows:
$$\epsilon = \frac{\left(a \circledast b \right) - \left(a \ast b\right)}{\left(a \ast b \right)}$$
\begin{lstlisting}[caption={Double-double computations in accurate phase},firstnumber=73,label={list:doubledouble}]
  Mul12(&t1h, &t1l, zh, highPoly);
  Add22(&t2h, &t2l, accPolyC9h, accPolyC9l, t1h, t1l);
  Mul22(&t3h, &t3l, zh, zl, t2h, t2l);
  Add22(&t4h, &t4l, accPolyC8h, accPolyC8l, t3h, t3l);
  Mul22(&t5h, &t5l, zh, zl, t4h, t4l);
  Add22(&t6h, &t6l, accPolyC7h, accPolyC7l, t5h, t5l);
  Mul22(&t7h, &t7l, zh, zl, t6h, t6l);
  Add22(&t8h, &t8l, accPolyC6h, accPolyC6l, t7h, t7l);
  Mul22(&t9h, &t9l, zh, zl, t8h, t8l);
  Add22(&t10h, &t10l, accPolyC5h, accPolyC5l, t9h, t9l);
  Mul22(&t11h, &t11l, zh, zl, t10h, t10l);
  Add22(&t12h, &t12l, accPolyC4h, accPolyC4l, t11h, t11l);
  Mul22(&t13h, &t13l, zh, zl, t12h, t12l);
  Add22(&t14h, &t14l, accPolyC3h, accPolyC3l, t13h, t13l);
\end{lstlisting}
\begin{center}
\begin{tabular}{|l|l|l|l|}
\hline
Result & Line(s) & Operation & Relative error $\epsilon$ \\
\hline
$T1hl$ through $T14hl$ & 73 - 86 & \AddDD~ / \MulDD &
$\left \vert \epsilon \right \vert \leq 2^{-103}$ / $\left \vert \epsilon \right \vert \leq 2^{-102}$ \\
\hline
$ZSquarehml$ & 90 & \MulDT & $\left \vert \epsilon \right \vert \leq 2^{-149}$ \\
\hline
$ZCubehml$ & 91 & \MulDTT & $\left \vert \epsilon \right \vert \leq 2^{-144}$ \\
\hline
$HigherPolyMultZhml$ & 92 & \MulDTT & $\left \vert \epsilon \right \vert \leq 2^{-141}$ \\
\hline
$PolyWithSquarehml$ & 96-98 & \AddTT & $\left \vert \epsilon \right \vert \leq 2^{-137}$ \\
\hline
$Polyhml$ & 99 & \AddDTT & $\left \vert \epsilon \right \vert \leq 2^{-134}$ \\
\hline
$Logyhml$ & 100 & \AddTT & $\left \vert \epsilon \right \vert \leq 2^{-128}$ \\
\hline
$Logoverhml$ & 107 & \AddTT & $\left \vert \epsilon \right \vert \leq 2^{-123}$ \\
\hline
\end{tabular}
\end{center}
Let us just explicitely check the bound for one of the operations for sake of an example. Let us take
therefore the \AddTT~ sequence at lines 96-98 computing $PolyWithSquarehml$ out of $ZSquareHalfhml$ and
$HigherPolyMultZhml$. We have already obtained to following overlap bounds:
\begin{eqnarray*}
\left \vert zSquareHalf_\mi \right \vert & \leq & 2^{-48} \cdot \left \vert zSquareHalf_\hi \right \vert \\
\left \vert zSquareHalf_\lo \right \vert & \leq & 2^{-53} \cdot \left \vert zSquareHalf_\mi \right \vert \\
\left \vert higherPolyMultZ_\mi \right \vert & \leq & 2^{-40} \cdot \left \vert higherPolyMultZ_\hi \right \vert \\
\left \vert higherPolyMultZ_\lo \right \vert & \leq & 2^{-53} \cdot \left \vert higherPolyMultZ_\mi \right \vert
\end{eqnarray*}
Feeding now this bounds into the theorem on the accuracy of \AddTT, we get
$$\left \vert \epsilon \right \vert \leq 2^{-\min\left( 48 + 53, 40 + 53 \right) - 47} + 2^{-\min \left( 48, 40 \right) - 98} \leq 2^{-140} + 2^{-138} \leq 2^{-137}$$
All other error bounds can be verified in a similar way. They are finally expressed in Gappa syntax as follows:
\begin{lstlisting}[caption={Relative error bounds in Gappa code},firstnumber=139]
(T2hl - T2) / T2 in [-1b-103,1b-103]
/\ (T3hl - T3) / T3 in [-1b-102,1b-102]
/\ (T4hl - T4) / T4 in [-1b-103,1b-103]
/\ (T5hl - T5) / T5 in [-1b-102,1b-102]
/\ (T6hl - T6) / T6 in [-1b-103,1b-103]
/\ (T7hl - T7) / T7 in [-1b-102,1b-102]
/\ (T8hl - T8) / T8 in [-1b-103,1b-103]
/\ (T9hl - T9) / T9 in [-1b-102,1b-102]
/\ (T10hl - T10) / T10 in [-1b-103,1b-103]
/\ (T11hl - T11) / T11 in [-1b-102,1b-102]
/\ (T12hl - T12) / T12 in [-1b-103,1b-103]
/\ (T13hl - T13) / T13 in [-1b-102,1b-102]
/\ (T14hl - T14) / T14 in [-1b-103,1b-103]
/\ (ZSquarehml - ZSquare) / ZSquare in [-1b-149,1b-149]
/\ (ZCubehml - ZCube) / ZCube in [-1b-144,1b-144]
/\ (HigherPolyMultZhml - HigherPolyMultZ) / HigherPolyMultZ in [-1b-141,1b-141]
/\ (PolyWithSquarehml - PolyWithSquare) / PolyWithSquare in [-1b-137,1b-137]
/\ (Polyhml - Poly) / Poly in [-1b-134,1b-134]
/\ (Logyhml - Logy) / Logy in [-1b-128,1b-128]
/\ (Loghml - Logover) / Logover in [-1b-123,1b-123]
/\ (Log2hml - MLog2) / MLog2 in [-1b-126,1b-126]
/\ (Logihml - MLogi) / MLogi in [-1b-159,1b-159]
/\ (MPoly - MLog1pZ) / MLog1pZ in [-_epsilonApproxAccurate,_epsilonApproxAccurate]
/\ Z in [_zmin,_zmax]
/\ ((logh + logm + logl) - Loghml) / Loghml in [-1b-159,1b-159]
\end{lstlisting}
Concerning the Gappa proofs for accurate phase, in a similar way as for the quick phase, three different proof
files are used. They reflect once again the three main cases for the argument of the logarithm function:
\begin{itemize}
\item For cases where after argument reduction $\left \vert E \right \vert \geq 1$,
the file {\tt log-td-accurate.gappa} is used. In this case, absolute error computations are sufficient
for the final relative error bound to be calculable because
$\left \vert \log\left( x \right) \right \vert \geq \frac{1}{2} \log\left( 2 \right)$ in this case.
\item For the case where after argument reduction, $E = 0$ and $i \not = 0$, the file
{\tt log-td-accurate-E0.gappa} is used. The same way here, we have a preponderant constant term so absolute
error computations suffice.
\item For the other case, where $E=0$ and $i=0$ the file {\tt log-td-accurate-E0-logir0.gappa} provides the
accuracy bound proof. In contrast to the other cases, we obliged to relative error estimations since the beginning
since the function $\log\left( x \right)$ has a zero in this intervall.
\end{itemize}
Once again, several term rewriting hints are needed in the Gappa proof files for enabling the Gappa tool to
generate a proof for the accuracy bounds. In a similar way, the hints which cannot directly be checked for their
mathematical correctness by the tool itself are verified by semi-automatic Maple scripts.\par
By the existence of an accuracy proof for a final relative error of $\left \vert \epsilon_{\mbox{\tiny accurate}} \right \vert \leq 2^{-119.5}$ and
by the use of the critical accuracy of the double precision natural logarithm function which is
$118$ bits\cite{DinDefLau2004LIP}, we can consider the implementation to be correctly rounding under the hypothesis
that the final rounding sequence used is exact and correct. Since we suppose this -- a correctness proof can be
found in \cite{Lauter2005LIP:tripledouble} -- the correctly rounding property is verified.


\section{Proof of correctness of the double-extended implementation \label{sec:logdeproof}}


\section{Performance results\label{sec:logperf}}
The given implementation of the natural logarithm function aims at
being both portable and more performant than the previous
implementations using the SCS libary for the accurate phase.  This
goal is acheived in terms of memory consumption (if the code sizes for
{\tt scslib} are taken into account) and in terms of speed
performance. The reason for this is mainly the possibility of reusing
all values computed in the argument reduction phase and the tables for
the accurate phase directly.

\subsection{Memory requirements}



\subsection{Timings}


%%% Local Variables:
%%% mode: latex
%%% TeX-master: "crlibm"
%%% End: