Page MenuHomeFreeBSD

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

Authored by jwalden_mit.edu on Nov 13 2019, 10:29 PM.
Tags
None
Referenced Files
Unknown Object (File)
Tue, Dec 10, 9:48 PM
Unknown Object (File)
Mon, Nov 25, 9:20 PM
Unknown Object (File)
Sun, Nov 24, 6:21 PM
Unknown Object (File)
Sat, Nov 23, 9:39 PM
Unknown Object (File)
Sat, Nov 23, 8:33 PM
Unknown Object (File)
Nov 5 2024, 5:24 PM
Unknown Object (File)
Oct 21 2024, 8:18 PM
Unknown Object (File)
Oct 21 2024, 8:18 PM

Details

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.

Diff Detail

Repository
rS FreeBSD src repository - subversion
Lint
Lint Passed
Unit
No Test Coverage
Build Status
Buildable 27630
Build 25837: arc lint + arc unit

Event Timeline

This revision is now accepted and ready to land.Nov 16 2019, 2:41 PM
lwhsu requested changes to this revision.Nov 16 2019, 3:00 PM

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
This revision now requires changes to proceed.Nov 16 2019, 3:00 PM

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.

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.

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.

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

All tests under lib/msun pass.

This revision is now accepted and ready to land.Nov 24 2019, 7:34 PM