port-vax/59264: t_strtod:strtod_gherman_bug test is failing

[campbell+netbsd%mumble.net@localhost]

>Number:         59264
>Category:       port-vax
>Synopsis:       t_strtod:strtod_gherman_bug test is failing
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    port-vax-maintainer
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Mon Apr 07 00:50:00 +0000 2025
>Originator:     Taylor R Campbell
>Release:        current
>Organization:
The VaxBSD Floating-point Halfway House
>Environment:
>Description:
Termination reason

FAILED: 1 checks failed; see output for more details
Standard error stream

*** Check failed: /tmp/build/2025.04.04.21.52.19-vax/src/tests/lib/libc/stdlib/t_strtod.c:334: d != 0x1.d34fd8378ea83p+0: d=1.82544=0xe.9a7ec1bc7541cp-3

https://releng.netbsd.org/b5reports/vax/2025/2025.04.04.21.52.19/test.html#lib_libc_stdlib_t_strtod_strtod_gherman_bug

Relevant information about the bug in parsing the decimal notation for a number verrrrrrrrrrry close to halfway between two binary floating-point numbers:

https://www.exploringbinary.com/a-bug-in-the-bigcomp-function-of-david-gays-strtod/
https://web.archive.org/web/20241212101326/https://www.exploringbinary.com/a-bug-in-the-bigcomp-function-of-david-gays-strtod/

The notation is confusing because printf on VAX uses 0xe as the leading hex digit instead of 0x1, but once you shift it, you see the result is off by a single bit in the least significant place:

Expected: 0x1.d34fd8378ea830p+0
VAXctual: 0x1.d34fd8378ea838p+0
                           ^

Before the Geza Herman bug fix chronicled at https://www.exploringbinary.com/a-bug-in-the-bigcomp-function-of-david-gays-strtod/, strtod would erroneously return 0x1.d34fd8378ea84p+0 _in IEEE 754 binary64 arithmetic_.

But IEEE 754 binary64 arithmetic uses 53-bit precision, whereas VAX D (double) arithmetic uses 56-bit precision.

I used Sollya and MIT Scheme to verify that the answer we compute on VAX is correct:

;;; Sollya, with 1000-bit floating-point approximation
> prec=1000;
The precision has been set to 1000 bits.
> display = hexadecimal;
Display mode is hexadecimal numbers.
> round(1.8254370818746402660437411213933955878019332885742187, double, RN);
Warning: Rounding occurred when converting the constant "1.8254370818746402660437411213933955878019332885742187" to floating-point with 1000 bits.
If safe computation is needed, try to increase the precision.
0x1.d34fd8378ea83p0
> round(1.8254370818746402660437411213933955878019332885742187, 53, RN);
Warning: Rounding occurred when converting the constant "1.8254370818746402660437411213933955878019332885742187" to floating-point with 1000 bits.
If safe computation is needed, try to increase the precision.
0x1.d34fd8378ea83p0
> round(1.8254370818746402660437411213933955878019332885742187, 56, RN);
Warning: Rounding occurred when converting the constant "1.8254370818746402660437411213933955878019332885742187" to floating-point with 1000 bits.
If safe computation is needed, try to increase the precision.
0x1.d34fd8378ea838p0

;;; MIT Scheme, entirely with arbitrary-precision bignum integer arithmetic -- note that we start with a leading 1 bit, so the inner shift is by p-1 (p=53 for binary64, p=55 for VAX D), not by p; then the outer shift is to adjust the leading hexadecimal digit to match.
(number->string (* 16 (round (* #e1.8254370818746402660437411213933955878019332885742187 (expt 2 52)))) #x10)
;Value: "1d34fd8378ea830"

(number->string (* 2 (round (* #e1.8254370818746402660437411213933955878019332885742187 (expt 2 55)))) #x10)
;Value: "1d34fd8378ea838"
>How-To-Repeat:
cd /usr/tests/lib/libc/stdlib
atf-run t_strtod | atf-report
>Fix:
#if DBL_MANT_DIG == 53
	ATF_CHECK_EQ_MSG(d, 0x1.d34fd8378ea83p+0, "d=%g=%a", d, d);
#elif DBL_MANT_DIG == 56
	ATF_CHECK_EQ_MSG(d, 0x1.d34fd8378ea838p+0, "d=%g=%a", d, d);
#else
#  error Fill in the correct answer for this machine's DBL_MANT_DIG!
#endif