$|\circ(r_h \cdot \text{INVLN2})| \le 709.79 \cdot (2^{12}/\log 2 + 2^{-43.447}) \cdot (1 + 2^{-52}) < 4194347.07$, so the integer $k$ nearest $\circ(r_h \cdot \text{INVLN2})$ satisfies $|k| \le 4194347$.

Furthermore, $|k - (r_h + r_\ell) \cdot 2^{12}/\log 2| \le D_1 + D_2 + D_3 + D_4$ with $D_1 := |k - \circ(r_h \text{INVLN2})|$, $D_2 := |\circ(r_h \text{INVLN2}) - r_h \text{INVLN2}|$, $D_3 := |(r_h + r_\ell)(\text{INVLN2} - 2^{12}/\log 2)|$, and $D_4 := |r_\ell \text{INVLN2}|$. Clearly, $D_1 \le 1/2$. In addition, recalling from Lemma 5 that $|r_h| \le 709.79$, $|r_\ell| \le 2^{-14.4187}$, and $|r_h + r_\ell| \le 709.79$, we deduce that $|r_h \text{INVLN2}| < 2^{22.1}$, $D_2 \le \text{ulp}(2^{22.1}) \le 2^{-30}$, $D_3 \le 709.79 \cdot 2^{-43.447} \le 2^{-33.975}$, and $D_4 \le 2^{-14.4187}$. $\text{INVLN2} \le 0.2698195$. All this gives $|k - (r_h + r_\ell) \cdot 2^{12}/\log 2| \le 0.7698196$ and, therefore,

$$|k \log 2 \cdot 2^{-12} - (r_h + r_\ell)| \le 2^{-12.906174}.$$

**The subtraction $\circ(r_h - k\text{LN2H})$ is exact.** Recalling from Lemma 5 that $|r_h| \in [2^{-970}, 709.79] \cap \mathbb{F}$ and since $\text{LN2H} \in (0, 1) \cap 2^{-65}\mathbb{Z}$ and $|k| \le 4194347$, we have $\Omega \ge |r_h - k\text{LN2H}| \in 2^{-1022}\mathbb{Z} \subset \alpha\mathbb{Z}$. Hence there is no underflow/overflow when computing $z_h = \circ(r_h - k\text{LN2H})$.

Let us now show that this operation is exact. We have $|k - r_h/\text{LN2H}| \le D_1 + D_2 + D'_3$ with $D_1$, $D_2$ as before, and $D'_3 := |r_h(\text{INVLN2} - 1/\text{LN2H})|$. As before, $D_1 \le 1/2$ and $D_2 \le 2^{-30}$. In addition, one can check that $|\text{INVLN2} - 1/\text{LN2H}| < 2^{-41.694}$, which leads to $D'_3 \le 709.79 \cdot 2^{-41.694} \le 2^{-32.222}$. Hence $|r_h - k\text{LN2H}| \le (D_1 + D_2 + D'_3)\text{LN2H} \le (1/2 + 2^{-30} + 2^{-32.222}) \cdot \text{LN2H} \le 2^{-13.528766}$.

Assume first that $|r_h| \ge 2^{-13}$. Then $r_h \in 2^{-65}\mathbb{Z}$ and so $r_h - k\text{LN2H} \in 2^{-65}\mathbb{Z}$. Hence $r_h - k\text{LN2H} = m \cdot 2^{-65}$ for some integer $m$ such that $|m| \le 2^{65} \cdot 2^{-13.528766} \le 2^{53}$, and thus $r_h - k\text{LN2H}$ is in $\mathbb{F}$.

Assume now that $|r_h| < 2^{-13}$. Then $|r_h \cdot \text{INVLN2}| < 0.73$ and so $k := \lfloor \circ(r_h \cdot \text{INVLN2}) \rfloor \in \{-1, 0, 1\}$. If $k=0$, then clearly $r_h - k\text{LN2H} \in \mathbb{F}$. If $k=1$, then $0 \le r_h < 2^{-13} < \text{LN2H} = 2^{-12.52...}$ and, on the other hand, $1/2 \le \circ(r_h \cdot \text{INVLN2})$. Hence $r_h \cdot \text{INVLN2} > 1/2 - 2^{-54}$ (the predecessor of $1/2$ in $\mathbb{F}$) and thus $r_h \ge \text{RU}((1/2 - 2^{-54})/\text{INVLN2}) = \text{LN2H}/2$. It follows from Sterbenz' theorem that $r_h - \text{LN2H}$ is in $\mathbb{F}$. Similarly, if $k=-1$ then $-2^{-13} < r_h \le 0$ and $\circ(r_h \cdot \text{INVLN2}) \le -1/2$, which requires $r_h \cdot \text{INVLN2} < -1/2 + 2^{-54}$ and so $r_h \le -\text{LN2H}/2$ and $r_h - k\text{LN2H} = -|r_h| + \text{LN2H} \in \mathbb{F}$ by Sterbenz' theorem.

We have thus shown that the value $z_h \in \mathbb{F}$ satisfies

$$z_h = r_h - k\text{LN2H}.$$

**About the value $z_\ell$.** Since $|r_\ell| < 2^{-14.4187}$ (from Lemma 5) and $|k| \le 4194347$, we have $|r_\ell - k\text{LN2L}| \le 2^{-14.4187} + 4194347 \cdot \text{LN2L} < 2^{-14.418}$, so overflow cannot occur when computing $z_\ell := \circ(r_\ell - k\text{LN2L})$. If underflow occurs, then $|z_\ell - (r_\ell - k\text{LN2L})| \le 2^{-1074}$; otherwise, $|z_\ell - (r_\ell - k\text{LN2L})| \le \text{ulp}(r_\ell - k\text{LN2L}) \le \text{ulp}(2^{-14.418}) = 2^{-67}$. Hence

in both cases the rounding error when computing $z_\ell$ satisfies

$$|z_\ell - (r_\ell - k\text{LN2L})| \le 2^{-67}.$$

**About the value $z$.** We have $|z_h + z_\ell| \le |z_h| + |z_\ell| = |r_h - k\text{LN2H}| + |z_\ell| \le 2^{-13.528766} + 2^{-14.418} + 2^{-67} \le 2^{-12.9059}$. Hence overflow cannot occur when computing $z := \circ(z_h + z_\ell)$. Furthermore, $|z - (z_h + z_\ell)|$ is at most $2^{-1074}$ if underflow occurs, and at most $\text{ulp}(2^{-12.9059}) = 2^{-65}$ otherwise. Hence in all cases

$$|z - (z_h + z_\ell)| \le 2^{-65}.$$

Since $|k(\text{LN2H} + \text{LN2L}) - k \log 2 \cdot 2^{-12}| \le 4194347 \cdot 2^{-122.43} \le 2^{-100.429}$, we deduce that

$$|z - (r_h + r_\ell - k \log 2 \cdot 2^{-12})| \le 2^{-65} + 2^{-67} + 2^{-100.429} \le 2^{-64.67807}.$$

Therefore, $r_h + r_\ell = k \log 2 \cdot 2^{-12} + z + \epsilon$ with $|\epsilon| \le 2^{-64.67807}$, which implies further $\exp(r_h + r_\ell) = 2^{k/2^{12}} \cdot \exp(z) \cdot \exp(\epsilon)$. Defining $\epsilon_1$ by $\exp(\epsilon) =: 1 + \epsilon_1$, we have $|\epsilon_1| \le 2^{-64.67806}$, that is,

$$\left| \frac{\exp(r_h + r_\ell)}{2^{k/2^{12}} \cdot \exp(z)} - 1 \right| \le 2^{-64.67806}.$$

Finally, recalling that $|k \log 2 \cdot 2^{-12} - (r_h + r_\ell)| \le 2^{-12.906174}$, we arrive at $|z| \le 2^{-64.67807} + 2^{-12.906174} \le 2^{-12.905}$.

**About the values $h_1$, $\ell_1$, $h_2$, and $\ell_2$.** Given $i_1, i_2 \in \{0, \dots, 63\}$, consider the numbers $2^{i_1/2^{12}}$ and $2^{i_2/64}$. These numbers are in $[1, 2)$ and their ulp is thus $2^{-52}$.

Let $h_1$ be the binary64 number nearest $2^{i_1/2^{12}}$, and let $\ell_1$ be the binary64 number nearest $2^{i_1/2^{12}} - h_1$. (These two numbers are normal binary64 numbers.) We have $h_1 = 2^{i_1/2^{12}} + \epsilon_1$ with $|\epsilon_1| \le 2^{-53}$, and $\ell_1$ is the binary64 number nearest $-\epsilon_1$, in particular, $|\ell_1| \le 2^{-53}$. If $|\epsilon_1| = 2^{-53}$, then $\ell_1 = -\epsilon_1$ and thus $h_1 + \ell_1 = 2^{i_1/2^{12}}$; if $|\epsilon_1| < 2^{-53}$, then $\ell_1 = -\epsilon_1 + \epsilon'_1$ with $|\epsilon'_1| \le 0.5\text{ulp}(\epsilon_1) \le 2^{-107}$ and it follows that $|h_1 + \ell_1 - 2^{i_1/2^{12}}| \le 2^{-107}$. If $i_1 = 0$, then $2^{i_1/2^{12}} = 1$, $h_1 = 1$, and $\ell_1 = 0$, so $h_1 + \ell_1$ is exactly $2^{i_1/2^{12}}$; if $i_1 \ge 1$, then $2^{i_1/2^{12}} \ge 2^{1/2^{12}} > 2^{0.000244}$. Hence in both cases, the relative error of $h_1 + \ell_1$ satisfies

$$\left| \frac{h_1 + \ell_1}{2^{i_1/2^{12}}} - 1 \right| < 2^{-107.000244}.$$

Note also that an exhaustive search over the 64 possible values gives a relative error always less than $2^{-107.0228}$. Furthermore, $h_1 \in \mathbb{F} \cap [1, 2)$ and $\ell_1 \in \mathbb{F}$ is either 0 or such that $|\ell_1| > 2^{-57.53} \ge 2^{-58}$, so that $h_1 \in 2^{-52}\mathbb{Z}$ and $\ell_1 \in 2^{-110}\mathbb{Z}$.

Similarly, with $h_2$ the binary64 number nearest $2^{i_2/64}$ and $\ell_2$ the binary64 number nearest $h_2 - 2^{i_2/64}$, one can check that $|\ell_2| \le 2^{-53}$ and (using now $2^{1/64} = 2^{0.015625}$) that

$$\left| \frac{h_2 + \ell_2}{2^{i_2/64}} - 1 \right| \le 2^{-107.015625}.$$
Note also that an exhaustive search over the 64 possible values gives a relative error always less than $2^{-107.57149}$. Furthermore, $h_2 \in \mathbb{F} \cap [1, 2)$ and $\ell_2 \in \mathbb{F}$ is either 0 or such that $|\ell_2| > 2^{-58.98} \ge 2^{-59}$, so that $h_2 \in 2^{-52}\mathbb{Z}$ and $\ell_2 \in 2^{-111}\mathbb{Z}$.

Finally, by definition, $1 \le h_1 \le 2^{63/2^{12}} \cdot (1 + 2^{-53}) < 2^{0.015381}$ and $1 \le h_2 \le 2^{63/64} \cdot (1 + 2^{-53}) < 2^{0.984376}$, which implies that the exact product $h_1h_2$ satisfies $1 \le h_1h_2 < 2^{0.999757} < 2$.

**About the values $p_h$ and $p_\ell$.** Since $h_1, h_2 \in \mathbb{F} \cap [1, 2)$, there is no underflow/overflow when computing $p_h := \circ(h_1h_2)$, and $h_1h_2 - p_h \in \mathbb{F}$. Hence, recalling that $1 \le h_1h_2 < 2$, we have $p_h = h_1h_2 + \epsilon_1$ with $|\epsilon_1| \le \text{ulp}(h_1h_2) \le 2^{-52}$. Since $h_1, h_2, p_h \in 2^{-52}\mathbb{Z}$, we have furthermore $|h_1h_2 - p_h| \in 2^{-104}\mathbb{Z}$, so there is no underflow/overflow when computing $s := \circ(h_1h_2 - p_h)$, and we thus have $s = h_1h_2 - p_h$ and $2^{-52} \ge |s| \in 2^{-104}\mathbb{Z}$.

Next consider the computation $t := \circ(\ell_1h_2 + s)$. Since $\ell_1 \in 2^{-110}\mathbb{Z}$ and $h_2 \in 2^{-52}\mathbb{Z}$, we have $\ell_1h_2 + s \in 2^{-162}\mathbb{Z}$; since $|\ell_1| \le 2^{-53}$ and $1 \le h_2 \le 2^{63/64} \cdot (1 + 2^{-53}) < 1.98$, we have also $|\ell_1h_2 + s| < 2^{-53} \cdot 1.98 + 2^{-52} < 2^{-51.007}$. Hence there is no underflow/overflow when computing $t := \circ(\ell_1h_2 + s)$, and we have $t \in 2^{-162}\mathbb{Z}$; furthermore, $t = \ell_1h_2 + s + \epsilon_2$ with $|\epsilon_2| \le \text{ulp}(\ell_1h_2 + s) \le 2^{-104}$, which implies that $|t| \le 2^{-51.007} + 2^{-104} < 2^{-51.00699}$.

Let us finally consider $p_\ell := \circ(h_1\ell_2 + t)$. Using $h_1 \in 2^{-52}\mathbb{Z}$, $\ell_2 \in 2^{-111}\mathbb{Z}$, and $t \in 2^{-162}\mathbb{Z}$ gives $h_1\ell_2 + t \in 2^{-163}\mathbb{Z}$. From $1 \le h_1 \le 2^{63/2^{12}} \cdot (1 + 2^{-53}) < 1.010719$ and $|\ell_2| \le 2^{-53}$, we have also $|h_1\ell_2 + t| \le 1.010719 \cdot 2^{-53} + 2^{-51.00699} < 2^{-50.6805}$. Hence there is no overflow/underflow when computing $p_\ell := \circ(h_1\ell_2 + t)$, and we have $p_\ell \in 2^{-163}\mathbb{Z}$; furthermore, $p_\ell = h_1\ell_2 + t + \epsilon_3$ with $|\epsilon_3| \le \text{ulp}(h_1\ell_2 + t) \le 2^{-103}$, which implies that $|p_\ell| \le 2^{-50.6805} + 2^{-103} \le 2^{-50.680499}$.

Let us now show that $p_h + p_\ell \in [1, 2)$. First, if $i_1 = i_2 = 0$, then $h_1 = h_2 = 1$, $\ell_1 = \ell_2 = 0$, $p_h = 1$, $s = t = p_\ell = 0$. Hence in this case $p_h + p_\ell = 1 \in [1, 2)$. Assume now that $i_1 \ge 1$ or $i_2 \ge 1$. Then either $h_1 \ge 2^{1/2^{12}} \cdot (1 - 2^{-53})$ or $h_2 \ge 2^{1/64} \cdot (1 - 2^{-53})$, and so $p_h := \circ(h_1h_2) \ge 2^{1/2^{12}} \cdot (1 - 2^{-53})(1 - 2^{-52}) > 1.0001$. Using the bound $|p_\ell| \le 2^{-50.680499}$ seen just before, we deduce that $p_h + p_\ell \ge p_h - |p_\ell| \ge 1.0001 - 2^{-50.680499} \ge 1.00009 > 1$. On the other hand, $p_h + p_\ell \le h_1h_2(1 + 2^{-52}) + |p_\ell| \le 2^{0.999757} \cdot (1 + 2^{-52}) + 2^{-50.680499} < 1.9997 < 2$.

Finally, let us bound the relative error of $p_h + p_\ell$ with respect to $(h_1+\ell_1)(h_2+\ell_2)$. We have $p_h + p_\ell = (h_1h_2+\epsilon_1)+(h_1\ell_2+t+\epsilon_3)$ with $t = \ell_1h_2 - \epsilon_1 + \epsilon_2$, and so $p_h + p_\ell - (h_1+\ell_1)(h_2+\ell_2) = \epsilon_2 + \epsilon_3 - \ell_1\ell_2$. Hence $|p_h + p_\ell - (h_1+\ell_1)(h_2+\ell_2)| \le |\epsilon_2| + |\epsilon_3| + |\ell_1| \cdot |\ell_2| \le 2^{-104} + 2^{-103} + 2^{-53} \cdot 2^{-53} < 2^{-102.299}$.

In order to deduce a relative error bound, let us provide a lower bound on the exact product $(h_1 + \ell_1)(h_2+\ell_2)$. Note first that if $i_1 = 0$, then $(h_1, \ell_1) =$

(1,0), which implies that $p_h := \circ(h_1h_2) = h_2$, $s := \circ(h_1h_2 - p_h) = 0$, $t := \circ(\ell_1h_2 + s) = 0$, and $p_\ell := \circ(h_1\ell_2 + t) = \ell_2$. Hence, in this case, $p_h + p_\ell = h_2 + \ell_2$ is the exact product. Similarly, one can check that if $i_2 = 0$, then $p_h + p_\ell = h_1 + \ell_1$ is also the exact product. Hence it remains to consider the case where $i_1 \ge 1$ and $i_2 \ge 1$. In this case, we have $h_1 \ge 2^{1/2^{12}} \cdot (1 - 2^{-53})$ and $h_2 \ge 2^{1/64} \cdot (1 - 2^{-53})$, which together with the bounds $|\ell_1|, |\ell_2| \le 2^{-53}$ leads to $(h_1 + \ell_1)(h_2 + \ell_2) > 1.0110603$. Hence, using $2^{-102.299}/1.0110603 < 2^{-102.314869}$, we arrive at the relative error bound

$$ \left| \frac{p_h + p_\ell}{(h_1 + \ell_1)(h_2 + \ell_2)} - 1 \right| \le 2^{-102.314869}. $$

By combining this bound with the relative error bounds $2^{-107.0228}$ and $2^{-107.57149}$ seen before for $h_1 + \ell_1$ and $h_2 + \ell_2$, we deduce that

$$ \left| \frac{p_h + p_\ell}{2^{i_2/64 + i_1/2^{12}}} - 1 \right| \le 2^{-102.2248}. $$

**About the values $h$ and $\ell$.** These values are defined by the four operations $h := \circ(p_hq_h)$, $s := \circ(p_hq_h - h)$, $t := \circ(p_\ell q_h + s)$, and $\ell := \circ(p_hq_\ell + t)$.

From Lemma 6, we know that $|q_\ell| \le 2^{-51.999}$ and that $q_h = \exp(z) \cdot (1 + \delta) - q_\ell$ with $|z| \le 2^{-12.905}$ and $|\delta| \le 2^{-64.90263}$. It follows that $q_h \in [0.99986, 1.000131]$. Hence, since $q_h \in \mathbb{F}$, this implies $q_h \in 2^{-53}\mathbb{Z}$. Together with $p_h = \circ(h_1h_2) \in \mathbb{F} \cap [1, 2)$, this also implies that there is no underflow/overflow when computing $h := \circ(p_hq_h)$ and that $p_hq_h - h$ is in $\mathbb{F}$. Furthermore, $h \ge 0.99986 \cdot (1 - 2^{-52}) \ge 0.999859$, so that $h \in 2^{-53}\mathbb{Z}$, and it then follows from $p_h \in 2^{-52}\mathbb{Z}$ and $q_h \in 2^{-53}\mathbb{Z}$ that $p_hq_h - h \in 2^{-105}\mathbb{Z}$. Hence $s = p_hq_h - h \in 2^{-105}\mathbb{Z}$ and, since $p_hq_h \le h_1h_2 \cdot (1 + 2^{-52}) \cdot q_h \le 2^{0.999757} \cdot (1 + 2^{-52}) \cdot 1.000131 < 2$, we have $h \le 2$ and $|s| \le \text{ulp}(p_hq_h) \le 2^{-52}$.

Using $p_\ell \in 2^{-163}\mathbb{Z}$, $q_h \in 2^{-53}\mathbb{Z}$, and $s \in 2^{-105}\mathbb{Z}$ gives $p_\ell q_h + s \in 2^{-216}\mathbb{Z}$. On the other hand, $|p_\ell q_h + s| \le 2^{-50.680499} \cdot 1.000131 + 2^{-52} \le 2^{-50.19424}$. Hence there is no underflow/overflow when computing $t := \circ(p_\ell q_h + s)$, and $t = p_\ell q_h + s + \epsilon_2$ with $|\epsilon_2| \le \text{ulp}(p_\ell q_h + s) \le 2^{-103}$. It follows that $|t| \le 2^{-50.19424} + 2^{-103} \le 2^{-50.194239}$.

Consider finally the operation $\ell := \circ(p_hq_\ell + t)$. We have $|p_hq_\ell + t| \le 2 \cdot 2^{-51.999} + 2^{-50.194239} \le 2^{-49.541218}$, so there is no overflow. Hence $\ell = p_hq_\ell + t + \epsilon_3$ with $|\epsilon_3| \le \text{ulp}(p_hq_\ell + t) \le 2^{-102}$ (which handles the possibility of underflow). It follows that $|\ell| \le 2^{-49.541218} + 2^{-102} \le 2^{-49.5412179}$, which together with $h \ge 0.999859$ shows that the ratio $\ell/h$ satisfies $|\ell/h| \le 2^{-49.541}$.

We also deduce that $h+\ell = (p_hq_h - s) + (p_hq_h + t + \epsilon_3)$ with $t = p_\ell q_h + s + \epsilon_2$, so that $h + \ell - (p_h + p_\ell)(q_h + q_\ell) = \epsilon_2 + \epsilon_3 - p_\ell q_\ell$. Consequently, $|h+\ell - (p_h + p_\ell)(q_h + q_\ell)| \le |\epsilon_2| + |\epsilon_3| + |p_\ell| \cdot |q_\ell| \le 2^{-103} + 2^{-102} + 2^{-50.680499} \cdot 2^{-51.999} \le 2^{-100.9129}$.
Note also that an exhaustive search over the 64 possible values gives a relative error always less than $2^{-107.57149}$. Furthermore, $h_2 \in \mathbb{F} \cap [1, 2)$ and $\ell_2 \in \mathbb{F}$ is either 0 or such that $|\ell_2| > 2^{-58.98} \ge 2^{-59}$, so that $h_2 \in 2^{-52}\mathbb{Z}$ and $\ell_2 \in 2^{-111}\mathbb{Z}$.

Finally, by definition, $1 \le h_1 \le 2^{63/2^{12}} \cdot (1 + 2^{-53}) < 2^{0.015381}$ and $1 \le h_2 \le 2^{63/64} \cdot (1 + 2^{-53}) < 2^{0.984376}$, which implies that the exact product $h_1h_2$ satisfies $1 \le h_1h_2 < 2^{0.999757} < 2$.

**About the values $p_h$ and $p_\ell$.** Since $h_1, h_2 \in \mathbb{F} \cap [1, 2)$, there is no underflow/overflow when computing $p_h := \circ(h_1h_2)$, and $h_1h_2 - p_h \in \mathbb{F}$. Hence, recalling that $1 \le h_1h_2 < 2$, we have $p_h = h_1h_2 + \epsilon_1$ with $|\epsilon_1| \le \text{ulp}(h_1h_2) \le 2^{-52}$. Since $h_1, h_2, p_h \in 2^{-52}\mathbb{Z}$, we have furthermore $|h_1h_2 - p_h| \in 2^{-104}\mathbb{Z}$, so there is no underflow/overflow when computing $s := \circ(h_1h_2 - p_h)$, and we thus have $s = h_1h_2 - p_h$ and $2^{-52} \ge |s| \in 2^{-104}\mathbb{Z}$.

Next consider the computation $t := \circ(\ell_1h_2 + s)$. Since $\ell_1 \in 2^{-110}\mathbb{Z}$ and $h_2 \in 2^{-52}\mathbb{Z}$, we have $\ell_1h_2 + s \in 2^{-162}\mathbb{Z}$; since $|\ell_1| \le 2^{-53}$ and $1 \le h_2 \le 2^{63/64} \cdot (1 + 2^{-53}) < 1.98$, we have also $|\ell_1h_2 + s| < 2^{-53} \cdot 1.98 + 2^{-52} < 2^{-51.007}$. Hence there is no underflow/overflow when computing $t := \circ(\ell_1h_2 + s)$, and we have $t \in 2^{-162}\mathbb{Z}$; furthermore, $t = \ell_1h_2 + s + \epsilon_2$ with $|\epsilon_2| \le \text{ulp}(\ell_1h_2 + s) \le 2^{-104}$, which implies that $|t| \le 2^{-51.007} + 2^{-104} < 2^{-51.00699}$.

Let us finally consider $p_\ell := \circ(h_1\ell_2 + t)$. Using $h_1 \in 2^{-52}\mathbb{Z}$, $\ell_2 \in 2^{-111}\mathbb{Z}$, and $t \in 2^{-162}\mathbb{Z}$ gives $h_1\ell_2 + t \in 2^{-163}\mathbb{Z}$. From $1 \le h_1 \le 2^{63/2^{12}} \cdot (1 + 2^{-53}) < 1.010719$ and $|\ell_2| \le 2^{-53}$, we have also $|h_1\ell_2 + t| \le 1.010719 \cdot 2^{-53} + 2^{-51.00699} < 2^{-50.6805}$. Hence there is no overflow/underflow when computing $p_\ell := \circ(h_1\ell_2 + t)$, and we have $p_\ell \in 2^{-163}\mathbb{Z}$; furthermore, $p_\ell = h_1\ell_2 + t + \epsilon_3$ with $|\epsilon_3| \le \text{ulp}(h_1\ell_2 + t) \le 2^{-103}$, which implies that $|p_\ell| \le 2^{-50.6805} + 2^{-103} \le 2^{-50.680499}$.

Let us now show that $p_h + p_\ell \in [1, 2)$. First, if $i_1 = i_2 = 0$, then $h_1 = h_2 = 1$, $\ell_1 = \ell_2 = 0$, $p_h = 1$, $s = t = p_\ell = 0$. Hence in this case $p_h + p_\ell = 1 \in [1, 2)$. Assume now that $i_1 \ge 1$ or $i_2 \ge 1$. Then either $h_1 \ge 2^{1/2^{12}} \cdot (1 - 2^{-53})$ or $h_2 \ge 2^{1/64} \cdot (1 - 2^{-53})$, and so $p_h := \circ(h_1h_2) \ge 2^{1/2^{12}} \cdot (1 - 2^{-53})(1 - 2^{-52}) > 1.0001$. Using the bound $|p_\ell| \le 2^{-50.680499}$ seen just before, we deduce that $p_h + p_\ell \ge p_h - |p_\ell| \ge 1.0001 - 2^{-50.680499} \ge 1.00009 > 1$. On the other hand, $p_h + p_\ell \le h_1h_2(1 + 2^{-52}) + |p_\ell| \le 2^{0.999757} \cdot (1 + 2^{-52}) + 2^{-50.680499} < 1.9997 < 2$.

Finally, let us bound the relative error of $p_h + p_\ell$ with respect to $(h_1+\ell_1)(h_2+\ell_2)$. We have $p_h + p_\ell = (h_1h_2+\epsilon_1)+(h_1\ell_2+t+\epsilon_3)$ with $t = \ell_1h_2 - \epsilon_1 + \epsilon_2$, and so $p_h + p_\ell - (h_1+\ell_1)(h_2+\ell_2) = \epsilon_2 + \epsilon_3 - \ell_1\ell_2$. Hence $|p_h + p_\ell - (h_1+\ell_1)(h_2+\ell_2)| \le |\epsilon_2| + |\epsilon_3| + |\ell_1| \cdot |\ell_2| \le 2^{-104} + 2^{-103} + 2^{-53} \cdot 2^{-53} < 2^{-102.299}$.

In order to deduce a relative error bound, let us provide a lower bound on the exact product $(h_1 + \ell_1)(h_2+\ell_2)$. Note first that if $i_1 = 0$, then $(h_1, \ell_1) =$

(1, 0), which implies that $p_h := \circ(h_1h_2) = h_2$, $s := \circ(h_1h_2 - p_h) = 0$, $t := \circ(\ell_1h_2 + s) = 0$, and $p_\ell := \circ(h_1\ell_2 + t) = \ell_2$. Hence, in this case, $p_h + p_\ell = h_2 + \ell_2$ is the exact product. Similarly, one can check that if $i_2 = 0$, then $p_h + p_\ell = h_1 + \ell_1$ is also the exact product. Hence it remains to consider the case where $i_1 \ge 1$ and $i_2 \ge 1$. In this case, we have $h_1 \ge 2^{1/2^{12}} \cdot (1 - 2^{-53})$ and $h_2 \ge 2^{1/64} \cdot (1 - 2^{-53})$, which together with the bounds $|\ell_1|, |\ell_2| \le 2^{-53}$ leads to $(h_1 + \ell_1)(h_2 + \ell_2) > 1.0110603$. Hence, using $2^{-102.299}/1.0110603 < 2^{-102.314869}$, we arrive at the relative error bound

$$ \left| \frac{p_h + p_\ell}{(h_1 + \ell_1)(h_2 + \ell_2)} - 1 \right| \le 2^{-102.314869}. $$

By combining this bound with the relative error bounds $2^{-107.0228}$ and $2^{-107.57149}$ seen before for $h_1 + \ell_1$ and $h_2 + \ell_2$, we deduce that

$$ \left| \frac{p_h + p_\ell}{2^{i_2/64 + i_1/2^{12}}} - 1 \right| \le 2^{-102.2248}. $$

**About the values $h$ and $\ell$.** These values are defined by the four operations $h := \circ(p_hq_h)$, $s := \circ(p_hq_h - h)$, $t := \circ(p_\ell q_h + s)$, and $\ell := \circ(p_hq_\ell + t)$.

From Lemma 6, we know that $|q_\ell| \le 2^{-51.999}$ and that $q_h = \exp(z) \cdot (1 + \delta) - q_\ell$ with $|z| \le 2^{-12.905}$ and $|\delta| \le 2^{-64.90263}$. It follows that $q_h \in [0.99986, 1.000131]$. Hence, since $q_h \in \mathbb{F}$, this implies $q_h \in 2^{-53}\mathbb{Z}$. Together with $p_h = \circ(h_1h_2) \in \mathbb{F} \cap [1, 2)$, this also implies that there is no underflow/overflow when computing $h := \circ(p_hq_h)$ and that $p_hq_h - h$ is in $\mathbb{F}$. Furthermore, $h \ge 0.99986 \cdot (1 - 2^{-52}) \ge 0.999859$, so that $h \in 2^{-53}\mathbb{Z}$, and it then follows from $p_h \in 2^{-52}\mathbb{Z}$ and $q_h \in 2^{-53}\mathbb{Z}$ that $p_hq_h - h \in 2^{-105}\mathbb{Z}$. Hence $s = p_hq_h - h \in 2^{-105}\mathbb{Z}$ and, since $p_hq_h \le h_1h_2 \cdot (1 + 2^{-52}) \cdot q_h \le 2^{0.999757} \cdot (1 + 2^{-52}) \cdot 1.000131 < 2$, we have $h \le 2$ and $|s| \le \text{ulp}(p_hq_h) \le 2^{-52}$.

Using $p_\ell \in 2^{-163}\mathbb{Z}$, $q_h \in 2^{-53}\mathbb{Z}$, and $s \in 2^{-105}\mathbb{Z}$ gives $p_\ell q_h + s \in 2^{-216}\mathbb{Z}$. On the other hand, $|p_\ell q_h + s| \le 2^{-50.680499} \cdot 1.000131 + 2^{-52} \le 2^{-50.19424}$. Hence there is no underflow/overflow when computing $t := \circ(p_\ell q_h + s)$, and $t = p_\ell q_h + s + \epsilon_2$ with $|\epsilon_2| \le \text{ulp}(p_\ell q_h + s) \le 2^{-103}$. It follows that $|t| \le 2^{-50.19424} + 2^{-103} \le 2^{-50.194239}$.

Consider finally the operation $\ell := \circ(p_hq_\ell + t)$. We have $|p_hq_\ell + t| \le 2 \cdot 2^{-51.999} + 2^{-50.194239} \le 2^{-49.541218}$, so there is no overflow. Hence $\ell = p_hq_\ell + t + \epsilon_3$ with $|\epsilon_3| \le \text{ulp}(p_hq_\ell + t) \le 2^{-102}$ (which handles the possibility of underflow). It follows that $|\ell| \le 2^{-49.541218} + 2^{-102} \le 2^{-49.5412179}$, which together with $h \ge 0.999859$ shows that the ratio $\ell/h$ satisfies $|\ell/h| \le 2^{-49.541}$.

We also deduce that $h+\ell = (p_hq_h - s) + (p_hq_h + t + \epsilon_3)$ with $t = p_\ell q_h + s + \epsilon_2$, so that $h + \ell - (p_h + p_\ell)(q_h + q_\ell) = \epsilon_2 + \epsilon_3 - p_\ell q_\ell$. Consequently, $|h+\ell - (p_h + p_\ell)(q_h + q_\ell)| \le |\epsilon_2| + |\epsilon_3| + |p_\ell| \cdot |q_\ell| \le 2^{-103} + 2^{-102} + 2^{-50.680499} \cdot 2^{-51.999} \le 2^{-100.9129}$.
In order to deduce a relative error bound, we use the following lower bound on the exact product: with $p_h \ge 1$, $q_h \ge 0.99986$, $|p_\ell| \le 2^{-50.680499}$, and $|q_\ell| \le 2^{-51.999}$, we have $(p_h + p_\ell)(q_h + q_\ell) \ge (p_h - |p_\ell|)(q_h - |q_\ell|) \ge 0.999859 \ge 2^{-0.000204}$. Therefore,

$$ \left| \frac{h + \ell}{(p_h + p_\ell)(q_h + q_\ell)} - 1 \right| \le 2^{-100.912696} $$

**About the values $e_h$ and $e_\ell$.** These values are defined by $e_h = \circ(2^e \cdot h)$ and $e_\ell = \circ(2^e \cdot \ell)$. Appendix A-J shows that when $\rho_1 \le r_h \le \rho_2$, we have $2^e \cdot h \in [2^{-991}, \Omega]$, thus the operation $e_h = \circ(2^e \cdot h)$ is exact. Since $|\ell/h| < 2^{-49.541}$ and $h \in [0.999859, 2]$, we have $|\ell| < \Omega$, thus no overflow can occur in $\circ(2^e \cdot \ell)$, but some underflow can occur. We thus have $e_\ell = 2^e \cdot \ell + \varepsilon$ with $|\varepsilon| \le \alpha$. Using $2^e \cdot h \ge 2^{-1022}$, it follows that

$$ \begin{aligned} |e_\ell/e_h| &\le |\ell/h| + |\varepsilon|/(2^e \cdot h) \\ &\le 2^{-49.541} + \alpha/2^{-1022} \le 2^{-49.2999}. \end{aligned} $$

**J. About $2^e \cdot h$ for $\rho_1 \le r_h \le \rho_2$**

We prove here that if $\rho_1 \le r_h \le \rho_2$, no underflow or overflow happens at the end of Algorithm exp_1. More precisely the approximation $2^e \cdot (h + \ell)$ of $\exp(r_h + r_\ell)$ obtained at line 25 of Algorithm exp_1 satisfies:

$$ 2^{-991} \le 2^e \cdot (h + \ell) \le \Omega, $$

moreover the same inequality holds for $2^e \cdot h$ alone, and in turn for $e_h = \circ(2^e \cdot h)$ since $2^{-991}$ and $\Omega$ are in $\mathbb{F}$.

In Appendices A-H and A-I we have analyzed Algorithms 7 and 8 (q_1 and exp_1) and shown that $|e_\ell/e_h| \le 2^{-49.2999}$ and that the following relative error bounds hold:

$$ \begin{aligned} \exp(r_h + r_\ell) &= 2^{k/2^{12}} \cdot \exp(z) \cdot (1 + \varepsilon), \\ p_h + p_\ell &= 2^{i_2/64 + i_1/2^{12}} \cdot (1 + \varepsilon'), \\ q_h + q_\ell &= \exp(z) \cdot (1 + \varepsilon''), \\ h + \ell &= (p_h + p_\ell)(q_h + q_\ell)(1 + \varepsilon''') \end{aligned} $$

where $|\varepsilon| \le 2^{-64.67806}$, $|\varepsilon'| \le 2^{-102.2248}$, $|\varepsilon''| \le 2^{-64.902632}$, and $|\varepsilon'''| \le 2^{-100.912696}$. Recalling that $k/2^{12} = e + i_2/64 + i_1/2^{12}$, we deduce that

$$ 2^e \cdot (h + \ell) = \exp(r_h + r_\ell) \cdot (1 + \delta), \quad (2) $$

where

$$ \delta := \frac{(1 + \varepsilon')(1 + \varepsilon'')(1 + \varepsilon''')}{1 + \varepsilon} - 1. $$

It can then be checked that $|\delta| \le \phi := 2^{-63.78598}$.

The smallest possible value of $r_h + r_\ell$ is obtained for $r_h = \rho_1$ and $r_\ell$ the smallest value compatible with $|r_\ell/r_h| < 2^{-23.8899}$, namely $r_\ell = -0x1.72b0feb06bbe9p-15$. For these values we get $\exp(r_h + r_\ell)(1 - \phi) > 2^{-991}$, which proves through Eq. (2) that $2^e \cdot (h + \ell) \ge 2^{-991}$. If

$\ell \le 0$, we have $2^e \cdot h \ge 2^e \cdot (h + \ell) \ge 2^{-991}$. If $\ell > 0$, we also know that $|\ell/h| \le 2^{-49.541}$, thus $2^e \cdot h \ge 2^e \cdot (h + \ell)(1 - 2^{-49.541})$, which yields

$$ 2^e \cdot h \ge \exp(r_h + r_\ell)(1 + \phi)(1 - 2^{-49.541}). $$

Again with the smallest possible value of $r_h + r_\ell$, we get a right-hand side larger than $2^{-991}$.

The largest possible value of $r_h + r_\ell$ is obtained for $r_h = \rho_2$ and $r_\ell$ the largest value compatible with $|r_\ell/r_h| < 2^{-23.8899}$, namely $r_\ell = 0x1.7f09093c9fe5bp-15$. For these values we get $\exp(r_h + r_\ell)(1 + \phi) < \Omega$, which proves that $2^e \cdot (h+\ell) \le \Omega$. If $\ell \ge 0$, we have $2^e \cdot h \le 2^e \cdot (h+\ell) \le \Omega$. If $\ell < 0$, we also know that $|\ell/h| \le 2^{-49.541}$, which yields

$$ 2^e \cdot h \le 2^e \cdot (h + \ell)(1 + 2^{-49.541}/(1 - 2^{-49.541})). $$

Again with the largest possible value of $r_h + r_\ell$, we get a right-hand side smaller than $\Omega$.

### K. Proof of Lemma 7

We have seen in Appendix A-J that under the hypothesis of Lemma 7, the values $2^e \cdot h$ and $2^e (h+\ell)$ at the end of Algorithm exp_1 do not underflow nor overflow, and that the operation $e_h = \circ(2^e \cdot h)$ is exact.

Still in Appendix A-J, we have proven:

$$ 2^e \cdot (h + \ell) = \exp(r_h + r_\ell) \cdot (1 + \delta), $$

with $|\delta| \le 2^{-63.78598}$.

However, we need a relative error bound in terms of $e_h + e_\ell$, not in terms of $2^e \cdot (h + \ell)$. While we know that $2^e \cdot h$ does not underflow nor overflow, thus $e_h = \circ(2^e \cdot h)$ is exact, and that $2^e \cdot \ell$ does not overflow, $2^e \cdot \ell$ might underflow, thus the operation $e_\ell = \circ(2^e \cdot \ell)$ at line 25 of Algorithm exp_1 might be inexact. It might give a maximal absolute error of $\alpha$, and since $2^e \cdot h \ge 2^{-991}$, $|\ell/h| < 2^{-49.541}$, this corresponds to an additional relative error:

$$ e_h + e_\ell = 2^e \cdot (h + \ell) \cdot (1 + \mu), $$

with $|\mu| < \alpha \cdot 2^{991}/(1 - 2^{-49.541}) \le 2^{-82.9}$. Thus

$$ e_h + e_\ell = \exp(r_h + r_\ell) \cdot (1 + \varepsilon_{\exp}), $$

with $\varepsilon_{\exp} := (1+\delta)(1+\mu)-1$. It can then be checked that $|\varepsilon_{\exp}| < 2^{-63.78597}$, and this concludes the proof of Lemma 7.

### L. Proof of Theorem 1 when $r_h < \rho_1$.

This can happen for example for $x = \Omega$ and $y = -\Omega$, where we get $y \log x \approx -2^{1033.47}$. The threshold $\rho_0$ is such that Eq. (3) yields $x^y < \alpha/2$ for $r_h < \rho_0$ (since $y \log x$ is negative, we need a lower bound of $|y \log x|$). In such a case, the correct rounding is +0, except for rounding up where it is $\alpha$. Note that in this case, $yh$ is in the normal range. Indeed, if $yh$ is in the subnormal range, i.e., $|yh| < 2^{-1022}$,
In order to deduce a relative error bound, we use the following lower bound on the exact product: with $p_h \ge 1$, $q_h \ge 0.99986$, $|p_\ell| \le 2^{-50.680499}$, and $|q_\ell| \le 2^{-51.999}$, we have $(p_h + p_\ell)(q_h + q_\ell) \ge (p_h - |p_\ell|)(q_h - |q_\ell|) \ge 0.999859 \ge 2^{-0.000204}$. Therefore,

$$ \left| \frac{h + \ell}{(p_h + p_\ell)(q_h + q_\ell)} - 1 \right| \le 2^{-100.912696}. $$

**About the values $e_h$ and $e_\ell$.** These values are defined by $e_h = \circ(2^e \cdot h)$ and $e_\ell = \circ(2^e \cdot \ell)$. Appendix A-J shows that when $\rho_1 \le r_h \le \rho_2$, we have $2^e \cdot h \in [2^{-991}, \Omega]$, thus the operation $e_h = \circ(2^e \cdot h)$ is exact. Since $|\ell/h| < 2^{-49.541}$ and $h \in [0.999859, 2]$, we have $|\ell| < \Omega$, thus no overflow can occur in $\circ(2^e \cdot \ell)$, but some underflow can occur. We thus have $e_\ell = 2^e \cdot \ell + \varepsilon$ with $|\varepsilon| \le \alpha$. Using $2^e \cdot h \ge 2^{-1022}$, it follows that

$$ \begin{aligned} |e_\ell/e_h| & \le |\ell/h| + |\varepsilon|(2^e \cdot h) \\ & \le 2^{-49.541} + \alpha/2^{-1022} \le 2^{-49.2999}. \end{aligned} $$

**J. About $2^e \cdot h$ for $\rho_1 \le r_h \le \rho_2$**.

We prove here that if $\rho_1 \le r_h \le \rho_2$, no underflow or overflow happens at the end of Algorithm exp_1. More precisely the approximation $2^e \cdot (h + \ell)$ of $\exp(r_h + r_\ell)$ obtained at line 25 of Algorithm exp_1 satisfies:

$$ 2^{-991} \le 2^e \cdot (h + \ell) \le \Omega, $$

moreover the same inequality holds for $2^e \cdot h$ alone, and in turn for $e_h = \circ(2^e \cdot h)$ since $2^{-991}$ and $\Omega$ are in $\mathbb{F}$.

In Appendices A-H and A-I we have analyzed Algorithms 7 and 8 (q_1 and exp_1) and shown that $|e_\ell/e_h| \le 2^{-49.2999}$ and that the following relative error bounds hold:

$$ \begin{aligned} \exp(r_h + r_\ell) &= 2^{k/2^{12}} \cdot \exp(z) \cdot (1 + \varepsilon), \\ p_h + p_\ell &= 2^{i_2/64 + i_1/2^{12}} \cdot (1 + \varepsilon'), \\ q_h + q_\ell &= \exp(z) \cdot (1 + \varepsilon''), \\ h + \ell &= (p_h + p_\ell)(q_h + q_\ell)(1 + \varepsilon'''). \end{aligned} $$

where $|\varepsilon| \le 2^{-64.67806}$, $|\varepsilon'| \le 2^{-102.2248}$, $|\varepsilon''| \le 2^{-64.902632}$, and $|\varepsilon'''| \le 2^{-100.912696}$. Recalling that $k/2^{12} = e + i_2/64 + i_1/2^{12}$, we deduce that

$$ 2^e \cdot (h + \ell) = \exp(r_h + r_\ell) \cdot (1 + \delta), \quad (2) $$

where

$$ \delta := \frac{(1 + \varepsilon')(1 + \varepsilon'')(1 + \varepsilon''')}{1 + \varepsilon} - 1. $$

It can then be checked that $|\delta| \le \phi := 2^{-63.78598}$.

The smallest possible value of $r_h + r_\ell$ is obtained for $r_h = \rho_1$ and $r_\ell$ the smallest value compatible with $|r_\ell/r_h| < 2^{-23.8899}$, namely $r_\ell = -0x1.72b0feb06bbe9p-15$. For these values we get $\exp(r_h + r_\ell)(1 - \phi) > 2^{-991}$, which proves through Eq. (2) that $2^e \cdot (h + \ell) \ge 2^{-991}$. If

$\ell \le 0$, we have $2^e \cdot h \ge 2^e \cdot (h + \ell) \ge 2^{-991}$. If $\ell > 0$, we also know that $|\ell/h| \le 2^{-49.541}$, thus $2^e \cdot h \ge 2^e \cdot (h + \ell)(1 - 2^{-49.541})$, which yields

$$ 2^e \cdot h \ge \exp(r_h + r_\ell)(1 + \phi)(1 - 2^{-49.541}). $$

Again with the smallest possible value of $r_h + r_\ell$, we get a right-hand side larger than $2^{-991}$.

The largest possible value of $r_h + r_\ell$ is obtained for $r_h = \rho_2$ and $r_\ell$ the largest value compatible with $|r_\ell/r_h| < 2^{-23.8899}$, namely $r_\ell = 0x1.7f09093c9fe5bp-15$. For these values we get $\exp(r_h + r_\ell)(1 + \phi) < \Omega$, which proves that $2^e \cdot (h+\ell) \le \Omega$. If $\ell \ge 0$, we have $2^e \cdot h \le 2^e \cdot (h+\ell) \le \Omega$. If $\ell < 0$, we also know that $|\ell/h| \le 2^{-49.541}$, which yields

$$ 2^e \cdot h \le 2^e \cdot (h + \ell)(1 + 2^{-49.541}/(1 - 2^{-49.541})). $$

Again with the largest possible value of $r_h + r_\ell$, we get a right-hand side smaller than $\Omega$.

### K. Proof of Lemma 7

We have seen in Appendix A-J that under the hypothesis of Lemma 7, the values $2^e \cdot h$ and $2^e (h+\ell)$ at the end of Algorithm exp_1 do not underflow nor overflow, and that the operation $e_h = \circ(2^e \cdot h)$ is exact.

Still in Appendix A-J, we have proven:

$$ 2^e \cdot (h + \ell) = \exp(r_h + r_\ell) \cdot (1 + \delta), $$

with $|\delta| \le 2^{-63.78598}$.

However, we need a relative error bound in terms of $e_h + e_\ell$, not in terms of $2^e \cdot (h + \ell)$. While we know that $2^e \cdot h$ does not underflow nor overflow, thus $e_h = \circ(2^e \cdot h)$ is exact, and that $2^e \cdot \ell$ does not overflow, $2^e \cdot \ell$ might underflow, thus the operation $e_\ell = \circ(2^e \cdot \ell)$ at line 25 of Algorithm exp_1 might be inexact. It might give a maximal absolute error of $\alpha$, and since $2^e \cdot h \ge 2^{-991}$, $|\ell/h| < 2^{-49.541}$, this corresponds to an additional relative error:

$$ e_h + e_\ell = 2^e \cdot (h + \ell) \cdot (1 + \mu), $$

with $|\mu| < \alpha \cdot 2^{991}/(1 - 2^{-49.541}) \le 2^{-82.9}$. Thus

$$ e_h + e_\ell = \exp(r_h + r_\ell) \cdot (1 + \varepsilon_{\exp}), $$

with $\varepsilon_{\exp} := (1+\delta)(1+\mu)-1$. It can then be checked that $|\varepsilon_{\exp}| < 2^{-63.78597}$, and this concludes the proof of Lemma 7.

### L. Proof of Theorem 1 when $r_h < \rho_1$.

This can happen for example for $x = \Omega$ and $y = -\Omega$, where we get $y \log x \approx -2^{1033.47}$. The threshold $\rho_0$ is such that Eq. (3) yields $x^y < \alpha/2$ for $r_h < \rho_0$ (since $y \log x$ is negative, we need a lower bound of $|y \log x|$). In such a case, the correct rounding is +0, except for rounding up where it is $\alpha$. Note that in this case, $yh$ is in the normal range. Indeed, if $yh$ is in the subnormal range, i.e., $|yh| < 2^{-1022}$,
since $r_h = \circ(yh)$, we have $|r_h| \le 2^{-1022} + \alpha$, which is incompatible with $r_h < \rho_1 \approx -686.909$.

For $r_h < \rho_0$, Algorithm exp_1 returns $e_h = \alpha$ and $e_\ell = -\alpha$. For rounding to nearest, this yields $\circ(e_\ell \pm \varepsilon e_h) = -\alpha$, thus $u = v = +0$ (remember that $(+x) + (-x)$ is $+0$ for rounding to nearest). For rounding toward zero or upwards, this yields $\circ(e_\ell - \varepsilon e_h) = -\alpha$ and $\circ(e_\ell + \varepsilon e_h) = -0$, thus $u = +0$ and $v = \alpha$, and the rounding test fails. For rounding down, this yields $\circ(e_\ell - \varepsilon e_h) = -2\alpha$ and $\circ(e_\ell + \varepsilon e_h) = -\alpha$, thus $u = -\alpha$ and $v = -0$, and the rounding test fails. Thus for rounding to nearest, the correct rounding is returned, and for directed rounding modes, the correct rounding is delegated to the second phase.

The threshold $\rho_1$ is computed such that for $r_h \ge \rho_1$, the value $2^e \cdot h$ at the end of Algorithm exp_1 is larger or equal to $2^{-991}$, thus $e_h = \circ(2^e \cdot h)$ does not underflow (see Appendix A-J). In the small range $\rho_0 \le r_h < \rho_1$, Algorithm exp_1 returns $e_h = e_\ell = \text{NaN}$, and like in the case $\rho_2 < r_h \le \rho_3$, these corner cases are deferred to the second phase.

### M. Proof of Theorem 1 when $\rho_2 < r_h$.

This can happen for example for $x = \Omega$ and $y = \Omega$, where we get $y \log x \approx 2^{1033.47}$. Note that in this case, $yh$ is in the normal range. Indeed, if $yh$ is in the subnormal range, i.e., $|yh| < 2^{-1022}$, since $r_h = \circ(yh)$, we have $|r_h| \le 2^{-1022} + \alpha$, which is incompatible with $r_h > \rho_2 \approx 709.78267$.

From Eq. (1), we have

$$|\log x| \ge |h + \ell|/(1 + \varepsilon_{\log}),$$

and since $|\ell| \le 2^{-23.89}|h|$ from Lemma 4,

$$|y \log x| \ge |yh|(1 - 2^{-23.89})/(1 + \varepsilon_{\log}).$$

Since $r_h = \circ(yh)$, we have $|yh| \ge |r_h|(1 - 2^{-52})$, thus:

$$|y \log x| \ge |r_h| \frac{(1 - 2^{-23.89})(1 - 2^{-52})}{1 + \varepsilon_{\log}}. \quad (3)$$

For $r_h > \rho_3$, we have $r_h \ge \rho_3 + 2^{-43}$, and it then follows from Eq. (3) that $x^y > 2^{1024}$, thus $x^y$ overflows for all rounding modes.

In that case Algorithm exp_1 returns $e_h = e_\ell = \Omega$. In Algorithm phase_1, for rounding down or toward zero, this will yield $u = v = \Omega$, and since $u = v$ the algorithm returns $u = \Omega$, which is the correct rounding. For rounding up or to nearest, we obtain $u = v = +\infty$, and again since $u = v$ the algorithm returns $u = +\infty$ which is the correct rounding.

On the other hand, we also have:

$$|y \log x| \le |r_h| \frac{(1 + 2^{-23.89})(1 + 2^{-52})}{1 - \varepsilon_{\log}}, \quad (4)$$

and we can check that for $r_h \le \rho_2$, $x^y < \Omega$, thus no overflow happens.

The threshold $\rho_2$ is computed such that for $r_h \le \rho_2$, the value $2^e \cdot h$ at the end of Algorithm exp_1

is smaller than or equal to $\Omega$, thus $e_h = \circ(2^e \cdot h)$ does not overflow (see Appendix A-J). In the small region $\rho_2 < r_h \le \rho_3$, Algorithm exp_1 returns $e_h = e_\ell = \text{NaN}$, thus in Algorithm phase_1 we get $u = v = \text{NaN}$, and by IEEE 754, the comparison $u = v$ is false, FAIL is returned, and these corner cases are deferred to the second phase.

### N. Proof of Theorem 1 when $\rho_1 \le r_h \le \rho_2$.

Since we will use here Lemma 5, we have to check that the assumption $2^{-969} \le |yh| \le 709.7827$ from that lemma holds, or exclude the cases $|yh| < 2^{-969}$ and $709.7827 < |yh|$. The latter case is easy to exclude: if $|yh| > 709.7827$, since $r_h = \circ(yh)$, we deduce that $|r_h| \ge 709.7827 \cdot (1 - 2^{-52}) \ge 709.78269$, which is incompatible with $\rho_1 \le r_h \le \rho_2$. The case $|yh| < 2^{-969}$ is dealt with separately at the end of this subsection, where we show that Theorem 1 holds in this case too.

Step 1 of Algorithm phase_1 computes $h, \ell \leftarrow \log_1(x)$. According to Lemma 4, we have $|\ell| \le 2^{-23.89}|h|$ and

$$\left| \frac{h + \ell}{\log x} - 1 \right| \le \varepsilon_{\log}$$

with $\varepsilon_{\log} = 2^{-73.527}$ if $x \notin (1/\sqrt{2}, \sqrt{2})$, and $\varepsilon_{\log} = 2^{-67.0544}$ otherwise.

Step 2 of Algorithm phase_1 computes $r_h, r_\ell \leftarrow \text{mul}_1(h, \ell, y)$. According to Lemma 5, since we assumed $2^{-969} \le |yh| \le 709.7827$, then $|r_h| \in [2^{-970}, 709.79]$, $|r_\ell| \le 2^{-14.4187}$, $|r_\ell/r_h| \le 2^{-23.8899}$, $|r_h + r_\ell| \le 709.79$, and

$$|r_h + r_\ell - y \log x| \le \varepsilon_{\text{mul}} \quad (5)$$

with $\varepsilon_{\text{mul}} = 2^{-63.799}$ if $x \notin (1/\sqrt{2}, \sqrt{2})$, and $\varepsilon_{\text{mul}} = 2^{-57.580}$ otherwise.

Step 3 of Algorithm phase_1 computes $e_h, e_\ell \leftarrow \text{exp}_1(r_h, r_\ell)$. The assumptions $|r_\ell/r_h| < 2^{-23.8899}$ and $|r_\ell| < 2^{-14.4187}$ of Lemma 7 are fulfilled (by Lemma 5), underflow and overflow cases are detected within Algorithm exp_1 (see Appendices A-L and A-M), thus Lemma 7 applies, and we have

$$\left| \frac{e_h + e_\ell}{\exp(r_h + r_\ell)} - 1 \right| < 2^{-63.78597}.$$

Let us define $r := r_h + r_\ell$. Thus

$$|e_h + e_\ell - e^r| < 2^{-63.78597} \cdot e^r. \quad (6)$$

By Eq. (5), we can write $r = y \log x + \sigma$ with $|\sigma| \le \varepsilon_{\text{mul}}$. Taking the exponential yields

$$e^r = x^y \cdot e^\sigma,$$

which gives

$$|e^r - x^y| \le |e^{-\sigma} - 1| \cdot e^r. \quad (7)$$

Combining Eq. (6) and Eq. (7) yields

$$|e_h + e_\ell - x^y| \le (2^{-63.78597} + |e^{-\sigma} - 1|) \cdot e^r.$$
since $r_h = \circ(yh)$, we have $|r_h| \le 2^{-1022} + \alpha$, which is incompatible with $r_h < \rho_1 \approx -686.909$.

For $r_h < \rho_0$, Algorithm exp_1 returns $e_h = \alpha$ and $e_\ell = -\alpha$. For rounding to nearest, this yields $\circ(e_\ell \pm \varepsilon e_h) = -\alpha$, thus $u = v = +0$ (remember that $(+x) + (-x)$ is $+0$ for rounding to nearest). For rounding toward zero or upwards, this yields $\circ(e_\ell - \varepsilon e_h) = -\alpha$ and $\circ(e_\ell + \varepsilon e_h) = -0$, thus $u = +0$ and $v = \alpha$, and the rounding test fails. For rounding down, this yields $\circ(e_\ell - \varepsilon e_h) = -2\alpha$ and $\circ(e_\ell + \varepsilon e_h) = -\alpha$, thus $u = -\alpha$ and $v = -0$, and the rounding test fails. Thus for rounding to nearest, the correct rounding is returned, and for directed rounding modes, the correct rounding is delegated to the second phase.

The threshold $\rho_1$ is computed such that for $r_h \ge \rho_1$, the value $2^e \cdot h$ at the end of Algorithm exp_1 is larger or equal to $2^{-991}$, thus $e_h = \circ(2^e \cdot h)$ does not underflow (see Appendix A-J). In the small range $\rho_0 \le r_h < \rho_1$, Algorithm exp_1 returns $e_h = e_\ell = \text{NaN}$, and like in the case $\rho_2 < r_h \le \rho_3$, these corner cases are deferred to the second phase.

### M. Proof of Theorem 1 when $\rho_2 < r_h$.

This can happen for example for $x = \Omega$ and $y = \Omega$, where we get $y \log x \approx 2^{1033.47}$. Note that in this case, $yh$ is in the normal range. Indeed, if $yh$ is in the subnormal range, i.e., $|yh| < 2^{-1022}$, since $r_h = \circ(yh)$, we have $|r_h| \le 2^{-1022} + \alpha$, which is incompatible with $r_h > \rho_2 \approx 709.78267$.

From Eq. (1), we have

$$|\log x| \ge |h + \ell|/(1 + \varepsilon_{\log}),$$

and since $|\ell| \le 2^{-23.89}|h|$ from Lemma 4,

$$|y \log x| \ge |yh|(1 - 2^{-23.89})/(1 + \varepsilon_{\log}).$$

Since $r_h = \circ(yh)$, we have $|yh| \ge |r_h|(1 - 2^{-52})$, thus:

$$|y \log x| \ge |r_h| \frac{(1 - 2^{-23.89})(1 - 2^{-52})}{1 + \varepsilon_{\log}}. \quad (3)$$

For $r_h > \rho_3$, we have $r_h \ge \rho_3 + 2^{-43}$, and it then follows from Eq. (3) that $x^y > 2^{1024}$, thus $x^y$ overflows for all rounding modes.

In that case Algorithm exp_1 returns $e_h = e_\ell = \Omega$. In Algorithm phase_1, for rounding down or toward zero, this will yield $u = v = \Omega$, and since $u = v$ the algorithm returns $u = \Omega$, which is the correct rounding. For rounding up or to nearest, we obtain $u = v = +\infty$, and again since $u = v$ the algorithm returns $u = +\infty$ which is the correct rounding.

On the other hand, we also have:

$$|y \log x| \le |r_h| \frac{(1 + 2^{-23.89})(1 + 2^{-52})}{1 - \varepsilon_{\log}}, \quad (4)$$

and we can check that for $r_h \le \rho_2$, $x^y < \Omega$, thus no overflow happens.

The threshold $\rho_2$ is computed such that for $r_h \le \rho_2$, the value $2^e \cdot h$ at the end of Algorithm exp_1

is smaller than or equal to $\Omega$, thus $e_h = \circ(2^e \cdot h)$ does not overflow (see Appendix A-J). In the small region $\rho_2 < r_h \le \rho_3$, Algorithm exp_1 returns $e_h = e_\ell = \text{NaN}$, thus in Algorithm phase_1 we get $u = v = \text{NaN}$, and by IEEE 754, the comparison $u = v$ is false, FAIL is returned, and these corner cases are deferred to the second phase.

### N. Proof of Theorem 1 when $\rho_1 \le r_h \le \rho_2$.

Since we will use here Lemma 5, we have to check that the assumption $2^{-969} \le |yh| \le 709.7827$ from that lemma holds, or exclude the cases $|yh| < 2^{-969}$ and $709.7827 < |yh|$. The latter case is easy to exclude: if $|yh| > 709.7827$, since $r_h = \circ(yh)$, we deduce that $|r_h| \ge 709.7827 \cdot (1 - 2^{-52}) \ge 709.78269$, which is incompatible with $\rho_1 \le r_h \le \rho_2$. The case $|yh| < 2^{-969}$ is dealt with separately at the end of this subsection, where we show that Theorem 1 holds in this case too.

Step 1 of Algorithm phase_1 computes $h, \ell \leftarrow \log_1(x)$. According to Lemma 4, we have $|\ell| \le 2^{-23.89}|h|$ and

$$\left| \frac{h + \ell}{\log x} - 1 \right| \le \varepsilon_{\log}$$

with $\varepsilon_{\log} = 2^{-73.527}$ if $x \notin (1/\sqrt{2}, \sqrt{2})$, and $\varepsilon_{\log} = 2^{-67.0544}$ otherwise.

Step 2 of Algorithm phase_1 computes $r_h, r_\ell \leftarrow \text{mul}_1(h, \ell, y)$. According to Lemma 5, since we assumed $2^{-969} \le |yh| \le 709.7827$, then $|r_h| \in [2^{-970}, 709.79]$, $|r_\ell| \le 2^{-14.4187}$, $|r_\ell/r_h| \le 2^{-23.8899}$, $|r_h + r_\ell| \le 709.79$, and

$$|r_h + r_\ell - y \log x| \le \varepsilon_{\text{mul}} \quad (5)$$

with $\varepsilon_{\text{mul}} = 2^{-63.799}$ if $x \notin (1/\sqrt{2}, \sqrt{2})$, and $\varepsilon_{\text{mul}} = 2^{-57.580}$ otherwise.

Step 3 of Algorithm phase_1 computes $e_h, e_\ell \leftarrow \text{exp}_1(r_h, r_\ell)$. The assumptions $|r_\ell/r_h| < 2^{-23.8899}$ and $|r_\ell| < 2^{-14.4187}$ of Lemma 7 are fulfilled (by Lemma 5), underflow and overflow cases are detected within Algorithm exp_1 (see Appendices A-L and A-M), thus Lemma 7 applies, and we have

$$\left| \frac{e_h + e_\ell}{\exp(r_h + r_\ell)} - 1 \right| < 2^{-63.78597}.$$

Let us define $r := r_h + r_\ell$. Thus

$$|e_h + e_\ell - e^r| < 2^{-63.78597} \cdot e^r. \quad (6)$$

By Eq. (5), we can write $r = y \log x + \sigma$ with $|\sigma| \le \varepsilon_{\text{mul}}$. Taking the exponential yields

$$e^r = x^y \cdot e^\sigma,$$

which gives

$$|e^r - x^y| \le |e^{-\sigma} - 1| \cdot e^r. \quad (7)$$

Combining Eq. (6) and Eq. (7) yields

$$|e_h + e_\ell - x^y| \le (2^{-63.78597} + |e^{-\sigma} - 1|) \cdot e^r.$$
The largest value of $|e^{-\sigma} - 1|$ is obtained for $\sigma = -\epsilon_{\text{mul}}$, and we obtain

$$|e_h + e_\ell - x^y| \le \mu \cdot e^r \quad (8)$$

with $\mu = 2^{-62.7924}$ if $x \notin (1/\sqrt{2}, \sqrt{2})$, and $\mu = 2^{-57.5605}$ otherwise. Now we want a bound in terms of $e_h$, $e_\ell$, which are known, not in terms of $e^r$, which is unknown. Eq. (6) rewrites as $e^r < |e_h + e_\ell| + 2^{-63.78597}e^r$, which yields $e^r < 1/(1 - 2^{-63.78597}) \cdot |e_h + e_\ell|$. Moreover we want a bound in terms of $e_h$ alone, since the left and right bounds $u, v$ computed in steps 6-7 of Algorithm phase_1 use an error bound in terms of $\epsilon e_h$. Using the hypothesis $|e_\ell/e_h| \le 2^{-49.2999}$, it follows that

$$e^r < (1 + 2^{-49.2999}) / (1 - 2^{-63.78597}) \cdot |e_h|.$$

Combining this inequality with Eq. (8), we obtain

$$|e_h + e_\ell - x^y| \le \tau \cdot e_h, \quad (9)$$

with $\tau = 2^{-62.7923}$ if $x \notin (1/\sqrt{2}, \sqrt{2})$, and $\tau = 2^{-57.5604}$ otherwise.

**About the case** $|yh| < 2^{-969}$. In this case we analyze separately Algorithms mul_1 and exp_1.

First in mul_1, since $|yh| < 2^{-969}$ and $r_h = \circ(yh)$, we have $|r_h| \le 2^{-969}$ since $2^{-969} \in \mathbb{F}$. Since $r_h$ is the rounding of $yh$, by definition $|r_h - yh| < \text{ulp}(yh) \le 2^{-1022}$, thus $|s| \le 2^{-1022}$. Now since $|\ell| \le 2^{-23.89}|h|$, we have $|y\ell| < 2^{-992.89}$, thus $|y\ell + s| < 2^{-992.89} + 2^{-1022} \le 2^{-992.88}$, and $|r_\ell| \le 2^{-992.88} \cdot (1+2^{-52}) \le 2^{-992}$. In summary we have in this case $|r_h| \le 2^{-969}$ and $|r_\ell| \le 2^{-992}$.

Now if we analyze Algorithm exp_1 with such tiny values of $r_h$, $r_\ell$, we see that $k=0$, which yields $z_h = r_h$, $z_\ell = r_\ell$, $z = \circ(r_h + r_\ell)$ and thus $|z| \le (2^{-969} + 2^{-992})(1+2^{-52}) \le 2^{-968.9}$. In addition, $k=0$ implies $e = i_2 = i_1 = 0$, $h_2 = h_1 = 1$, $\ell_2 = \ell_1 = 0$, $p_h = 1$, $s = 0$, $t = p_\ell = 0$, and at the end of the algorithm, $h = q_h$, $s = 0$, $t = 0$, $\ell = q_\ell$, thus $e_h = q_h$ and $e_\ell = q_\ell$.

It thus remains to analyze Algorithm q_1 for tiny inputs, namely, $|z| \le 2^{-968.9}$. First, $q = \circ(Q_4z+Q_3)$ is rounded either to $Q_3 \approx 0.17$, or to one of the two adjacent binary64 numbers, thus $|q| \le 2^{-2}$. Similarly $q' = \circ(qz+Q_2)$ is rounded either to $Q_2 = 1/2$, or to one of the two adjacent binary64 numbers, thus $|q'| \le 1$. In $h_0 = \circ(q'z+Q_1)$, $q'z$ is so tiny that $h_0$ is either $Q_1 = 1$ or one of the two adjacent binary64 numbers, thus $|h_0| \le 2$. Then in $h_1$, $\ell_1 = \text{ExactMul}(z, h_0)$, we have $|h_1|, |\ell_1| \le 2^{-967.9}$. We thus perform the last operation $q_h$, $q_\ell = \text{FastSum}(Q_0, h_1, \ell_1)$, with $Q_0 = 1$ and $|h_1|, |\ell_1| \le 2^{-967.9}$. For rounding to nearest-even, it necessarily gives $q_h = 1$ and $|q_\ell| \le 2^{-965}$. Since we have $e_h = q_h = 1$ and $e_\ell = q_\ell$ at the end of exp_1, it is easy to check that the rounding test from Algorithm phase_1 will succeed, since $\circ(e_\ell \pm \epsilon e_h)$ rounds to $\pm \epsilon$, and $\circ(e_h \pm \epsilon)$ rounds to 1, which is the expected correct rounding for $\exp(y \log x)$ with tiny $y \log x$. For directed roundings, we will either

have $e_h = 1$ and tiny $e_\ell$, or $e_h = 1 - 2^{-53}$ and $e_\ell \approx 2^{-53}$, or $e_h = 1 + 2^{-52}$ and $e_\ell \approx -2^{-52}$. In the first case ($e_h = 1$ and tiny $e_\ell$), the rounding test will fail, since $\circ(e_\ell \pm \epsilon e_h)$ rounds to about $\pm \epsilon$, and $1-\epsilon, 1+\epsilon$ round to different values. In the second case ($e_h = 1 - 2^{-53}$ and $e_\ell \approx 2^{-53}$), $\circ(e_\ell \pm \epsilon e_h)$ rounds to about $2^{-53} \pm \epsilon$, thus we have approximately $u = \circ(1-\epsilon)$ and $v = \circ(1+\epsilon)$, which again round to different values. The third case is similar.

Note that this reasoning also holds when $z=0$, which can happen when $x=1$, for example. In this case we have $q_h = 1$, $q_\ell = 0$, thus $e_h = 1$, $e_\ell = 0$, $\circ(e_\ell \pm \epsilon e_h) = \pm \epsilon$, and the rounding test simplifies to $\circ(1-\epsilon) = \circ(1+\epsilon)$, which only holds for rounding to nearest.

In summary, when $|yh| < 2^{-969}$, Algorithm phase_1 returns the correct rounding for rounding to nearest-even, and returns FAIL for directed rounding modes, thus Theorem 1 holds.

**About the rounding test.** We know from (9) that

$$e_h + e_\ell - \tau e_h \le x^y \le e_h + e_\ell + \tau e_h,$$

from which we deduce by monotonicity of rounding:

$$\circ(e_h + e_\ell - \tau e_h) \le \circ(x^y) \le \circ(e_h + e_\ell + \tau e_h).$$

But since there is no instruction that directly computes $\circ(e_h + e_\ell \pm \tau e_h)$, in Algorithm phase_1, we instead compute $u := \circ(e_h + \circ(e_\ell - \epsilon e_h))$ and, similarly, $v := \circ(e_h + \circ(e_\ell + \epsilon e_h))$. We thus have an additional rounding error in $\circ(e_\ell \pm \epsilon e_h)$. Here, $\epsilon$ will be chosen slightly larger than $\tau$ (and thus a fortiori in $[\tau, 2\tau]$, say) in order to ensure that $u \le \circ(x^y) \le v$ despite these two extra rounding errors.

Let us now determine $\epsilon$. It suffices to have $\circ(e_\ell - \epsilon e_h) \le e_\ell - \tau e_h$ and $e_\ell + \tau e_h \le \circ(e_\ell + \epsilon e_h)$, since by monotonicity of rounding, this will ensure that $u \le \circ(x^y) \le v$. Now, defining for simplicity $u' := e_\ell - \epsilon e_h$ and $v' := e_\ell + \epsilon e_h$, we have $\circ(u') \le u' + \text{ulp}(u')$ and $v' - \text{ulp}(v') \le \circ(v')$, and so it suffices that $\epsilon$ be such that $-\epsilon e_h + \text{ulp}(u') \le -\tau e_h$ and $\tau e_h \le \epsilon e_h - \text{ulp}(v')$. Now, $|u'|, |v'| \le |e_\ell| + \epsilon e_h \le (2^{-49.2999} + 2\tau)e_h < 2^{-49.2905}e_h$ for $\epsilon \le 2\tau \le 2 \cdot 2^{-57.5604}$. If $|u'| \ge 2^{-1022}$ then $\text{ulp}(u') = 2^{-52}|u'| \le 2^{-101.2905}e_h$; if $|u'| < 2^{-1022}$, then $\text{ulp}(u') = 2^{-1074} \le 2^{-83}e_h$, since $e_h \ge 2^{-991}$. Consequently, even if underflow occurs, we have $\text{ulp}(u') \le 2^{-83}e_h$ and, similarly, $\text{ulp}(v') \le 2^{-83}e_h$. Hence it suffices that $\epsilon$ be such that $\epsilon \ge \tau + 2^{-83}$. It can then be checked that the values RU$(2^{-62.792})$ and RU$(2^{-57.560})$ taken for $\epsilon$ in Algorithm phase_1 satisfy this sufficient condition as well as the constraint $\epsilon \le 2\tau$.
The largest value of $|e^{-\sigma} - 1|$ is obtained for $\sigma = -\epsilon_{\text{mul}}$, and we obtain

$$|e_h + e_\ell - x^y| \le \mu \cdot e^r \quad (8)$$

with $\mu = 2^{-62.7924}$ if $x \notin (1/\sqrt{2}, \sqrt{2})$, and $\mu = 2^{-57.5605}$ otherwise. Now we want a bound in terms of $e_h$, $e_\ell$, which are known, not in terms of $e^r$, which is unknown. Eq. (6) rewrites as $e^r < |e_h + e_\ell| + 2^{-63.78597}e^r$, which yields $e^r < 1/(1 - 2^{-63.78597}) \cdot |e_h + e_\ell|$. Moreover we want a bound in terms of $e_h$ alone, since the left and right bounds $u, v$ computed in steps 6-7 of Algorithm phase_1 use an error bound in terms of $\epsilon e_h$. Using the hypothesis $|e_\ell/e_h| \le 2^{-49.2999}$, it follows that

$$e^r < (1 + 2^{-49.2999})/(1 - 2^{-63.78597}) \cdot |e_h|.$$

Combining this inequality with Eq. (8), we obtain

$$|e_h + e_\ell - x^y| \le \tau \cdot e_h, \quad (9)$$

with $\tau = 2^{-62.7923}$ if $x \notin (1/\sqrt{2}, \sqrt{2})$, and $\tau = 2^{-57.5604}$ otherwise.

**About the case** $|yh| < 2^{-969}$. In this case we analyze separately Algorithms mul_1 and exp_1.

First in mul_1, since $|yh| < 2^{-969}$ and $r_h = \circ(yh)$, we have $|r_h| \le 2^{-969}$ since $2^{-969} \in \mathbb{F}$. Since $r_h$ is the rounding of $yh$, by definition $|r_h - yh| < \text{ulp}(yh) \le 2^{-1022}$, thus $|s| \le 2^{-1022}$. Now since $|\ell| \le 2^{-23.89}|h|$, we have $|y\ell| < 2^{-992.89}$, thus $|y\ell + s| < 2^{-992.89} + 2^{-1022} \le 2^{-992.88}$, and $|r_\ell| \le 2^{-992.88} \cdot (1+2^{-52}) \le 2^{-992}$. In summary we have in this case $|r_h| \le 2^{-969}$ and $|r_\ell| \le 2^{-992}$.

Now if we analyze Algorithm exp_1 with such tiny values of $r_h$, $r_\ell$, we see that $k=0$, which yields $z_h = r_h$, $z_\ell = r_\ell$, $z = \circ(r_h + r_\ell)$ and thus $|z| \le (2^{-969} + 2^{-992})(1 + 2^{-52}) \le 2^{-968.9}$. In addition, $k=0$ implies $e = i_2 = i_1 = 0$, $h_2 = h_1 = 1$, $\ell_2 = \ell_1 = 0$, $p_h = 1$, $s = 0$, $t = p_\ell = 0$, and at the end of the algorithm, $h = q_h$, $s = 0$, $t = 0$, $\ell = q_\ell$, thus $e_h = q_h$ and $e_\ell = q_\ell$.

It thus remains to analyze Algorithm q_1 for tiny inputs, namely, $|z| \le 2^{-968.9}$. First, $q = \circ(Q_4z+Q_3)$ is rounded either to $Q_3 \approx 0.17$, or to one of the two adjacent binary64 numbers, thus $|q| \le 2^{-2}$. Similarly $q' = \circ(qz+Q_2)$ is rounded either to $Q_2 = 1/2$, or to one of the two adjacent binary64 numbers, thus $|q'| \le 1$. In $h_0 = \circ(q'z+Q_1)$, $q'z$ is so tiny that $h_0$ is either $Q_1 = 1$ or one of the two adjacent binary64 numbers, thus $|h_0| \le 2$. Then in $h_1$, $\ell_1 = \text{ExactMul}(z, h_0)$, we have $|h_1|, |\ell_1| \le 2^{-967.9}$. We thus perform the last operation $q_h, q_\ell = \text{FastSum}(Q_0, h_1, \ell_1)$, with $Q_0 = 1$ and $|h_1|, |\ell_1| \le 2^{-967.9}$. For rounding to nearest-even, it necessarily gives $q_h = 1$ and $|q_\ell| \le 2^{-965}$. Since we have $e_h = q_h = 1$ and $e_\ell = q_\ell$ at the end of exp_1, it is easy to check that the rounding test from Algorithm phase_1 will succeed, since $\circ(e_\ell \pm \epsilon e_h)$ rounds to $\pm \epsilon$, and $\circ(e_h \pm \epsilon)$ rounds to 1, which is the expected correct rounding for $\exp(y \log x)$ with tiny $y \log x$. For directed roundings, we will either

have $e_h = 1$ and tiny $e_\ell$, or $e_h = 1 - 2^{-53}$ and $e_\ell \approx 2^{-53}$, or $e_h = 1 + 2^{-52}$ and $e_\ell \approx -2^{-52}$. In the first case ($e_h = 1$ and tiny $e_\ell$), the rounding test will fail, since $\circ(e_\ell \pm \epsilon e_h)$ rounds to about $\pm \epsilon$, and $1 - \epsilon, 1 + \epsilon$ round to different values. In the second case ($e_h = 1 - 2^{-53}$ and $e_\ell \approx 2^{-53}$), $\circ(e_\ell \pm \epsilon e_h)$ rounds to about $2^{-53} \pm \epsilon$, thus we have approximately $u = \circ(1 - \epsilon)$ and $v = \circ(1 + \epsilon)$, which again round to different values. The third case is similar.

Note that this reasoning also holds when $z = 0$, which can happen when $x = 1$, for example. In this case we have $q_h = 1$, $q_\ell = 0$, thus $e_h = 1$, $e_\ell = 0$, $\circ(e_\ell \pm \epsilon e_h) = \pm \epsilon$, and the rounding test simplifies to $\circ(1 - \epsilon) = \circ(1 + \epsilon)$, which only holds for rounding to nearest.

In summary, when $|yh| < 2^{-969}$, Algorithm phase_1 returns the correct rounding for rounding to nearest-even, and returns FAIL for directed rounding modes, thus Theorem 1 holds.

**About the rounding test.** We know from (9) that

$$e_h + e_\ell - \tau e_h \le x^y \le e_h + e_\ell + \tau e_h,$$

from which we deduce by monotonicity of rounding:

$$\circ(e_h + e_\ell - \tau e_h) \le \circ(x^y) \le \circ(e_h + e_\ell + \tau e_h).$$

But since there is no instruction that directly computes $\circ(e_h + e_\ell \pm \tau e_h)$, in Algorithm phase_1, we instead compute $u := \circ(e_h + \circ(e_\ell - \epsilon e_h))$ and, similarly, $v := \circ(e_h + \circ(e_\ell + \epsilon e_h))$. We thus have an additional rounding error in $\circ(e_\ell \pm \epsilon e_h)$. Here, $\epsilon$ will be chosen slightly larger than $\tau$ (and thus a fortiori in $[\tau, 2\tau]$, say) in order to ensure that $u \le \circ(x^y) \le v$ despite these two extra rounding errors.

Let us now determine $\epsilon$. It suffices to have $\circ(e_\ell - \epsilon e_h) \le e_\ell - \tau e_h$ and $e_\ell + \tau e_h \le \circ(e_\ell + \epsilon e_h)$, since by monotonicity of rounding, this will ensure that $u \le \circ(x^y) \le v$. Now, defining for simplicity $u' := e_\ell - \epsilon e_h$ and $v' := e_\ell + \epsilon e_h$, we have $\circ(u') \le u' + \text{ulp}(u')$ and $v' - \text{ulp}(v') \le \circ(v')$, and so it suffices that $\epsilon$ be such that $-\epsilon e_h + \text{ulp}(u') \le -\tau e_h$ and $\tau e_h \le \epsilon e_h - \text{ulp}(v')$. Now, $|u'|, |v'| \le |e_\ell| + \epsilon e_h \le (2^{-49.2999} + 2\tau)e_h < 2^{-49.2905}e_h$ for $\epsilon \le 2\tau \le 2 \cdot 2^{-57.5604}$. If $|u'| \ge 2^{-1022}$ then $\text{ulp}(u') = 2^{-52}|u'| \le 2^{-101.2905}e_h$; if $|u'| < 2^{-1022}$, then $\text{ulp}(u') = 2^{-1074} \le 2^{-83}e_h$, since $e_h \ge 2^{-991}$. Consequently, even if underflow occurs, we have $\text{ulp}(u') \le 2^{-83}e_h$ and, similarly, $\text{ulp}(v') \le 2^{-83}e_h$. Hence it suffices that $\epsilon$ be such that $\epsilon \ge \tau + 2^{-83}$. It can then be checked that the values RU$(2^{-62.792})$ and RU$(2^{-57.560})$ taken for $\epsilon$ in Algorithm phase_1 satisfy this sufficient condition as well as the constraint $\epsilon \le 2\tau$.