signed overflow in atan2
Bruce Evans
brde at optusnet.com.au
Wed Feb 14 05:34:44 UTC 2018
On Tue, 13 Feb 2018, Steve Kargl wrote:
> On Tue, Feb 13, 2018 at 08:23:01PM +1100, Bruce Evans wrote:
>> On Mon, 12 Feb 2018, Eitan Adler wrote:
>>
>>> source: https://github.com/freebsd/freebsd/pull/130
>>>
>>> `atan2(y, x)`'s special case for `x == 1.0` (in which case the result
>>> must be `atan(y)`) is unnecessarily complicated #130
>>>
>>> As a component of atan2(y, x), the case of x == 1.0 is farmed out to
>>> atan(y). The current implementation of this comparison is vulnerable
>>> to signed integer underflow (that is, undefined behavior), and it's
>>> performed in a somewhat more complicated way than it need be. Change
>>> it to not be quite so cute, rather directly comparing the high/low
>>> bits of x to the specific IEEE-754 bit pattern that encodes 1.0.
>>> ...
>>> -----
>>>
>>> patch:
>>>
>>> if(((ix|((lx|-lx)>>31))>0x7ff00000)||
>>> ((iy|((ly|-ly)>>31))>0x7ff00000)) /* x or y is NaN */
>>> return x+y;
>>> - if((hx-0x3ff00000|lx)==0) return atan(y); /* x=1.0 */
>>> + if(hx==0x3ff00000&&lx==0) return atan(y); /* x=1.0 */
>>> m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
>>> /* when y = 0 */
>>
>> Seems reasonable.
>>
>> Probably hx and hy should be uint32_t, but that would be a much larger
>> change (fdlibm+Cygnus fixes does it wrong in almost all places. Cygnus
>> changed lots of types to fixed-width, but didn't change many from signed
>> to unsigned). The next line that does ((hy>31&1) to extract the sign
>> bit is also unportable.
>
> I've had the following in my msun tree for awhile.
>
> --- math_private.h.orig 2017-05-30 09:22:20.565401000 -0700
> +++ math_private.h 2017-12-14 13:05:54.128219000 -0800
> @@ -84,7 +84,7 @@
>
> #endif
>
> -/* Get two 32 bit ints from a double. */
> +/* Get two unsigned 32 bit ints from a double. */
>
> #define EXTRACT_WORDS(ix0,ix1,d) \
> do { \
> @@ -166,8 +166,7 @@
> typedef union
> {
> float value;
> - /* FIXME: Assumes 32 bit int. */
> - unsigned int word;
> + uint32_t word;
> } ieee_float_shape_type;
>
> /* Get a 32 bit int from a float. */
So it is a type hack that is usually negatively useful to use signed
types for the output of these macros. (Usually, for EXTRACT_WORDS(hx,
lx, x), lx is u_int32_t but hx is int32_t so that "hx < 0" works right
but everything else tends to be unportable. The variable hx is
converted by assignment of the top word in the struct, so the sign
mismatch is hard to see. For ld80, the top bits in the struct are
only 16 bits and buggy code cloned from double precision might still
use int32_t for hx; then the conversion or rather the use of hx is more
magic.)
Bruce
More information about the freebsd-numerics
mailing list