Skip to content

Fpe fixes #25

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 2 commits into from
Closed

Fpe fixes #25

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
654f313
BUG: core: Some fixes and clean up of floating point exception code
mwiebe Dec 3, 2010
2835664
STY: Break some long lines
mwiebe Dec 3, 2010
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions numpy/core/include/numpy/npy_math.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -422,4 +422,14 @@ npy_clongdouble npy_csqrtl(npy_clongdouble z);
npy_clongdouble npy_ccosl(npy_clongdouble z);
npy_clongdouble npy_csinl(npy_clongdouble z);

/*
* Functions that set the floating point error
* status word.
*/

void npy_set_floatstatus_divbyzero(void);
void npy_set_floatstatus_overflow(void);
void npy_set_floatstatus_underflow(void);
void npy_set_floatstatus_invalid(void);

#endif
91 changes: 19 additions & 72 deletions numpy/core/include/numpy/ufuncobject.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,8 @@
#ifndef Py_UFUNCOBJECT_H
#define Py_UFUNCOBJECT_H

#include <numpy/npy_math.h>

#ifdef __cplusplus
extern "C" {
#endif
Expand Down Expand Up @@ -289,7 +292,9 @@ typedef struct _loop1d_info {
/* Solaris --------------------------------------------------------*/
/* --------ignoring SunOS ieee_flags approach, someone else can
** deal with that! */
#elif defined(sun) || defined(__BSD__) || defined(__OpenBSD__) || (defined(__FreeBSD__) && (__FreeBSD_version < 502114)) || defined(__NetBSD__)
#elif defined(sun) || defined(__BSD__) || defined(__OpenBSD__) || \
(defined(__FreeBSD__) && (__FreeBSD_version < 502114)) || \
defined(__NetBSD__)
#include <ieeefp.h>

#define UFUNC_CHECK_STATUS(ret) { \
Expand All @@ -303,9 +308,12 @@ typedef struct _loop1d_info {
(void) fpsetsticky(0); \
}

#elif defined(__GLIBC__) || defined(__APPLE__) || defined(__CYGWIN__) || defined(__MINGW32__) || (defined(__FreeBSD__) && (__FreeBSD_version >= 502114))
#elif defined(__GLIBC__) || defined(__APPLE__) || \
defined(__CYGWIN__) || defined(__MINGW32__) || \
(defined(__FreeBSD__) && (__FreeBSD_version >= 502114))

#if defined(__GLIBC__) || defined(__APPLE__) || defined(__MINGW32__) || defined(__FreeBSD__)
#if defined(__GLIBC__) || defined(__APPLE__) || \
defined(__MINGW32__) || defined(__FreeBSD__)
#include <fenv.h>
#elif defined(__CYGWIN__)
#include "fenv/fenv.c"
Expand All @@ -322,19 +330,14 @@ typedef struct _loop1d_info {
FE_UNDERFLOW | FE_INVALID); \
}

#define generate_divbyzero_error() feraiseexcept(FE_DIVBYZERO)
#define generate_overflow_error() feraiseexcept(FE_OVERFLOW)
#define generate_underflow_error() feraiseexcept(FE_UNDERFLOW)
#define generate_invalid_error() feraiseexcept(FE_INVALID)

#elif defined(_AIX)

#include <float.h>
#include <fpxcp.h>

#define UFUNC_CHECK_STATUS(ret) { \
fpflag_t fpstatus; \
\
\
fpstatus = fp_read_flag(); \
ret = ((FP_DIV_BY_ZERO & fpstatus) ? UFUNC_FPE_DIVIDEBYZERO : 0) \
| ((FP_OVERFLOW & fpstatus) ? UFUNC_FPE_OVERFLOW : 0) \
Expand All @@ -343,77 +346,21 @@ typedef struct _loop1d_info {
fp_swap_flag(0); \
}

#define generate_divbyzero_error() fp_raise_xcp(FP_DIV_BY_ZERO)
#define generate_overflow_error() fp_raise_xcp(FP_OVERFLOW)
#define generate_underflow_error() fp_raise_xcp(FP_UNDERFLOW)
#define generate_invalid_error() fp_raise_xcp(FP_INVALID)

#else

#define NO_FLOATING_POINT_SUPPORT
#define UFUNC_CHECK_STATUS(ret) { \
ret = 0; \
ret = 0; \
}

#endif

/* These should really be altered to just set the corresponding bit
in the floating point status flag. Need to figure out how to do that
on all the platforms...
*/

#if !defined(generate_divbyzero_error)
static int numeric_zero2 = 0;
static void generate_divbyzero_error(void) {
double dummy;
dummy = 1./numeric_zero2;
if (dummy) /* to prevent optimizer from eliminating expression */
return;
else /* should never be called */
numeric_zero2 += 1;
return;
}
#endif

#if !defined(generate_overflow_error)
static double numeric_two = 2.0;
static void generate_overflow_error(void) {
double dummy;
dummy = pow(numeric_two,1000);
if (dummy)
return;
else
numeric_two += 0.1;
return;
return;
}
#endif

#if !defined(generate_underflow_error)
static double numeric_small = 1e-300;
static void generate_underflow_error(void) {
double dummy;
dummy = numeric_small * 1e-300;
if (!dummy)
return;
else
numeric_small += 1e-300;
return;
}
#endif

#if !defined(generate_invalid_error)
static double numeric_inv_inf = NPY_INF;
static void generate_invalid_error(void) {
double dummy;
dummy = numeric_inv_inf - NPY_INF;
if (!dummy)
return;
else
numeric_inv_inf += 1.0;
return;
}
#endif
/*
* THESE MACROS ARE DEPRECATED.
* Use npy_set_floatstatus_* in the npymath library.
*/
#define generate_divbyzero_error() npy_set_floatstatus_dividebyzero()
#define generate_overflow_error() npy_set_floatstatus_overflow()

/* Make sure it gets defined if it isn't already */
#ifndef UFUNC_NOFPE
Expand Down
31 changes: 18 additions & 13 deletions numpy/core/src/npymath/halffloat.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,11 +77,16 @@ npy_half npy_half_spacing(npy_half h)
npy_half ret;
npy_uint16 h_exp = h&0x7c00u;
npy_uint16 h_sig = h&0x03ffu;
if (h_exp == 0x7c00u || h == 0x7bffu) {
if (h_exp == 0x7c00u) {
#if NPY_HALF_GENERATE_INVALID
generate_invalid_error();
npy_set_floatstatus_invalid();
#endif
ret = NPY_HALF_NAN;
} else if (h == 0x7bffu) {
#if NPY_HALF_GENERATE_OVERFLOW
npy_set_floatstatus_overflow();
#endif
ret = NPY_HALF_PINF;
} else if ((h&0x8000u) && h_sig == 0) { /* Negative boundary case */
if (h_exp > 0x2c00u) { /* If result is normalized */
ret = h_exp - 0x2c00u;
Expand Down Expand Up @@ -112,7 +117,7 @@ npy_half npy_half_nextafter(npy_half x, npy_half y)

if (!npy_half_isfinite(x) || npy_half_isnan(y)) {
#if NPY_HALF_GENERATE_INVALID
generate_invalid_error();
npy_set_floatstatus_invalid();
#endif
ret = NPY_HALF_NAN;
} else if (npy_half_eq_nonan(x, y)) {
Expand All @@ -134,7 +139,7 @@ npy_half npy_half_nextafter(npy_half x, npy_half y)
}
#ifdef NPY_HALF_GENERATE_OVERFLOW
if (npy_half_isinf(ret)) {
generate_overflow_error();
npy_set_floatstatus_overflow();
}
#endif

Expand Down Expand Up @@ -255,7 +260,7 @@ npy_uint16 npy_floatbits_to_halfbits(npy_uint32 f)
} else {
/* overflow to signed inf */
#if NPY_HALF_GENERATE_OVERFLOW
generate_overflow_error();
npy_set_floatstatus_overflow();
#endif
return (npy_uint16) (h_sgn + 0x7c00u);
}
Expand All @@ -271,7 +276,7 @@ npy_uint16 npy_floatbits_to_halfbits(npy_uint32 f)
#if NPY_HALF_GENERATE_UNDERFLOW
/* If f != 0, it underflowed to 0 */
if ((f&0x7fffffff) != 0) {
generate_underflow_error();
npy_set_floatstatus_underflow();
}
#endif
return h_sgn;
Expand All @@ -282,7 +287,7 @@ npy_uint16 npy_floatbits_to_halfbits(npy_uint32 f)
#if NPY_HALF_GENERATE_UNDERFLOW
/* If it's not exactly represented, it underflowed */
if ((f_sig&(((npy_uint32)1 << (126 - f_exp)) - 1)) != 0) {
generate_underflow_error();
npy_set_floatstatus_underflow();
}
#endif
f_sig >>= (113 - f_exp);
Expand Down Expand Up @@ -334,7 +339,7 @@ npy_uint16 npy_floatbits_to_halfbits(npy_uint32 f)
#if NPY_HALF_GENERATE_OVERFLOW
h_sig += h_exp;
if (h_sig == 0x7c00u) {
generate_overflow_error();
npy_set_floatstatus_overflow();
}
return h_sgn + h_sig;
#else
Expand Down Expand Up @@ -370,7 +375,7 @@ npy_uint16 npy_doublebits_to_halfbits(npy_uint64 d)
} else {
/* overflow to signed inf */
#if NPY_HALF_GENERATE_OVERFLOW
generate_overflow_error();
npy_set_floatstatus_overflow();
#endif
return h_sgn + 0x7c00u;
}
Expand All @@ -385,8 +390,8 @@ npy_uint16 npy_doublebits_to_halfbits(npy_uint64 d)
if (d_exp < 0x3e60000000000000u) {
#if NPY_HALF_GENERATE_UNDERFLOW
/* If d != 0, it underflowed to 0 */
if ((d&0x7fffffffffffffff) != 0) {
generate_underflow_error();
if ((d&0x7fffffffffffffffu) != 0) {
npy_set_floatstatus_underflow();
}
#endif
return h_sgn;
Expand All @@ -397,7 +402,7 @@ npy_uint16 npy_doublebits_to_halfbits(npy_uint64 d)
#if NPY_HALF_GENERATE_UNDERFLOW
/* If it's not exactly represented, it underflowed */
if ((d_sig&(((npy_uint64)1 << (1051 - d_exp)) - 1)) != 0) {
generate_underflow_error();
npy_set_floatstatus_underflow();
}
#endif
d_sig >>= (1009 - d_exp);
Expand Down Expand Up @@ -450,7 +455,7 @@ npy_uint16 npy_doublebits_to_halfbits(npy_uint64 d)
#if NPY_HALF_GENERATE_OVERFLOW
h_sig += h_exp;
if (h_sig == 0x7c00u) {
generate_overflow_error();
npy_set_floatstatus_overflow();
}
return h_sgn + h_sig;
#else
Expand Down
124 changes: 124 additions & 0 deletions numpy/core/src/npymath/ieee754.c.src
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -548,3 +548,127 @@ npy_longdouble npy_nextafterl(npy_longdouble x, npy_longdouble y)
return nextafterl(x, y);
}
#endif

/*
* Functions to set the floating point status word.
*/

#if defined(sun) || defined(__BSD__) || defined(__OpenBSD__) || \
(defined(__FreeBSD__) && (__FreeBSD_version < 502114)) || \
defined(__NetBSD__)
#include <ieeefp.h>

void npy_set_floatstatus_divbyzero(void)
{
fpsetsticky(FP_X_DZ);
}

void npy_set_floatstatus_overflow(void)
{
fpsetsticky(FP_X_OFL);
}

void npy_set_floatstatus_underflow(void)
{
fpsetsticky(FP_X_UFL);
}

void npy_set_floatstatus_invalid(void)
{
fpsetsticky(FP_X_INV);
}


#elif defined(__GLIBC__) || defined(__APPLE__) || \
defined(__CYGWIN__) || defined(__MINGW32__) || \
(defined(__FreeBSD__) && (__FreeBSD_version >= 502114))

# if defined(__GLIBC__) || defined(__APPLE__) || \
defined(__MINGW32__) || defined(__FreeBSD__)
# include <fenv.h>
# elif defined(__CYGWIN__)
# include "fenv/fenv.c"
# endif

void npy_set_floatstatus_divbyzero(void)
{
feraiseexcept(FE_DIVBYZERO);
}

void npy_set_floatstatus_overflow(void)
{
feraiseexcept(FE_OVERFLOW);
}

void npy_set_floatstatus_underflow(void)
{
feraiseexcept(FE_UNDERFLOW);
}

void npy_set_floatstatus_invalid(void)
{
feraiseexcept(FE_INVALID);
}

#elif defined(_AIX)
#include <float.h>
#include <fpxcp.h>

void npy_set_floatstatus_divbyzero(void)
{
fp_raise_xcp(FP_DIV_BY_ZERO);
}

void npy_set_floatstatus_overflow(void)
{
fp_raise_xcp(FP_OVERFLOW);
}

void npy_set_floatstatus_underflow(void)
{
fp_raise_xcp(FP_UNDERFLOW);
}

void npy_set_floatstatus_invalid(void)
{
fp_raise_xcp(FP_INVALID);
}

#else

/*
* By using a volatile floating point value,
* the compiler is forced to actually do the requested
* operations because of potential concurrency.
*
* We shouldn't write multiple values to a single
* global here, because that would cause
* a race condition.
*/
static volatile double _npy_floatstatus_x,
_npy_floatstatus_zero = 0.0, _npy_floatstatus_big = 1e300,
_npy_floatstatus_small = 1e-300, _npy_floatstatus_inf;

void npy_set_floatstatus_divbyzero(void)
{
_npy_floatstatus_x = 1.0 / _npy_floatstatus_zero;
}

void npy_set_floatstatus_overflow(void)
{
_npy_floatstatus_x = _npy_floatstatus_big * 1e300;
}

void npy_set_floatstatus_underflow(void)
{
_npy_floatstatus_x = _npy_floatstatus_small * 1e-300;
}

void npy_set_floatstatus_invalid(void)
{
_npy_floatstatus_inf = NPY_INFINITY;
_npy_floatstatus_x = _npy_floatstatus_inf - NPY_INFINITY;
}

#endif

Loading