Proofs: Clinger's correct decimal input

Clinger's fast path and AlgorithmM return the correctly rounded binary64 number full

This proof shows: correctness.

Setting. Binary64: precision $P = 53$, the smallest subnormal $2^{K_{\min}}$ with $K_{\min} = -1074$, the largest finite number $(2^{P} - 1)\,2^{K_{\max}}$ with $K_{\max} = 971$. $\mathrm{RN}$ is rounding to nearest, ties to even, with the overflow rule of IEEE 754-2019 (Clause 7.4): the result is $\infty$ when the value rounded with an unbounded exponent range is at least $2^{1024}$. The input is a pair of integers $f > 0$ and $e$; the goal is $\mathrm{RN}(f \cdot 10^{e})$. Two facts about the format are used. (F1) If $k \ge K_{\min}$, then the floats in $[2^{k + P - 1}, 2^{k + P}]$ are exactly the integer multiples of $2^{k}$ in that interval (for $k + P - 1 \ge -1022$ this is the binade of exponent $k + P - 1$ and its upper end). (F2) The floats in $[0, 2^{-1022}]$ are exactly the integer multiples of $2^{K_{\min}}$ in that interval (the subnormals, zero and the smallest normal).

Lemma 1 (fast path, lines L2 and L3). If $f < 2^{53}$ and $|e| \le 22$, then L3 returns $\mathrm{RN}(f \cdot 10^{e})$.

Proof. $f$ is an integer below $2^{53}$, so $\mathrm{float}(f) = f$ exactly. For $0 \le j \le 22$, $10^{j} = 5^{j} 2^{j}$ and $5^{j} \le 5^{22} < 2^{53}$, so $10^{j}$ is a binary64 number and the table entry is exact. IEEE 754 multiplication and division return the correctly rounded value of the exact product or quotient of their operands (Clause 5.4.1), so the single operation returns $\mathrm{RN}(f \cdot 10^{j})$ or $\mathrm{RN}(f / 10^{j})$. No overflow or underflow occurs: $f \cdot 10^{22} < 2^{53} \cdot 10^{22} < 2^{127}$ and $f / 10^{22} \ge 10^{-22} > 2^{-74}$. $\square$

Lemma 2 (scale, lines L5 and L6). Let $u > 0$ and $v > 0$ be integers, $\ell_u$ and $\ell_v$ their lengths in binary, and $x = u/v$. After L6, $2^{P - 1} \le x / 2^{k} < 2^{P}$.

Proof. From $2^{\ell_u - 1} \le u < 2^{\ell_u}$ and $2^{\ell_v - 1} \le v < 2^{\ell_v}$ we get $2^{\ell_u - \ell_v - 1} < x < 2^{\ell_u - \ell_v + 1}$, so with $k = \ell_u - \ell_v - P$ (L5), $2^{P - 1} < x/2^{k} < 2^{P + 1}$. L6 computes $\lfloor x / 2^{k} \rfloor$. If it is at least $2^{P}$, then $x / 2^{k} \ge 2^{P}$ and after $k \leftarrow k + 1$ we have $2^{P - 1} \le x / 2^{k} < 2^{P}$. Otherwise $x / 2^{k} < 2^{P}$ already, since $x/2^k < \lfloor x/2^k \rfloor + 1 \le 2^{P}$. $\square$

Lemma 3 (the neighbors of $x$). After L7, set $q_0 = \lfloor x / 2^{k} \rfloor$. Then $q_0 2^{k}$ and $(q_0 + 1) 2^{k}$ are consecutive numbers of the binary64 grid with an unbounded exponent range, and $q_0 2^{k} \le x < (q_0 + 1) 2^{k}$.

Proof. The inequality holds by the definition of the floor. If L7 did not change $k$, then $k \ge K_{\min}$ and $2^{P - 1} \le q_0 < 2^{P}$ by Lemma 2, so both points lie in $[2^{k + P - 1}, 2^{k + P}]$, and by (F1), applied to the format without an upper exponent limit, the grid points there are exactly the multiples of $2^{k}$. If L7 raised $k$ to $K_{\min}$, then by Lemma 2 the old scale $k' < K_{\min}$ gives $x < 2^{k' + P} \le 2^{K_{\min} - 1 + P} = 2^{-1022}$, so $q_0 < 2^{P - 1}$ and both points lie in $[0, 2^{-1022}]$, where by (F2) the floats are exactly the multiples of $2^{K_{\min}}$. In both cases there is no grid point strictly between them. $\square$

Lemma 4 (rounding, lines L8 and L9). After L9, $q 2^{k}$ is the grid point nearest to $x$, with ties broken to the even significand.

Proof. L8 divides the integers $\mathit{num}$ and $\mathit{den}$ with $\mathit{num}/\mathit{den} = x / 2^{k}$, so $q_0$ is the quotient and $x / 2^{k} = q_0 + r/\mathit{den}$ with $0 \le r < \mathit{den}$. By Lemma 3 the nearest grid point is $q_0 2^{k}$ if $r/\mathit{den} < 1/2$, that is $2r < \mathit{den}$, and $(q_0 + 1) 2^{k}$ if $2r > \mathit{den}$. If $2r = \mathit{den}$, then $x$ is the midpoint, and IEEE ties to even picks the neighbor whose integral significand is even. For $q_0 \ge 2^{P - 1}$ (normal case) and for $q_0 < 2^{P - 1}$ (subnormal case) the integral significand of $q_0 2^{k}$ is $q_0$ itself; if $q_0 + 1 = 2^{P}$, the upper neighbor is the power of two $2^{P - 1} \cdot 2^{k + 1}$ with even significand $2^{P - 1}$, and $q_0 = 2^P - 1$ is odd; if $q_0 + 1 = 2^{P - 1}$ in the subnormal case, the upper neighbor is $2^{-1022}$, whose significand $2^{P-1}$ is even, and $q_0$ is odd. In every case "the even one" is $q_0$ when $q_0$ is even and $q_0 + 1$ otherwise, which is the test of L9. The carry step rewrites $2^{P} \cdot 2^{k}$ as $2^{P - 1} \cdot 2^{k + 1}$ without changing the value. $\square$

Theorem. For integers $f > 0$ and $e$, clinger_read(f, e) returns $\mathrm{RN}(f \cdot 10^{e})$, including $+\infty$ on overflow and subnormal or zero results on underflow. With the sign applied afterwards, read_float returns the correctly rounded value of any decimal string it accepts, since $\mathrm{RN}(-y) = -\mathrm{RN}(y)$.

Proof. If L2 holds, Lemma 1 applies. Otherwise L4 forms $u$ and $v$ with $u / v = f \cdot 10^{e}$ exactly: $u = f \cdot 10^{e}$ and $v = 1$ for $e \ge 0$, $u = f$ and $v = 10^{-e}$ for $e < 0$. By Lemma 4, after L9 the value $q 2^{k}$ is the rounding of $x$ to the binary64 grid with an unbounded exponent range, with $q \le 2^{P} - 1$ after the carry step and $k \ge K_{\min}$. If $k \le K_{\max}$, then $q 2^{k} \le (2^{P} - 1) 2^{K_{\max}}$ is a finite float, and ldexp returns it exactly. If $k > K_{\max}$, then $q \ge 2^{P - 1}$ (the case $q < 2^{P-1}$ occurs only with $k = K_{\min}$), so $q 2^{k} \ge 2^{P - 1 + K_{\max} + 1} = 2^{1024}$, and the IEEE overflow rule for rounding to nearest gives $+\infty$, which L10 returns. $\blacksquare$

Remarks. The proof uses only exact integer arithmetic in AlgorithmM, so it holds for any number of digits and any exponent; the cost, not the correctness, grows with them. Lemma 1 is the only place where floating-point operations are used, and it needs the processor to round binary64 operations once, to binary64 precision; on x87 hardware with 64-bit internal precision the fast path can round twice (double rounding).