Do not left-shift a negative number (inducing undefined behavior in

C/C++) in exp(3), expf(3), expm1(3) and expm1f(3) during intermediate
computations that compute the IEEE-754 bit pattern for |2**k| for
integer |k|.

The implementations of exp(3), expf(3), expm1(3) and expm1f(3) need to
compute IEEE-754 bit patterns for 2**k in certain places.  (k is an
integer and 2**k is exactly representable in IEEE-754.)

Currently they do things like 0x3FF0'0000+(k<<20), which is to say they
take the bit pattern representing 1 and then add directly to the
exponent field to get the desired power of two.  This is fine when k is
non-negative.

But when k<0 (and certain classes of input trigger this), this
left-shifts a negative number -- an operation with undefined behavior in
C and C++.

The desired semantics can be achieved by instead adding the
possibly-negative k to the IEEE-754 exponent bias to get the desired
exponent field, _then_ shifting that into its proper overall position.

(Note that in case of s_expm1.c and s_expm1f.c, there are SET_HIGH_WORD
and SET_FLOAT_WORD uses further down in each of these files that perform
shift operations involving k, but by these points k's range has been
restricted to 2 < k <= 56, and the shift operations under those
circumstances can't do anything that would be UB.)

Submitted by:	Jeff Walden, https://github.com/jswalden
Obtained from:	https://github.com/freebsd/freebsd/pull/411
Obtained from:	https://github.com/freebsd/freebsd/pull/412
MFC after:	3 days
This commit is contained in:
Dimitry Andric 2019-09-25 18:50:57 +00:00
parent 618d66a56f
commit 5763a8cf06
4 changed files with 6 additions and 6 deletions

View File

@ -145,9 +145,9 @@ __ieee754_exp(double x) /* default IEEE double exp */
/* x is now in primary range */
t = x*x;
if(k >= -1021)
INSERT_WORDS(twopk,0x3ff00000+(k<<20), 0);
INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20, 0);
else
INSERT_WORDS(twopk,0x3ff00000+((k+1000)<<20), 0);
INSERT_WORDS(twopk,((u_int32_t)(0x3ff+(k+1000)))<<20, 0);
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
if(k==0) return one-((x*c)/(c-2.0)-x);
else y = one-((lo-(x*c)/(2.0-c))-hi);

View File

@ -83,9 +83,9 @@ __ieee754_expf(float x)
/* x is now in primary range */
t = x*x;
if(k >= -125)
SET_FLOAT_WORD(twopk,0x3f800000+(k<<23));
SET_FLOAT_WORD(twopk,((u_int32_t)(0x7f+k))<<23);
else
SET_FLOAT_WORD(twopk,0x3f800000+((k+100)<<23));
SET_FLOAT_WORD(twopk,((u_int32_t)(0x7f+(k+100)))<<23);
c = x - t*(P1+t*P2);
if(k==0) return one-((x*c)/(c-(float)2.0)-x);
else y = one-((lo-(x*c)/((float)2.0-c))-hi);

View File

@ -188,7 +188,7 @@ expm1(double x)
e = hxs*((r1-t)/(6.0 - x*t));
if(k==0) return x - (x*e-hxs); /* c is 0 */
else {
INSERT_WORDS(twopk,0x3ff00000+(k<<20),0); /* 2^k */
INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20,0); /* 2^k */
e = (x*(e-c)-c);
e -= hxs;
if(k== -1) return 0.5*(x-e)-0.5;

View File

@ -94,7 +94,7 @@ expm1f(float x)
e = hxs*((r1-t)/((float)6.0 - x*t));
if(k==0) return x - (x*e-hxs); /* c is 0 */
else {
SET_FLOAT_WORD(twopk,0x3f800000+(k<<23)); /* 2^k */
SET_FLOAT_WORD(twopk,((u_int32_t)(0x7f+k))<<23); /* 2^k */
e = (x*(e-c)-c);
e -= hxs;
if(k== -1) return (float)0.5*(x-e)-(float)0.5;