From 6cf865dba69bab6346dc268d9173609af36b984e Mon Sep 17 00:00:00 2001 From: nsz Date: Sun, 6 May 2012 13:08:59 +0200 Subject: [PATCH] math: nextafter and nexttoward cleanup make nexttoward, nexttowardf independent of long double representation. fix nextafterl: it did not raise underflow flag when the result was 0. --- src/math/nextafter.c | 97 +++++++++----------------------- src/math/nextafterf.c | 84 +++++++++------------------- src/math/nextafterl.c | 123 ++++++++++++++++++++++------------------- src/math/nexttoward.c | 82 ++++++++++----------------- src/math/nexttowardf.c | 81 ++++++++++----------------- 5 files changed, 176 insertions(+), 291 deletions(-) diff --git a/src/math/nextafter.c b/src/math/nextafter.c index 5e53654a..d9e29236 100644 --- a/src/math/nextafter.c +++ b/src/math/nextafter.c @@ -1,79 +1,36 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_nextafter.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ -/* IEEE functions - * nextafter(x,y) - * return the next machine floating-point number of x in the - * direction toward y. - * Special cases: - */ - #include "libm.h" +#define SIGN ((uint64_t)1<<63) + double nextafter(double x, double y) { - volatile double t; - int32_t hx,hy,ix,iy; - uint32_t lx,ly; - - EXTRACT_WORDS(hx, lx, x); - EXTRACT_WORDS(hy, ly, y); - ix = hx & 0x7fffffff; /* |x| */ - iy = hy & 0x7fffffff; /* |y| */ + union dshape ux, uy; + uint64_t ax, ay; + int e; - if ((ix >= 0x7ff00000 && (ix-0x7ff00000)|lx) != 0 || /* x is nan */ - (iy >= 0x7ff00000 && (iy-0x7ff00000)|ly) != 0) /* y is nan */ - return x+y; - if (x == y) /* x == y */ + if (isnan(x) || isnan(y)) + return x + y; + ux.value = x; + uy.value = y; + if (ux.bits == uy.bits) return y; - if ((ix|lx) == 0) { /* x == 0 */ - INSERT_WORDS(x, hy&0x80000000, 1); /* return +-minsubnormal */ - /* raise underflow flag */ - t = x*x; - if (t == x) - return t; - return x; - } - if (hx >= 0) { /* x > 0 */ - if (hx > hy || (hx == hy && lx > ly)) { /* x > y, x -= ulp */ - if (lx == 0) - hx--; - lx--; - } else { /* x < y, x += ulp */ - lx++; - if (lx == 0) - hx++; - } - } else { /* x < 0 */ - if (hy >= 0 || hx > hy || (hx == hy && lx > ly)) { /* x < y, x -= ulp */ - if (lx == 0) - hx--; - lx--; - } else { /* x > y, x += ulp */ - lx++; - if (lx == 0) - hx++; - } - } - hy = hx & 0x7ff00000; - if (hy >= 0x7ff00000) /* overflow */ - return x+x; - if (hy < 0x00100000) { /* underflow */ - /* raise underflow flag */ - t = x*x; - if (t != x) { - INSERT_WORDS(y, hx, lx); + ax = ux.bits & ~SIGN; + ay = uy.bits & ~SIGN; + if (ax == 0) { + if (ay == 0) return y; - } + ux.bits = (uy.bits & SIGN) | 1; + } else if (ax > ay || ((ux.bits ^ uy.bits) & SIGN)) + ux.bits--; + else + ux.bits++; + e = ux.bits >> 52 & 0x7ff; + /* raise overflow if ux.value is infinite and x is finite */ + if (e == 0x7ff) + return x + x; + /* raise underflow if ux.value is subnormal or zero */ + if (e == 0) { + volatile double z = x*x + ux.value*ux.value; } - INSERT_WORDS(x, hx, lx); - return x; + return ux.value; } diff --git a/src/math/nextafterf.c b/src/math/nextafterf.c index bdc88cab..4727f7b0 100644 --- a/src/math/nextafterf.c +++ b/src/math/nextafterf.c @@ -1,67 +1,35 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_nextafterf.c */ -/* - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - #include "libm.h" +#define SIGN 0x80000000 + float nextafterf(float x, float y) { - volatile float t; - int32_t hx,hy,ix,iy; - - GET_FLOAT_WORD(hx, x); - GET_FLOAT_WORD(hy, y); - ix = hx & 0x7fffffff; /* |x| */ - iy = hy & 0x7fffffff; /* |y| */ + union fshape ux, uy; + uint32_t ax, ay, e; - if (ix > 0x7f800000 || /* x is nan */ - iy > 0x7f800000) /* y is nan */ - return x+y; - if (x == y) /* x == y */ + if (isnan(x) || isnan(y)) + return x + y; + ux.value = x; + uy.value = y; + if (ux.bits == uy.bits) return y; - if (ix == 0) { /* x == 0 */ - SET_FLOAT_WORD(x, (hy&0x80000000)|1); /* return +-minsubnormal */ - /* raise underflow flag */ - t = x*x; - if (t == x) - return t; - return x; - } - if (hx >= 0) { /* x > 0 */ - if (hx > hy) { /* x > y, x -= ulp */ - hx--; - } else { /* x < y, x += ulp */ - hx++; - } - } else { /* x < 0 */ - if (hy >= 0 || hx > hy) { /* x < y, x -= ulp */ - hx--; - } else { /* x > y, x += ulp */ - hx++; - } - } - hy = hx & 0x7f800000; - if (hy >= 0x7f800000) /* overflow */ - return x+x; - if (hy < 0x00800000) { /* underflow */ - /* raise underflow flag */ - t = x*x; - if (t != x) { - SET_FLOAT_WORD(y, hx); + ax = ux.bits & ~SIGN; + ay = uy.bits & ~SIGN; + if (ax == 0) { + if (ay == 0) return y; - } + ux.bits = (uy.bits & SIGN) | 1; + } else if (ax > ay || ((ux.bits ^ uy.bits) & SIGN)) + ux.bits--; + else + ux.bits++; + e = ux.bits & 0x7f800000; + /* raise overflow if ux.value is infinite and x is finite */ + if (e == 0x7f800000) + return x + x; + /* raise underflow if ux.value is subnormal or zero */ + if (e == 0) { + volatile float z = x*x + ux.value*ux.value; } - SET_FLOAT_WORD(x, hx); - return x; + return ux.value; } diff --git a/src/math/nextafterl.c b/src/math/nextafterl.c index aec8ab40..611ea53d 100644 --- a/src/math/nextafterl.c +++ b/src/math/nextafterl.c @@ -1,21 +1,3 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_nextafterl.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ -/* IEEE functions - * nextafter(x,y) - * return the next machine floating-point number of x in the - * direction toward y. - * Special cases: - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -23,55 +5,80 @@ long double nextafterl(long double x, long double y) { return nextafter(x, y); } -#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 +#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 +#define MSB ((uint64_t)1<<63) long double nextafterl(long double x, long double y) { - volatile long double t; - union IEEEl2bits ux, uy; - - ux.e = x; - uy.e = y; + union ldshape ux, uy; - if ((ux.bits.exp == 0x7fff && ((ux.bits.manh&~LDBL_NBIT)|ux.bits.manl) != 0) || - (uy.bits.exp == 0x7fff && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0)) - return x+y; /* x or y is nan */ + if (isnan(x) || isnan(y)) + return x + y; if (x == y) - return y; /* x=y, return y */ - if (x == 0.0) { - /* return +-minsubnormal */ - ux.bits.manh = 0; - ux.bits.manl = 1; + return y; + ux.value = x; + if (x == 0) { + uy.value = y; + ux.bits.m = 1; ux.bits.sign = uy.bits.sign; - /* raise underflow flag */ - t = ux.e*ux.e; - if (t == ux.e) - return t; - return ux.e; - } - if(x > 0.0 ^ x < y) { /* x -= ulp */ - if (ux.bits.manl == 0) { - if ((ux.bits.manh&~LDBL_NBIT) == 0) - ux.bits.exp--; - ux.bits.manh = (ux.bits.manh - 1) | (ux.bits.manh & LDBL_NBIT); + } else if (x < y ^ ux.bits.sign) { + ux.bits.m++; + if ((ux.bits.m & ~MSB) == 0) { + ux.bits.m = MSB; + ux.bits.exp++; } - ux.bits.manl--; - } else { /* x += ulp */ - ux.bits.manl++; - if (ux.bits.manl == 0) { - ux.bits.manh = (ux.bits.manh + 1) | (ux.bits.manh & LDBL_NBIT); - if ((ux.bits.manh&~LDBL_NBIT)==0) + } else { + if ((ux.bits.m & ~MSB) == 0) { + ux.bits.exp--; + if (ux.bits.exp) + ux.bits.m = 0; + } + ux.bits.m--; + } + /* raise overflow if ux.value is infinite and x is finite */ + if (ux.bits.exp == 0x7fff) + return x + x; + /* raise underflow if ux.value is subnormal or zero */ + if (ux.bits.exp == 0) { + volatile float z = x*x + ux.value*ux.value; + } + return ux.value; +} +#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 +long double nextafterl(long double x, long double y) +{ + union ldshape ux, uy; + + if (isnan(x) || isnan(y)) + return x + y; + if (x == y) + return y; + ux.value = x; + if (x == 0) { + uy.value = y; + ux.bits.mlo = 1; + ux.bits.sign = uy.bits.sign; + } else if (x < y ^ ux.bits.sign) { + ux.bits.mlo++; + if (ux.bits.mlo == 0) { + ux.bits.mhi++; + if (ux.bits.mhi == 0) ux.bits.exp++; } + } else { + if (ux.bits.mlo == 0) { + if (ux.bits.mhi == 0) + ux.bits.exp--; + ux.bits.mhi--; + } + ux.bits.mlo--; } - if (ux.bits.exp == 0x7fff) /* overflow */ - return x+x; - if (ux.bits.exp == 0) { /* underflow */ - mask_nbit_l(ux); - /* raise underflow flag */ - t = ux.e * ux.e; - if (t != ux.e) - return ux.e; + /* raise overflow if ux.value is infinite and x is finite */ + if (ux.bits.exp == 0x7fff) + return x + x; + /* raise underflow if ux.value is subnormal or zero */ + if (ux.bits.exp == 0) { + volatile float z = x*x + ux.value*ux.value; } - return ux.e; + return ux.value; } #endif diff --git a/src/math/nexttoward.c b/src/math/nexttoward.c index e150eaad..741b6b5f 100644 --- a/src/math/nexttoward.c +++ b/src/math/nexttoward.c @@ -1,15 +1,3 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_nexttoward.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -17,52 +5,42 @@ double nexttoward(double x, long double y) { return nextafter(x, y); } -#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 +#else +#define SIGN ((uint64_t)1<<63) + double nexttoward(double x, long double y) { - union IEEEl2bits uy; - volatile double t; - int32_t hx,ix; - uint32_t lx; + union dshape ux; + int e; - EXTRACT_WORDS(hx, lx, x); - ix = hx & 0x7fffffff; - uy.e = y; - - if ((ix >= 0x7ff00000 && ((ix-0x7ff00000)|lx) != 0) || - (uy.bits.exp == 0x7fff && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0)) - return x + y; /* x or y is nan */ + if (isnan(x) || isnan(y)) + return x + y; if (x == y) - return (double)y; - if (x == 0.0) { - INSERT_WORDS(x, uy.bits.sign<<31, 1); /* return +-minsubnormal */ - /* raise underflow */ - t = x * x; - if (t == x) - return t; - return x; - } - if (hx > 0.0 ^ x < y) { /* x -= ulp */ - if (lx == 0) - hx--; - lx--; - } else { /* x += ulp */ - lx++; - if (lx == 0) - hx++; + return y; + ux.value = x; + if (x == 0) { + ux.bits = 1; + if (signbit(y)) + ux.bits |= SIGN; + } else if (x < y) { + if (signbit(x)) + ux.bits--; + else + ux.bits++; + } else { + if (signbit(x)) + ux.bits++; + else + ux.bits--; } - ix = hx & 0x7ff00000; - if (ix >= 0x7ff00000) /* overflow */ + e = ux.bits>>52 & 0x7ff; + /* raise overflow if ux.value is infinite and x is finite */ + if (e == 0x7ff) return x + x; - if (ix < 0x00100000) { /* underflow */ - /* raise underflow flag */ - t = x * x; - if (t != x) { - INSERT_WORDS(x, hx, lx); - return x; - } + /* raise underflow if ux.value is subnormal or zero */ + if (e == 0) { + volatile float z = x*x + ux.value*ux.value; } - INSERT_WORDS(x, hx, lx); - return x; + return ux.value; } #endif diff --git a/src/math/nexttowardf.c b/src/math/nexttowardf.c index c52ef3aa..821f72a5 100644 --- a/src/math/nexttowardf.c +++ b/src/math/nexttowardf.c @@ -1,62 +1,37 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_nexttowardf.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - #include "libm.h" -#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 -// FIXME -#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 -#define LDBL_INFNAN_EXP (LDBL_MAX_EXP * 2 - 1) - float nexttowardf(float x, long double y) { - union IEEEl2bits uy; - volatile float t; - int32_t hx,ix; - - GET_FLOAT_WORD(hx, x); - ix = hx & 0x7fffffff; /* |x| */ - uy.e = y; + union fshape ux; + uint32_t e; - if (ix > 0x7f800000 || - (uy.bits.exp == LDBL_INFNAN_EXP && - ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl) != 0)) - return x + y; /* x or y is nan */ + if (isnan(x) || isnan(y)) + return x + y; if (x == y) - return (float)y; /* x=y, return y */ - if (ix == 0) { /* x == 0 */ - SET_FLOAT_WORD(x, (uy.bits.sign<<31)|1); /* return +-minsubnormal */ - /* raise underflow flag */ - t = x*x; - if (t == x) - return t; - return x; + return y; + ux.value = x; + if (x == 0) { + ux.bits = 1; + if (signbit(y)) + ux.bits |= 0x80000000; + } else if (x < y) { + if (signbit(x)) + ux.bits--; + else + ux.bits++; + } else { + if (signbit(x)) + ux.bits++; + else + ux.bits--; } - if (hx >= 0 ^ x < y) /* x -= ulp */ - hx--; - else /* x += ulp */ - hx++; - ix = hx & 0x7f800000; - if (ix >= 0x7f800000) /* overflow */ - return x+x; - if (ix < 0x00800000) { /* underflow */ - /* raise underflow flag */ - t = x*x; - if (t != x) { - SET_FLOAT_WORD(x, hx); - return x; - } + e = ux.bits & 0x7f800000; + /* raise overflow if ux.value is infinite and x is finite */ + if (e == 0x7f800000) + return x + x; + /* raise underflow if ux.value is subnormal or zero */ + if (e == 0) { + volatile float z = x*x + ux.value*ux.value; } - SET_FLOAT_WORD(x, hx); - return x; + return ux.value; } -#endif -- 2.20.1