diff options
| -rw-r--r-- | numpy/core/src/npymath/ieee754.c.src | 90 | ||||
| -rw-r--r-- | numpy/core/src/npymath/npy_math_private.h | 309 |
2 files changed, 211 insertions, 188 deletions
diff --git a/numpy/core/src/npymath/ieee754.c.src b/numpy/core/src/npymath/ieee754.c.src index 6a27ee58c..a8db2361b 100644 --- a/numpy/core/src/npymath/ieee754.c.src +++ b/numpy/core/src/npymath/ieee754.c.src @@ -1,4 +1,7 @@ +/* -*- c -*- */ /* + * vim:syntax=c + * * Low-level routines related to IEEE-754 format */ #include "npy_math_common.h" @@ -159,66 +162,67 @@ float npy_nextafterf(float x, float y) npy_longdouble npy_nextafterl(npy_longdouble x, npy_longdouble y) { volatile npy_longdouble t; - ldouble_man_t xmanl, xmanh, ymanl, ymanh; - ldouble_exp_t xexp, yexp; - ldouble_sign_t xsign, ysign; - - GET_LDOUBLE_SIGN(x, xsign); - GET_LDOUBLE_EXP(x, xexp); - GET_LDOUBLE_MANL(x, xmanl); - GET_LDOUBLE_MANH(x, xmanh); + union IEEEl2bitsrep ux; + union IEEEl2bitsrep uy; - GET_LDOUBLE_SIGN(y, ysign); - GET_LDOUBLE_EXP(y, yexp); - GET_LDOUBLE_MANL(y, ymanl); - GET_LDOUBLE_MANH(y, ymanh); + ux.e = x; + uy.e = y; - if ((xexp == 0x7fff && ((xmanh & ~LDBL_NBIT) | xmanl) != 0) || - (yexp == 0x7fff && ((ymanh & ~LDBL_NBIT) | ymanl) != 0)) { - return x + y; /* x or y is nan */ + if ((GET_LDOUBLE_EXP(ux) == 0x7fff && + ((GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT) | GET_LDOUBLE_MANL(ux)) != 0) || + (GET_LDOUBLE_EXP(uy) == 0x7fff && + ((GET_LDOUBLE_MANH(uy) & ~LDBL_NBIT) | GET_LDOUBLE_MANL(uy)) != 0)) { + return ux.e + uy.e; /* x or y is nan */ } - if (x == y) { - return y; /* x=y, return y */ + if (ux.e == uy.e) { + return uy.e; /* x=y, return y */ } - if (x == 0.0) { - xmanh = 0; /* return +-minsubnormal */ - xmanl = 1; - xsign = ysign; - t = x * x; - if (t == x) { + if (ux.e == 0.0) { + SET_LDOUBLE_MANH(ux, 0); /* return +-minsubnormal */ + SET_LDOUBLE_MANL(ux, 1); + SET_LDOUBLE_SIGN(ux, GET_LDOUBLE_SIGN(uy)); + t = ux.e * ux.e; + if (t == ux.e) { return t; } else { - return x; /* raise underflow flag */ + return ux.e; /* raise underflow flag */ } } - if (x > 0.0 ^ x < y) { /* x -= ulp */ - if (xmanl == 0) { - if ((xmanh & ~LDBL_NBIT) == 0) { - xexp -= 1; + if ((ux.e > 0.0) ^ (ux.e < uy.e)) { /* x -= ulp */ + if (GET_LDOUBLE_MANL(ux) == 0) { + if ((GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT) == 0) { + SET_LDOUBLE_EXP(ux, GET_LDOUBLE_EXP(ux) - 1); } - xmanh = (xmanh - 1) | (xmanh & LDBL_NBIT); + SET_LDOUBLE_MANH(ux, + (GET_LDOUBLE_MANH(ux) - 1) | + (GET_LDOUBLE_MANH(ux) & LDBL_NBIT)); } - xmanl -= 1; + SET_LDOUBLE_MANL(ux, GET_LDOUBLE_MANL(ux) - 1); } else { /* x += ulp */ - xmanl += 1; - if (xmanl == 0) { - xmanh = (xmanh + 1) | (xmanh & LDBL_NBIT); - if ((xmanh & ~LDBL_NBIT) == 0) { - xexp += 1; + SET_LDOUBLE_MANL(ux, GET_LDOUBLE_MANL(ux) + 1); + if (GET_LDOUBLE_MANL(ux) == 0) { + SET_LDOUBLE_MANH(ux, + (GET_LDOUBLE_MANH(ux) + 1) | + (GET_LDOUBLE_MANH(ux) & LDBL_NBIT)); + if ((GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT) == 0) { + SET_LDOUBLE_EXP(ux, GET_LDOUBLE_EXP(ux) + 1); } } } - if (xexp == 0x7fff) { - return x + x; /* overflow */ + if (GET_LDOUBLE_EXP(ux) == 0x7fff) { + return ux.e + ux.e; /* overflow */ } - if (xexp == 0) { /* underflow */ - mask_nbit_l(xmanh); - t = x * x; - if (t != x) { /* raise underflow flag */ - return x; + if (GET_LDOUBLE_EXP(ux) == 0) { /* underflow */ + if (LDBL_NBIT) { + SET_LDOUBLE_MANH(ux, GET_LDOUBLE_MANH(ux) & ~LDBL_NBIT); + } + t = ux.e * ux.e; + if (t != ux.e) { /* raise underflow flag */ + return ux.e; } } - return x; + + return ux.e; } #endif diff --git a/numpy/core/src/npymath/npy_math_private.h b/numpy/core/src/npymath/npy_math_private.h index 9880b5ed8..a5a46c514 100644 --- a/numpy/core/src/npymath/npy_math_private.h +++ b/numpy/core/src/npymath/npy_math_private.h @@ -129,12 +129,12 @@ do { \ /* Set a double from two 32 bit ints. */ -#define INSERT_WORDS(d,ix0,ix1) \ -do { \ - ieee_double_shape_type iw_u; \ - iw_u.parts.msw = (ix0); \ - iw_u.parts.lsw = (ix1); \ - (d) = iw_u.value; \ +#define INSERT_WORDS(d,ix0,ix1) \ +do { \ + ieee_double_shape_type iw_u; \ + iw_u.parts.msw = (ix0); \ + iw_u.parts.lsw = (ix1); \ + (d) = iw_u.value; \ } while (0) /* @@ -151,20 +151,20 @@ typedef union /* Get a 32 bit int from a float. */ -#define GET_FLOAT_WORD(i,d) \ -do { \ - ieee_float_shape_type gf_u; \ - gf_u.value = (d); \ - (i) = gf_u.word; \ +#define GET_FLOAT_WORD(i,d) \ +do { \ + ieee_float_shape_type gf_u; \ + gf_u.value = (d); \ + (i) = gf_u.word; \ } while (0) /* Set a float from a 32 bit int. */ -#define SET_FLOAT_WORD(d,i) \ -do { \ - ieee_float_shape_type sf_u; \ - sf_u.word = (i); \ - (d) = sf_u.value; \ +#define SET_FLOAT_WORD(d,i) \ +do { \ + ieee_float_shape_type sf_u; \ + sf_u.word = (i); \ + (d) = sf_u.value; \ } while (0) #ifdef NPY_USE_C99_COMPLEX @@ -183,65 +183,69 @@ do { \ * * 16 stronger bits of a[2] are junk */ + typedef npy_uint32 IEEEl2bitsrep_part; + +/* my machine */ + union IEEEl2bitsrep { - npy_longdouble e; - npy_uint32 a[3]; + npy_longdouble e; + IEEEl2bitsrep_part a[3]; }; #define LDBL_MANL_INDEX 0 #define LDBL_MANL_MASK 0xFFFFFFFF - #define LDBL_MANL_SHIFT(x) (x) + #define LDBL_MANL_SHIFT 0 #define LDBL_MANH_INDEX 1 #define LDBL_MANH_MASK 0xFFFFFFFF - #define LDBL_MANH_SHIFT(x) (x) + #define LDBL_MANH_SHIFT 0 #define LDBL_EXP_INDEX 2 #define LDBL_EXP_MASK 0x7FFF - #define LDBL_EXP_SHIFT(x) (x) + #define LDBL_EXP_SHIFT 0 #define LDBL_SIGN_INDEX 2 #define LDBL_SIGN_MASK 0x8000 - #define LDBL_SIGN_SHIFT(x) ((x) >> 15) + #define LDBL_SIGN_SHIFT 15 #define LDBL_NBIT 0x80000000 - #define mask_nbit_l(manh) ((manh) &= ~LDBL_NBIT) typedef npy_uint32 ldouble_man_t; typedef npy_uint32 ldouble_exp_t; typedef npy_uint32 ldouble_sign_t; #elif defined(HAVE_LDOUBLE_INTEL_EXTENDED_16_BYTES_LE) - /* + /* * Intel extended 80 bits precision, 16 bytes alignment.. Bit representation is - * | junk | s |eeeeeeeeeeeeeee|mmmmmmmm................mmmmmmm| - * | 16 bits| 1 bit | 15 bits | 64 bits | - * | a[2] | a[1] | a[0] | + * | junk | s |eeeeeeeeeeeeeee|mmmmmmmm................mmmmmmm| + * | 16 bits| 1 bit | 15 bits | 64 bits | + * | a[2] | a[1] | a[0] | * * a[3] and 16 stronger bits of a[2] are junk */ + typedef npy_uint32 IEEEl2bitsrep_part; + union IEEEl2bitsrep { - npy_longdouble e; - npy_uint32 a[4]; + npy_longdouble e; + IEEEl2bitsrep_part a[4]; }; #define LDBL_MANL_INDEX 0 - #define LDBL_MANL_MASK 0xFFFFFFFF - #define LDBL_MANL_SHIFT(x) ((x)) - + #define LDBL_MANL_MASK 0xFFFFFFFF + #define LDBL_MANL_SHIFT 0 + #define LDBL_MANH_INDEX 1 - #define LDBL_MANH_MASK 0xFFFFFFFF - #define LDBL_MANH_SHIFT(x) ((x)) - - #define LDBL_EXP_INDEX 2 - #define LDBL_EXP_MASK 0x7FFF - #define LDBL_EXP_SHIFT(x) (x) - + #define LDBL_MANH_MASK 0xFFFFFFFF + #define LDBL_MANH_SHIFT 0 + + #define LDBL_EXP_INDEX 2 + #define LDBL_EXP_MASK 0x7FFF + #define LDBL_EXP_SHIFT 0 + #define LDBL_SIGN_INDEX 2 - #define LDBL_SIGN_MASK 0x8000 - #define LDBL_SIGN_SHIFT(x) ((x) >> 15) + #define LDBL_SIGN_MASK 0x8000 + #define LDBL_SIGN_SHIFT 15 - #define LDBL_NBIT 0x80000000 - #define mask_nbit_l(manh) ((manh) &= ~LDBL_NBIT) + #define LDBL_NBIT 0x800000000 typedef npy_uint32 ldouble_man_t; typedef npy_uint32 ldouble_exp_t; @@ -251,35 +255,36 @@ do { \ /* 64 bits IEEE double precision aligned on 16 bytes: used by ppc arch on * Mac OS X */ - /* + /* * IEEE double precision. Bit representation is - * | s |eeeeeeeeeee|mmmmmmmm................mmmmmmm| - * |1 bit| 11 bits | 52 bits | - * | a[0] | a[1] | + * | s |eeeeeeeeeee|mmmmmmmm................mmmmmmm| + * |1 bit| 11 bits | 52 bits | + * | a[0] | a[1] | */ + typedef npy_uint32 IEEEl2bitsrep_part; + union IEEEl2bitsrep { - long double e; - npy_uint32 a[2]; + npy_longdouble e; + IEEEl2bitsrep_part a[2]; }; - + #define LDBL_MANL_INDEX 1 - #define LDBL_MANL_MASK 0xFFFFFFFF - #define LDBL_MANL_SHIFT(x) ((x)) - + #define LDBL_MANL_MASK 0xFFFFFFFF + #define LDBL_MANL_SHIFT(x) 0 + #define LDBL_MANH_INDEX 0 - #define LDBL_MANH_MASK 0x000FFFFF - #define LDBL_MANH_SHIFT(x) ((x)) - - #define LDBL_EXP_INDEX 0 - #define LDBL_EXP_MASK 0x7FF00000 - #define LDBL_EXP_SHIFT(x) ((x) >> 20) - + #define LDBL_MANH_MASK 0x000FFFFF + #define LDBL_MANH_SHIFT(x) 0 + + #define LDBL_EXP_INDEX 0 + #define LDBL_EXP_MASK 0x7FF00000 + #define LDBL_EXP_SHIFT(x) 20 + #define LDBL_SIGN_INDEX 0 - #define LDBL_SIGN_MASK 0x80000000 - #define LDBL_SIGN_SHIFT(x) ((x) >> 31) + #define LDBL_SIGN_MASK 0x80000000 + #define LDBL_SIGN_SHIFT(x) 31 #define LDBL_NBIT 0 - #define mask_nbit_l(u) ((void)0) typedef npy_uint32 ldouble_man_t; typedef npy_uint32 ldouble_exp_t; @@ -287,137 +292,151 @@ do { \ #elif defined(HAVE_LDOUBLE_IEEE_DOUBLE_LE) /* 64 bits IEEE double precision, Little Endian. */ - /* + /* * IEEE double precision. Bit representation is - * | s |eeeeeeeeeee|mmmmmmmm................mmmmmmm| - * |1 bit| 11 bits | 52 bits | - * | a[1] | a[0] | + * | s |eeeeeeeeeee|mmmmmmmm................mmmmmmm| + * |1 bit| 11 bits | 52 bits | + * | a[1] | a[0] | */ + typedef npy_uint32 IEEEl2bitsrep_part; + union IEEEl2bitsrep { - npy_longdouble e; - npy_uint32 a[2]; + npy_longdouble e; + IEEEl2bitsrep_part a[2]; }; - + #define LDBL_MANL_INDEX 0 - #define LDBL_MANL_MASK 0xFFFFFFFF - #define LDBL_MANL_SHIFT(x) ((x)) - + #define LDBL_MANL_MASK 0xFFFFFFFF + #define LDBL_MANL_SHIFT 0 + #define LDBL_MANH_INDEX 1 - #define LDBL_MANH_MASK 0x000FFFFF - #define LDBL_MANH_SHIFT(x) ((x)) - - #define LDBL_EXP_INDEX 1 - #define LDBL_EXP_MASK 0x7FF00000 - #define LDBL_EXP_SHIFT(x) ((x) >> 20) - + #define LDBL_MANH_MASK 0x000FFFFF + #define LDBL_MANH_SHIFT 0 + + #define LDBL_EXP_INDEX 1 + #define LDBL_EXP_MASK 0x7FF00000 + #define LDBL_EXP_SHIFT 20 + #define LDBL_SIGN_INDEX 1 - #define LDBL_SIGN_MASK 0x80000000 - #define LDBL_SIGN_SHIFT(x) ((x) >> 31) + #define LDBL_SIGN_MASK 0x80000000 + #define LDBL_SIGN_SHIFT 31 - #define LDBL_NBIT 0x80000000 - #define mask_nbit_l(manh) ((manh) &= ~LDBL_NBIT) + #define LDBL_NBIT 0x00000080 typedef npy_uint32 ldouble_man_t; typedef npy_uint32 ldouble_exp_t; typedef npy_uint32 ldouble_sign_t; #elif defined(HAVE_LDOUBLE_IEEE_QUAD_BE) - /* + /* * IEEE quad precision, Big Endian. Bit representation is - * | s |eeeeeeeeeee|mmmmmmmm................mmmmmmm| - * |1 bit| 15 bits | 112 bits | - * | a[0] | a[1] | + * | s |eeeeeeeeeee|mmmmmmmm................mmmmmmm| + * |1 bit| 15 bits | 112 bits | + * | a[0] | a[1] | */ + typedef npy_uint64 IEEEl2bitsrep_part; + union IEEEl2bitsrep { - npy_longdouble e; - npy_uint64 a[2]; + npy_longdouble e; + IEEEl2bitsrep_part a[2]; }; - + #define LDBL_MANL_INDEX 1 - #define LDBL_MANL_MASK 0xFFFFFFFFFFFFFFFF - #define LDBL_MANL_SHIFT(x) ((x)) - + #define LDBL_MANL_MASK 0xFFFFFFFFFFFFFFFF + #define LDBL_MANL_SHIFT 0 + #define LDBL_MANH_INDEX 0 - #define LDBL_MANH_MASK 0x0000FFFFFFFFFFFF - #define LDBL_MANH_SHIFT(x) ((x)) - - #define LDBL_EXP_INDEX 0 - #define LDBL_EXP_MASK 0x7FFF000000000000 - #define LDBL_EXP_SHIFT(x) ((x) >> 48) - + #define LDBL_MANH_MASK 0x0000FFFFFFFFFFFF + #define LDBL_MANH_SHIFT 0 + + #define LDBL_EXP_INDEX 0 + #define LDBL_EXP_MASK 0x7FFF000000000000 + #define LDBL_EXP_SHIFT 48 + #define LDBL_SIGN_INDEX 0 - #define LDBL_SIGN_MASK 0x8000000000000000 - #define LDBL_SIGN_SHIFT(x) ((x) >> 63) + #define LDBL_SIGN_MASK 0x8000000000000000 + #define LDBL_SIGN_SHIFT 63 - #define mask_nbit_l(u) ((void)0) + #define LDBL_NBIT 0 typedef npy_uint64 ldouble_man_t; typedef npy_uint64 ldouble_exp_t; typedef npy_uint32 ldouble_sign_t; #endif -/* Put the sign bit of x into i - i should be defined as ldouble_sign_t */ -#define GET_LDOUBLE_SIGN(x, i) \ -do { \ - const union IEEEl2bitsrep u = {x}; \ - (i) = LDBL_SIGN_SHIFT(u.a[LDBL_SIGN_INDEX] & LDBL_SIGN_MASK); \ -} while (0) - -/* Put the exp bits of x into i - i should be defined as ldouble_exp_t */ -#define GET_LDOUBLE_EXP(x, i) \ -do { \ - const union IEEEl2bitsrep u = {x}; \ - (i) = LDBL_EXP_SHIFT(u.a[LDBL_EXP_INDEX] & LDBL_EXP_MASK); \ -} while (0) - -/* Put the manl bits of x into i - i should be defined as ldouble_man_t */ -#define GET_LDOUBLE_MANL(x, i) \ -do { \ - const union IEEEl2bitsrep u = {x}; \ - (i) = LDBL_MANL_SHIFT(u.a[LDBL_MANL_INDEX] & LDBL_MANL_MASK); \ -} while (0) - -/* Put the manh bits (hist bits of mantisss) of x into i - i should be defined - * as ldouble_man_t */ -#define GET_LDOUBLE_MANH(x, i) \ -do { \ - const union IEEEl2bitsrep u = {x}; \ - (i) = LDBL_MANH_SHIFT(u.a[LDBL_MANH_INDEX] & LDBL_MANH_MASK); \ -} while (0) +/* Get the sign bit of x. x should be of type IEEEl2bitsrep */ +#define GET_LDOUBLE_SIGN(x) \ + (((x).a[LDBL_SIGN_INDEX] & LDBL_SIGN_MASK) >> LDBL_SIGN_SHIFT) + +/* Set the sign bit of x to v. x should be of type IEEEl2bitsrep */ +#define SET_LDOUBLE_SIGN(x, v) \ + ((x).a[LDBL_SIGN_INDEX] = \ + ((x).a[LDBL_SIGN_INDEX] & ~LDBL_SIGN_MASK) | \ + (((IEEEl2bitsrep_part)(v) << LDBL_SIGN_SHIFT) & LDBL_SIGN_MASK)) + +/* Get the exp bits of x. x should be of type IEEEl2bitsrep */ +#define GET_LDOUBLE_EXP(x) \ + (((x).a[LDBL_EXP_INDEX] & LDBL_EXP_MASK) >> LDBL_EXP_SHIFT) + +/* Set the exp bit of x to v. x should be of type IEEEl2bitsrep */ +#define SET_LDOUBLE_EXP(x, v) \ + ((x).a[LDBL_EXP_INDEX] = \ + ((x).a[LDBL_EXP_INDEX] & ~LDBL_EXP_MASK) | \ + (((IEEEl2bitsrep_part)(v) << LDBL_EXP_SHIFT) & LDBL_EXP_MASK)) + +/* Get the manl bits of x. x should be of type IEEEl2bitsrep */ +#define GET_LDOUBLE_MANL(x) \ + (((x).a[LDBL_MANL_INDEX] & LDBL_MANL_MASK) >> LDBL_MANL_SHIFT) + +/* Set the manl bit of x to v. x should be of type IEEEl2bitsrep */ +#define SET_LDOUBLE_MANL(x, v) \ + ((x).a[LDBL_MANL_INDEX] = \ + ((x).a[LDBL_MANL_INDEX] & ~LDBL_MANL_MASK) | \ + (((IEEEl2bitsrep_part)(v) << LDBL_MANL_SHIFT) & LDBL_MANL_MASK)) + +/* Get the manh bits of x. x should be of type IEEEl2bitsrep */ +#define GET_LDOUBLE_MANH(x) \ + (((x).a[LDBL_MANH_INDEX] & LDBL_MANH_MASK) >> LDBL_MANH_SHIFT) + +/* Set the manh bit of x to v. x should be of type IEEEl2bitsrep */ +#define SET_LDOUBLE_MANH(x, v) \ + ((x).a[LDBL_MANH_INDEX] = \ + ((x).a[LDBL_MANH_INDEX] & ~LDBL_MANH_MASK) | \ + (((IEEEl2bitsrep_part)(v) << LDBL_MANH_SHIFT) & LDBL_MANH_MASK)) /* * Those unions are used to convert a pointer of npy_cdouble to native C99 - * complex or our own complex type indenpendently on whether C99 complex + * complex or our own complex type independently on whether C99 complex * support is available */ #ifdef NPY_USE_C99_COMPLEX typedef union { - npy_cdouble npy_z; - complex c99_z; + npy_cdouble npy_z; + complex c99_z; } __npy_cdouble_to_c99_cast; typedef union { - npy_cfloat npy_z; - complex float c99_z; + npy_cfloat npy_z; + complex float c99_z; } __npy_cfloat_to_c99_cast; typedef union { - npy_clongdouble npy_z; - complex long double c99_z; + npy_clongdouble npy_z; + complex long double c99_z; } __npy_clongdouble_to_c99_cast; #else typedef union { - npy_cdouble npy_z; - npy_cdouble c99_z; + npy_cdouble npy_z; + npy_cdouble c99_z; } __npy_cdouble_to_c99_cast; typedef union { - npy_cfloat npy_z; - npy_cfloat c99_z; + npy_cfloat npy_z; + npy_cfloat c99_z; } __npy_cfloat_to_c99_cast; typedef union { - npy_clongdouble npy_z; - npy_clongdouble c99_z; + npy_clongdouble npy_z; + npy_clongdouble c99_z; } __npy_clongdouble_to_c99_cast; #endif |