From 3884948b7190448ba92eae205532a30637d0d810 Mon Sep 17 00:00:00 2001 From: dim Date: Wed, 4 Dec 2019 17:45:34 +0000 Subject: [PATCH] r355120 | dim | 2019-11-26 23:01:09 +0100 (Tue, 26 Nov 2019) | 32 lines The fdlibm hypot() implementations shouldn't potentially left-shift negative numbers (invoking undefined behavior) 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. Submitted by: Jeff Walden Obtained from: https://github.com/freebsd/freebsd/pull/414 Reviewed by: dim, lwhsu Differential Revision: https://reviews.freebsd.org/D22354 git-svn-id: svn://svn.freebsd.org/base/stable/10@355395 ccf9f872-aa2e-dd11-9fc8-001c23d0bc1f --- lib/msun/src/e_hypot.c | 6 ++---- lib/msun/src/e_hypotf.c | 2 +- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/lib/msun/src/e_hypot.c b/lib/msun/src/e_hypot.c index 2398e9837..812ba9b7d 100644 --- a/lib/msun/src/e_hypot.c +++ b/lib/msun/src/e_hypot.c @@ -118,10 +118,8 @@ __ieee754_hypot(double x, double y) w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b))); } if(k!=0) { - u_int32_t high; - t1 = 1.0; - GET_HIGH_WORD(high,t1); - SET_HIGH_WORD(t1,high+(k<<20)); + t1 = 0.0; + SET_HIGH_WORD(t1,(1023+k)<<20); return t1*w; } else return w; } diff --git a/lib/msun/src/e_hypotf.c b/lib/msun/src/e_hypotf.c index 6d083e4d2..efaf0dcd1 100644 --- a/lib/msun/src/e_hypotf.c +++ b/lib/msun/src/e_hypotf.c @@ -77,7 +77,7 @@ __ieee754_hypotf(float x, float y) w = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b))); } if(k!=0) { - SET_FLOAT_WORD(t1,0x3f800000+(k<<23)); + SET_FLOAT_WORD(t1,(127+k)<<23); return t1*w; } else return w; } -- 2.42.0