The fdlibm hypot() implementations shouldn't potentially left-shift negative numbers (invoking undefined behavior)

[jwalden_mit.edu]

Summary

Various paths through hypot(x, y) will multiply x and y by a power of two, perform the calculation in a range where IEEE-754 provides greater precision, then undo the multiplication to determine the true result. Undoing that multiplication is implemented as t1*w, where t1=2**k.

2**k is often computed by taking the high word of 1.0, then adding k<<20 (for doubles or long doubles) or k<<23 (for floats) to it, then overwriting that high word. But when k is negative this left-shifts a negative value -- and that's undefined behavior in many editions of C and C++.

This patch should fix all hypot implementations to compute 2**k without triggering this particular bit of undefined behavior.

Test Plan

I've only very lightly tested out the hypot(double, double) change, in SpiderMonkey's JavaScript engine, for consistency with prior behavior. The other functions' changes have more or less only been eyeballed. Careful examination appreciated! Do note, however, that an error in any of these changes would most likely produce a value that is incorrect by a factor of two, so any mistake would most likely be glaring if invoked.

[lwhsu]

https://lists.freebsd.org/pipermail/freebsd-numerics/2019-November/000923.html

This is reviewed by kargl@

[lwhsu]

One regression test fails:

lwhsu@x1c:/usr/tests/lib/msun > kyua debug csqrt_test:main
1..18
ok 1 - csqrt
ok 2 - csqrt
ok 3 - csqrt
ok 4 - csqrt
ok 5 - csqrt
ok 6 - csqrt
ok 7 - csqrt
ok 8 - csqrt
ok 9 - csqrt
ok 10 - csqrt
ok 11 - csqrt
ok 12 - csqrt
ok 13 - csqrt
ok 14 - csqrt
ok 15 - csqrt
ok 16 - csqrt
Assertion failed: (creall(result) == ldexpl(14 * 0x1p-4, exp / 2)), function test_overflow, file /usr/home/lwhsu/freebsd-src/lib/msun/tests/csqrt_test.c, line 236.
Process with PID 18700 exited with signal 6 and dumped core; attempting to gather stack trace
[New LWP 101179]
Core was generated by `/usr/tests/lib/msun/csqrt_test'.
Program terminated with signal SIGABRT, Aborted.
#0  thr_kill () at thr_kill.S:3
3       RSYSCALL(thr_kill)
#0  thr_kill () at thr_kill.S:3
#1  0x00000008004543a4 in __raise (s=6) at /usr/src/lib/libc/gen/raise.c:52
#2  0x00000008003c3029 in abort () at /usr/src/lib/libc/stdlib/abort.c:67
#3  0x0000000800440e71 in __assert (func=<optimized out>, file=<optimized out>, line=<optimized out>, failedexpr=<optimized out>) at /usr/src/lib/libc/gen/assert.c:51
#4  0x00000000002040bd in test_overflow (maxexp=16384) at /usr/home/lwhsu/freebsd-src/lib/msun/tests/csqrt_test.c:236
#5  0x0000000000202a91 in main () at /usr/home/lwhsu/freebsd-src/lib/msun/tests/csqrt_test.c:363
GDB exited successfully
Files left in work directory after failure: csqrt_test.core
csqrt_test:main  ->  broken: Received signal 6

[dim]

In D22354#489981, @lwhsu wrote:

Assertion failed: (creall(result) == ldexpl(14 * 0x1p-4, exp / 2)), function test_overflow, file /usr/home/lwhsu/freebsd-src/lib/msun/tests/csqrt_test.c, line 236.

Indeed, and this test passed before the patch was applied. I'm taking a look at the generated code.

[dim]

lib/msun/src/e_hypotl.c

118 ↗

(On Diff #64291)

Why was this changed from 1.0 to 0? Steve Kargl posted on https://lists.freebsd.org/pipermail/freebsd-numerics/2019-November/000925.html:

This is the original code

     u_int32_t high;
     t1 = 1.0;
     GET_HIGH_WORD(high,t1);
     SET_HIGH_WORD(t1,high+DESW(k));

high + DESW(k) = high + k
               = 16383 + k

and this is the code after the patch

     t1 = 0.0;
     SET_HIGH_WORD(t1,ESW(k));

ESW(k) = MAX_EXP - 1 + k
       = LDBL_MAX_EXP - 1 + k
       = 16384 - 1 + k
       = 16383 + k

So, in principle there is no functional change.

but changing t1 from 1 to 0 does constitute a functional change?

When I put it back to 1.0, but leave the rest of the patch in, all tests succeed, including the csqrt test.

[jwalden_mit.edu]

Removed the e_hypotl.c changes

My understanding -- which is evidently wrong? -- was that long double had an IEEE-754ish representation, and any power of two -- so long as it wouldn't have subnormal representation -- would only have bits in an exponent field. 1.0 isn't subnormal, so it wouldn't have any significand bits set.

And the line if(ha > ESW(MAX_EXP/2-12)) { /* a>2**(MAX_EXP/2-12) */, where ha was a high word -- indicated to me that ESW(k) would be the right thing to use to set *specifically the exponent bits*.

Seems I'm mistaken here...somehow. (I work in lands of IEEE-754 single- and double-precision, so I generally have no need to know much about long double.)

Anyway, if all that's not the case for some reason here, and my attempt to simplify things here (beyond the avoid-UB reason for the overall change) went awry for reasons unclear to me, let's just withdraw the long double portion of this change.

[dim]

Ok, if the long double part is dropped, the tests succeed. So LGTM.

[lwhsu]

All tests under lib/msun pass.