>From cb8397f38a8a99884ffa85e2e4f6ec21e5d6cb87 Mon Sep 17 00:00:00 2001 From: Sergey Dmitrouk Date: Wed, 10 Sep 2014 15:09:49 +0300 Subject: [PATCH] Raise some fenv exceptions explicitly --- src/math/ilogb.c | 12 ++++++++++++ src/math/ilogbf.c | 12 ++++++++++++ src/math/j0.c | 9 +++++++-- src/math/j0f.c | 9 +++++++-- src/math/j1.c | 9 +++++++-- src/math/j1f.c | 9 +++++++-- src/math/jn.c | 5 ++++- src/math/jnf.c | 5 ++++- src/math/llrint.c | 13 ++++++++++++- src/math/llrintf.c | 13 ++++++++++++- src/math/llround.c | 13 ++++++++++++- src/math/llroundf.c | 13 ++++++++++++- src/math/llroundl.c | 13 ++++++++++++- src/math/pow.c | 33 +++++++++++++++++++++++++-------- src/math/sqrtl.c | 2 +- src/math/tgamma.c | 5 ++++- 16 files changed, 150 insertions(+), 25 deletions(-) diff --git a/src/math/ilogb.c b/src/math/ilogb.c index 64d4015..12c3aec 100644 --- a/src/math/ilogb.c +++ b/src/math/ilogb.c @@ -1,3 +1,5 @@ +#include +#include #include #include "libm.h" @@ -8,6 +10,16 @@ int ilogb(double x) uint64_t i = u.i; int e = i>>52 & 0x7ff; + if (x == 0.0) { + feraiseexcept(FE_INVALID); + return INT_MIN; + } + + if (isinf(x)) { + feraiseexcept(FE_INVALID); + return INT_MAX; + } + if (!e) { i <<= 12; if (i == 0) { diff --git a/src/math/ilogbf.c b/src/math/ilogbf.c index e23ba20..59fb2ba 100644 --- a/src/math/ilogbf.c +++ b/src/math/ilogbf.c @@ -1,4 +1,6 @@ +#include #include +#include #include "libm.h" int ilogbf(float x) @@ -8,6 +10,16 @@ int ilogbf(float x) uint32_t i = u.i; int e = i>>23 & 0xff; + if (x == 0.0) { + feraiseexcept(FE_INVALID); + return INT_MIN; + } + + if (isinf(x)) { + feraiseexcept(FE_INVALID); + return INT_MAX; + } + if (!e) { i <<= 9; if (i == 0) { diff --git a/src/math/j0.c b/src/math/j0.c index d722d94..9214d52 100644 --- a/src/math/j0.c +++ b/src/math/j0.c @@ -54,6 +54,7 @@ * 3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0. */ +#include #include "libm.h" static double pzero(double), qzero(double); @@ -164,10 +165,14 @@ double y0(double x) EXTRACT_WORDS(ix, lx, x); /* y0(nan)=nan, y0(<0)=nan, y0(0)=-inf, y0(inf)=0 */ - if ((ix<<1 | lx) == 0) + if ((ix<<1 | lx) == 0) { + feraiseexcept(FE_DIVBYZERO); return -1/0.0; - if (ix>>31) + } + if (ix>>31) { + feraiseexcept(FE_INVALID); return 0/0.0; + } if (ix >= 0x7ff00000) return 1/x; diff --git a/src/math/j0f.c b/src/math/j0f.c index 45883dc..2ac017b 100644 --- a/src/math/j0f.c +++ b/src/math/j0f.c @@ -14,6 +14,7 @@ */ #define _GNU_SOURCE +#include #include "libm.h" static float pzerof(float), qzerof(float); @@ -107,10 +108,14 @@ float y0f(float x) uint32_t ix; GET_FLOAT_WORD(ix, x); - if ((ix & 0x7fffffff) == 0) + if ((ix & 0x7fffffff) == 0) { + feraiseexcept(FE_DIVBYZERO); return -1/0.0f; - if (ix>>31) + } + if (ix>>31) { + feraiseexcept(FE_INVALID); return 0/0.0f; + } if (ix >= 0x7f800000) return 1/x; if (ix >= 0x40000000) { /* |x| >= 2.0 */ diff --git a/src/math/j1.c b/src/math/j1.c index df724d1..c83c568 100644 --- a/src/math/j1.c +++ b/src/math/j1.c @@ -54,6 +54,7 @@ * by method mentioned above. */ +#include #include "libm.h" static double pone(double), qone(double); @@ -156,10 +157,14 @@ double y1(double x) EXTRACT_WORDS(ix, lx, x); /* y1(nan)=nan, y1(<0)=nan, y1(0)=-inf, y1(inf)=0 */ - if ((ix<<1 | lx) == 0) + if ((ix<<1 | lx) == 0) { + feraiseexcept(FE_DIVBYZERO); return -1/0.0; - if (ix>>31) + } + if (ix>>31) { + feraiseexcept(FE_INVALID); return 0/0.0; + } if (ix >= 0x7ff00000) return 1/x; diff --git a/src/math/j1f.c b/src/math/j1f.c index 58875af..3950d4d 100644 --- a/src/math/j1f.c +++ b/src/math/j1f.c @@ -14,6 +14,7 @@ */ #define _GNU_SOURCE +#include #include "libm.h" static float ponef(float), qonef(float); @@ -106,10 +107,14 @@ float y1f(float x) uint32_t ix; GET_FLOAT_WORD(ix, x); - if ((ix & 0x7fffffff) == 0) + if ((ix & 0x7fffffff) == 0) { + feraiseexcept(FE_DIVBYZERO); return -1/0.0f; - if (ix>>31) + } + if (ix>>31) { + feraiseexcept(FE_INVALID); return 0/0.0f; + } if (ix >= 0x7f800000) return 1/x; if (ix >= 0x40000000) /* |x| >= 2.0 */ diff --git a/src/math/jn.c b/src/math/jn.c index 4878a54..a9fa1cc 100644 --- a/src/math/jn.c +++ b/src/math/jn.c @@ -34,6 +34,7 @@ * values of n>1. */ +#include #include "libm.h" static const double invsqrtpi = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */ @@ -224,8 +225,10 @@ double yn(int n, double x) if ((ix | (lx|-lx)>>31) > 0x7ff00000) /* nan */ return x; - if (sign && (ix|lx)!=0) /* x < 0 */ + if (sign && (ix|lx)!=0) { /* x < 0 */ + feraiseexcept(FE_INVALID); return 0/0.0; + } if (ix == 0x7ff00000) return 0.0; diff --git a/src/math/jnf.c b/src/math/jnf.c index f63c062..68b9cad 100644 --- a/src/math/jnf.c +++ b/src/math/jnf.c @@ -14,6 +14,7 @@ */ #define _GNU_SOURCE +#include #include "libm.h" float jnf(int n, float x) @@ -170,8 +171,10 @@ float ynf(int n, float x) ix &= 0x7fffffff; if (ix > 0x7f800000) /* nan */ return x; - if (sign && ix != 0) /* x < 0 */ + if (sign && ix != 0) { /* x < 0 */ + feraiseexcept(FE_INVALID); return 0/0.0f; + } if (ix == 0x7f800000) return 0.0f; diff --git a/src/math/llrint.c b/src/math/llrint.c index 4f583ae..701df94 100644 --- a/src/math/llrint.c +++ b/src/math/llrint.c @@ -1,8 +1,19 @@ +#include +#include #include /* uses LLONG_MAX > 2^53, see comments in lrint.c */ long long llrint(double x) { - return rint(x); + if (isnan(x) || isinf(x)) { + feraiseexcept(FE_INVALID); + return x; + } + + x = rint(x); + if (x > LLONG_MAX || x < LLONG_MIN) { + feraiseexcept(FE_INVALID); + } + return x; } diff --git a/src/math/llrintf.c b/src/math/llrintf.c index 96949a0..81153d2 100644 --- a/src/math/llrintf.c +++ b/src/math/llrintf.c @@ -1,8 +1,19 @@ +#include +#include #include /* uses LLONG_MAX > 2^24, see comments in lrint.c */ long long llrintf(float x) { - return rintf(x); + if (isnan(x) || isinf(x)) { + feraiseexcept(FE_INVALID); + return x; + } + + x = rintf(x); + if (x > LLONG_MAX || x < LLONG_MIN) { + feraiseexcept(FE_INVALID); + } + return x; } diff --git a/src/math/llround.c b/src/math/llround.c index 4d94787..ff48d9b 100644 --- a/src/math/llround.c +++ b/src/math/llround.c @@ -1,6 +1,17 @@ +#include +#include #include long long llround(double x) { - return round(x); + if (isnan(x) || isinf(x)) { + feraiseexcept(FE_INVALID); + return 0; + } + + x = round(x); + if (x > LLONG_MAX || x < LLONG_MIN) { + feraiseexcept(FE_INVALID); + } + return x; } diff --git a/src/math/llroundf.c b/src/math/llroundf.c index 19eb77e..b1bc9ec 100644 --- a/src/math/llroundf.c +++ b/src/math/llroundf.c @@ -1,6 +1,17 @@ +#include +#include #include long long llroundf(float x) { - return roundf(x); + if (isnan(x) || isinf(x)) { + feraiseexcept(FE_INVALID); + return 0; + } + + x = roundf(x); + if (x > LLONG_MAX || x < LLONG_MIN) { + feraiseexcept(FE_INVALID); + } + return x; } diff --git a/src/math/llroundl.c b/src/math/llroundl.c index 2c2ee5e..9a4e1a2 100644 --- a/src/math/llroundl.c +++ b/src/math/llroundl.c @@ -1,6 +1,17 @@ +#include +#include #include long long llroundl(long double x) { - return roundl(x); + if (isnan(x) || isinf(x)) { + feraiseexcept(FE_INVALID); + return 0; + } + + x = roundl(x); + if (x > LLONG_MAX || x < LLONG_MIN) { + feraiseexcept(FE_INVALID); + } + return x; } diff --git a/src/math/pow.c b/src/math/pow.c index 82f684b..f19ea84 100644 --- a/src/math/pow.c +++ b/src/math/pow.c @@ -57,6 +57,7 @@ * to produce the hexadecimal values shown. */ +#include #include "libm.h" static const double @@ -196,16 +197,24 @@ double pow(double x, double y) /* |y| is huge */ if (iy > 0x41e00000) { /* if |y| > 2**31 */ if (iy > 0x43f00000) { /* if |y| > 2**64, must o/uflow */ - if (ix <= 0x3fefffff) + if (ix <= 0x3fefffff) { + feraiseexcept(FE_INEXACT| ((hy < 0) ? FE_OVERFLOW : FE_UNDERFLOW)); return hy < 0 ? huge*huge : tiny*tiny; - if (ix >= 0x3ff00000) + } + if (ix >= 0x3ff00000) { + feraiseexcept(FE_INEXACT| ((hy > 0) ? FE_OVERFLOW : FE_UNDERFLOW)); return hy > 0 ? huge*huge : tiny*tiny; + } } /* over/underflow if x is not close to one */ - if (ix < 0x3fefffff) + if (ix < 0x3fefffff) { + feraiseexcept(FE_INEXACT| ((hy < 0) ? FE_OVERFLOW : FE_UNDERFLOW)); return hy < 0 ? s*huge*huge : s*tiny*tiny; - if (ix > 0x3ff00000) + } + if (ix > 0x3ff00000) { + feraiseexcept(FE_INEXACT| ((hy > 0) ? FE_OVERFLOW : FE_UNDERFLOW)); return hy > 0 ? s*huge*huge : s*tiny*tiny; + } /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ t = ax - 1.0; /* t has 20 trailing zeros */ @@ -282,15 +291,23 @@ double pow(double x, double y) z = p_l + p_h; EXTRACT_WORDS(j, i, z); if (j >= 0x40900000) { /* z >= 1024 */ - if (((j-0x40900000)|i) != 0) /* if z > 1024 */ + if (((j-0x40900000)|i) != 0) { /* if z > 1024 */ + feraiseexcept(FE_INEXACT|FE_OVERFLOW); return s*huge*huge; /* overflow */ - if (p_l + ovt > z - p_h) + } + if (p_l + ovt > z - p_h) { + feraiseexcept(FE_INEXACT|FE_OVERFLOW); return s*huge*huge; /* overflow */ + } } else if ((j&0x7fffffff) >= 0x4090cc00) { /* z <= -1075 */ // FIXME: instead of abs(j) use unsigned j - if (((j-0xc090cc00)|i) != 0) /* z < -1075 */ + if (((j-0xc090cc00)|i) != 0) { /* z < -1075 */ + feraiseexcept(FE_INEXACT|FE_UNDERFLOW); return s*tiny*tiny; /* underflow */ - if (p_l <= z - p_h) + } + if (p_l <= z - p_h) { + feraiseexcept(FE_INEXACT|FE_UNDERFLOW); return s*tiny*tiny; /* underflow */ + } } /* * compute 2**(p_h+p_l) diff --git a/src/math/sqrtl.c b/src/math/sqrtl.c index 83a8f80..0872e15 100644 --- a/src/math/sqrtl.c +++ b/src/math/sqrtl.c @@ -3,5 +3,5 @@ long double sqrtl(long double x) { /* FIXME: implement in C, this is for LDBL_MANT_DIG == 64 only */ - return sqrt(x); + return isnan(x) ? 0.0l/0.0l : sqrt(x); } diff --git a/src/math/tgamma.c b/src/math/tgamma.c index 28f6e0f..51e7e20 100644 --- a/src/math/tgamma.c +++ b/src/math/tgamma.c @@ -22,6 +22,7 @@ Gamma(x)*Gamma(-x) = -pi/(x sin(pi x)) most ideas and constants are from boost and python */ +#include #include "libm.h" static const double pi = 3.141592653589793238462643383279502884; @@ -124,8 +125,10 @@ double tgamma(double x) /* integer arguments */ /* raise inexact when non-integer */ if (x == floor(x)) { - if (sign) + if (sign) { + feraiseexcept(FE_INVALID); return 0/0.0; + } if (x <= sizeof fact/sizeof *fact) return fact[(int)x - 1]; } -- 1.8.4