diff options
49 files changed, 581 insertions, 14 deletions
diff --git a/src/math/acos.c b/src/math/acos.c index 2f740d03..9a356100 100644 --- a/src/math/acos.c +++ b/src/math/acos.c @@ -1,12 +1,15 @@ # define TGSOURCE "acos.c" #include <math.h> +#include "fenv.h" #include "errno.h" #include "_tgmath.h" /** arc cosine **/ TYPE TGFN(acos)(TYPE x) { - if (x < -1 || x > 1) { + if (TGFN(fabs)(x) > 1) { + feraiseexcept(FE_INVALID); + return NAN; errno = EDOM; /* ARGUMENT(x) not in the range [-1, +1] */ return TGHUGE; } @@ -17,6 +20,10 @@ TYPE TGFN(acos)(TYPE x) return TGHUGE; } + if (x == 1.0) { + return 0.0; + } + /* RETURN_SUCCESS(a value in range `[0, PI()]'); */ return (TYPE)0.0; } diff --git a/src/math/acosh.c b/src/math/acosh.c index aeb15816..f557c875 100644 --- a/src/math/acosh.c +++ b/src/math/acosh.c @@ -1,9 +1,23 @@ # define TGSOURCE "acosh.c" #include "_tgmath.h" #include <math.h> +#include "fenv.h" TYPE TGFN(acosh)(TYPE x) { + if (x == 1.0) { + return 0.0; + } + + if (x < 1.0) { + feraiseexcept(FE_INVALID); + return NAN; + } + + if (fpclassify(x) == FP_INFINITE && !signbit(x)) { + return x; + } + return x; } diff --git a/src/math/asin.c b/src/math/asin.c index 92dae8c8..a94338f5 100644 --- a/src/math/asin.c +++ b/src/math/asin.c @@ -6,11 +6,15 @@ /** arc sine **/ TYPE TGFN(asin)(TYPE x) { - if (x < -1.0 || x > 1.0) { + if (TGFN(fabs)(x) > 1.0) { errno = EDOM; /* ARGUMENT(x) is not in the range [-1, +1] */ return TGHUGE; } + if (fpclassify(x) == FP_ZERO) { + return x; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/asinh.c b/src/math/asinh.c index 421f574d..9f4589e6 100644 --- a/src/math/asinh.c +++ b/src/math/asinh.c @@ -4,6 +4,12 @@ TYPE TGFN(asinh)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + default: break; + } + return x; } diff --git a/src/math/atan.c b/src/math/atan.c index 141367cb..6d9eee63 100644 --- a/src/math/atan.c +++ b/src/math/atan.c @@ -3,6 +3,8 @@ #include "errno.h" #include "_tgmath.h" +#include "M_PI_2.c" + /** arc tangent **/ TYPE TGFN(atan)(TYPE x) { @@ -11,6 +13,12 @@ TYPE TGFN(atan)(TYPE x) TYPE power = 1.0; int i; + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return signbit(x) ? - M_PI_2 : M_PI_2; + default: break; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/atan2.c b/src/math/atan2.c index c4c56a63..3b88b617 100644 --- a/src/math/atan2.c +++ b/src/math/atan2.c @@ -4,11 +4,60 @@ #include "errno.h" #include "nonstd/assert.h" +#include "M_PI.c" +#include "M_PI_2.c" + /** arc tangent **/ TYPE TGFN(atan2)(TYPE y, TYPE x) { + int classy = fpclassify(y); + int classx = fpclassify(x); ASSERT_NONZERO(x); + if (classy == FP_ZERO && classx == 0) { + if (signbit(x)) { + return copysign(M_PI, y); + } else { + return y; + } + } + + if (classy == FP_ZERO) { + if (x < 0) { + return copysign(M_PI, y); + } else if (x > 0) { + return y; + } + } + + if (classx == FP_ZERO) { + if (y < 0) { + return - M_PI_2; + } else if (y > 0) { + return M_PI_2; + } + } + + if (classx == FP_INFINITE) { + if (classy == FP_INFINITE) { + if (signbit(x)) { + return copysign(3.0 * M_PI / 4.0, y); + } else { + return copysign(M_PI / 4.0, y); + } + } else if (y > 0) { + if (signbit(x)) { + return copysign(M_PI, y); + } else { + return copysign(0.0, y); + } + } + } + + if (classy == FP_INFINITE && classx != FP_INFINITE) { + return copysign(M_PI_2, y); + } + if (y == 0 && x == 0) { errno = EDOM; /* ARGUMENT(y) and ARGUMENT(x) are both LITERAL(0)) */ return TGHUGE; diff --git a/src/math/atanh.c b/src/math/atanh.c index eef1cfc6..bfd4e57b 100644 --- a/src/math/atanh.c +++ b/src/math/atanh.c @@ -1,9 +1,24 @@ # define TGSOURCE "atanh.c" #include "_tgmath.h" #include <math.h> +#include "fenv.h" TYPE TGFN(atanh)(TYPE x) { + if (fpclassify(x) == FP_ZERO) { + return x; + } + + if (TGFN(fabs)(x) == 1.0) { + feraiseexcept(FE_DIVBYZERO); + return copysign(INFINITY, x); + } + + if (TGFN(fabs)(x) > 1.0) { + feraiseexcept(FE_INVALID); + return NAN; + } + return x; } diff --git a/src/math/cbrt.c b/src/math/cbrt.c index 156093ba..74c48d52 100644 --- a/src/math/cbrt.c +++ b/src/math/cbrt.c @@ -4,6 +4,11 @@ TYPE TGFN(cbrt)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + } + return x; } diff --git a/src/math/ceil.c b/src/math/ceil.c index 7af395ac..d90ba15f 100644 --- a/src/math/ceil.c +++ b/src/math/ceil.c @@ -6,6 +6,12 @@ /** round up to nearest integer **/ TYPE TGFN(ceil)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + default: break; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/copysign.c b/src/math/copysign.c index a462348b..2ace3c1d 100644 --- a/src/math/copysign.c +++ b/src/math/copysign.c @@ -4,7 +4,16 @@ TYPE TGFN(copysign)(TYPE x, TYPE y) { - return x - y; + if (isnan(x)) { + if (signbit(y)) { + /* return -NaN; */ + } else { + /* return NaN; */ + } + } + + x = TGFN(fabs)(x); + return signbit(y) ? -x : x; } /* diff --git a/src/math/cos.c b/src/math/cos.c index d10629a3..86a1c995 100644 --- a/src/math/cos.c +++ b/src/math/cos.c @@ -1,6 +1,7 @@ # define TGSOURCE "cos.c" #include <math.h> #include "_tgmath.h" +#include "fenv.h" #include "errno.h" /** cosine **/ @@ -12,6 +13,12 @@ TYPE TGFN(cos)(TYPE x) TYPE cosine = 1.0; TYPE power = 1.0; + switch (fpclassify(x)) { + case FP_ZERO: return 1.0; + case FP_INFINITE: feraiseexcept(FE_INVALID); return NAN; + default: break; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/cosh.c b/src/math/cosh.c index 9065cd03..0f981e36 100644 --- a/src/math/cosh.c +++ b/src/math/cosh.c @@ -6,6 +6,12 @@ /** hyperbolic cosine **/ TYPE TGFN(cosh)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return 1.0; + case FP_INFINITE: return TGFN(fabs)(x); + default: break; + } + if (0) { errno = ERANGE; /* The magnitude of ARGUMENT(x) is too large */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/erf.c b/src/math/erf.c index 33059480..2115982b 100644 --- a/src/math/erf.c +++ b/src/math/erf.c @@ -4,6 +4,12 @@ TYPE TGFN(erf)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return signbit(x) ? -1.0 : 1.0; + default: break; + } + return x; } diff --git a/src/math/erfc.c b/src/math/erfc.c index 84579e10..b937590d 100644 --- a/src/math/erfc.c +++ b/src/math/erfc.c @@ -4,6 +4,10 @@ TYPE TGFN(erfc)(TYPE x) { + if (fpclassify(x) == FP_INFINITE) { + return signbit(x) ? 2.0 : 0.0; + } + return x; } diff --git a/src/math/exp.c b/src/math/exp.c index f244a60f..412697b2 100644 --- a/src/math/exp.c +++ b/src/math/exp.c @@ -12,6 +12,12 @@ TYPE TGFN(exp)(TYPE x) TYPE exponent = 1.0; TYPE power = 1.0; + switch (fpclassify(x)) { + case FP_ZERO: return 1.0; + case FP_INFINITE: return signbit(x) ? 0.0 : x; + default: break; + } + if (0) { errno = ERANGE; /* The magnitude of ARGUMENT(x) is too large */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/exp2.c b/src/math/exp2.c index 599634e6..16ff5940 100644 --- a/src/math/exp2.c +++ b/src/math/exp2.c @@ -4,6 +4,12 @@ TYPE TGFN(exp2)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return 1.0; + case FP_INFINITE: return signbit(x) ? 0.0 : x; + default: break; + } + return x; } diff --git a/src/math/expm1.c b/src/math/expm1.c index 152c6c1e..40684b05 100644 --- a/src/math/expm1.c +++ b/src/math/expm1.c @@ -4,6 +4,12 @@ TYPE TGFN(expm1)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return signbit(x) ? -1.0 : x; + default: break; + } + return x; } diff --git a/src/math/fabs.c b/src/math/fabs.c index 9664b702..5fc42750 100644 --- a/src/math/fabs.c +++ b/src/math/fabs.c @@ -6,6 +6,12 @@ /** absolute value **/ TYPE TGFN(fabs)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return 0.0; + case FP_INFINITE: return INFINITY; + default: break; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ @@ -13,7 +19,7 @@ TYPE TGFN(fabs)(TYPE x) } /* RETURN_SUCCESS(ABS(ARGUMENT(x))); */ - return x < 0.0 ? -x : x; + return signbit(x) ? -x : x; } /*** diff --git a/src/math/fdim.c b/src/math/fdim.c index c39c6753..4f6b3896 100644 --- a/src/math/fdim.c +++ b/src/math/fdim.c @@ -4,7 +4,11 @@ TYPE TGFN(fdim)(TYPE x, TYPE y) { - return x - y; + if (x > y) { + return x - y; + } + + return 0; } /* diff --git a/src/math/floor.c b/src/math/floor.c index baa7d3a6..41350fe8 100644 --- a/src/math/floor.c +++ b/src/math/floor.c @@ -6,6 +6,12 @@ /** round down to nearest integer **/ TYPE TGFN(floor)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + default: break; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/fma.c b/src/math/fma.c index 18b059d1..845afea4 100644 --- a/src/math/fma.c +++ b/src/math/fma.c @@ -1,10 +1,32 @@ # define TGSOURCE "fma.c" #include "_tgmath.h" #include <math.h> +#include "fenv.h" TYPE TGFN(fma)(TYPE x, TYPE y, TYPE z) { - return x - y - z; + int classx = fpclassify(x); + int classy = fpclassify(y); + /* int classz = fpclassify(z); */ + + if (classx == FP_INFINITE && classy == FP_ZERO) { + feraiseexcept(FE_INVALID); + return NAN; + } + + if (classx == FP_ZERO && classy == FP_INFINITE) { + feraiseexcept(FE_INVALID); + return NAN; + } + + /* + if (x * y == -z && fpclassify(x*y) == classz == FP_INFINITE) { + feraiseexcept(FE_INVALID); + return NAN; + } + */ + + return x * y + z; } /* diff --git a/src/math/fmax.c b/src/math/fmax.c index dc354ca3..7e5d8a61 100644 --- a/src/math/fmax.c +++ b/src/math/fmax.c @@ -4,7 +4,13 @@ TYPE TGFN(fmax)(TYPE x, TYPE y) { - return x - y; + if (isnan(x) && !isnan(y)) { + return y; + } else if (!isnan(x) && isnan(y)) { + return x; + } + + return x > y ? x : y; } /* diff --git a/src/math/fmin.c b/src/math/fmin.c index 822f893b..4408f432 100644 --- a/src/math/fmin.c +++ b/src/math/fmin.c @@ -4,7 +4,13 @@ TYPE TGFN(fmin)(TYPE x, TYPE y) { - return x - y; + if (isnan(x) && !isnan(y)) { + return y; + } else if (!isnan(x) && isnan(y)) { + return x; + } + + return x < y ? x : y; } /* diff --git a/src/math/fmod.c b/src/math/fmod.c index d8254d40..4370830f 100644 --- a/src/math/fmod.c +++ b/src/math/fmod.c @@ -2,10 +2,27 @@ #include <math.h> #include "_tgmath.h" #include "errno.h" +#include "fenv.h" /** floating-point remainder **/ TYPE TGFN(fmod)(TYPE x, TYPE y) { + int classx = fpclassify(x); + int classy = fpclassify(y); + + if (classx == FP_ZERO && classy != FP_ZERO) { + return x; + } + + if (classx == FP_INFINITE && classy == FP_ZERO) { + feraiseexcept(FE_INVALID); + return NAN; + } + + if (classx != FP_INFINITE && classy == FP_INFINITE) { + return x; + } + if (y == 0) { return 0; } diff --git a/src/math/frexp.c b/src/math/frexp.c index 40bb59c8..b60dadc5 100644 --- a/src/math/frexp.c +++ b/src/math/frexp.c @@ -6,6 +6,13 @@ /** extract mantissa and exponent **/ TYPE TGFN(frexp)(TYPE value, int *exp) { + switch (fpclassify(value)) { + case FP_ZERO: *exp = 0; return value; + case FP_INFINITE: *exp = /* unspecified */0; return value; + case FP_NAN: *exp = /* unspecified */0; return value; + default: break; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/hypot.c b/src/math/hypot.c index 6ab2ecf7..ce902c59 100644 --- a/src/math/hypot.c +++ b/src/math/hypot.c @@ -4,6 +4,14 @@ TYPE TGFN(hypot)(TYPE x, TYPE y) { + if (fpclassify(y) == FP_ZERO) { + return fabs(x); + } + + if (fpclassify(x) == FP_INFINITE) { + return INFINITY; + } + return x - y; } diff --git a/src/math/ilogb.c b/src/math/ilogb.c index 0d826e04..d1295815 100644 --- a/src/math/ilogb.c +++ b/src/math/ilogb.c @@ -1,9 +1,18 @@ # define TGSOURCE "ilogb.c" #include "_tgmath.h" #include <math.h> +#include "limits.h" +#include "fenv.h" -TYPE TGFN(ilogb)(TYPE x) +int TGFN(ilogb)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: feraiseexcept(FE_INVALID); return FP_ILOGB0; + case FP_INFINITE: feraiseexcept(FE_INVALID); return INT_MAX; + case FP_NAN: feraiseexcept(FE_INVALID); return FP_ILOGBNAN; + default: break; + } + return x; } diff --git a/src/math/ldexp.c b/src/math/ldexp.c index ab143982..33e2f381 100644 --- a/src/math/ldexp.c +++ b/src/math/ldexp.c @@ -6,6 +6,16 @@ /** multiply by a power of 2 **/ TYPE TGFN(ldexp)(TYPE x, int exp) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + default: break; + } + + if (exp == 0) { + return x; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/lgamma.c b/src/math/lgamma.c index 3612c825..3f0799c2 100644 --- a/src/math/lgamma.c +++ b/src/math/lgamma.c @@ -1,9 +1,25 @@ # define TGSOURCE "lgamma.c" #include "_tgmath.h" #include <math.h> +#include "fenv.h" TYPE TGFN(lgamma)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: feraiseexcept(FE_DIVBYZERO); return INFINITY; + case FP_INFINITE: return INFINITY; + default: break; + } + + if (x == 1.0 || x == 2.0) { + return 0.0; + } + + if (signbit(x)) { + feraiseexcept(FE_DIVBYZERO); + return INFINITY; + } + return x; } diff --git a/src/math/llrint.c b/src/math/llrint.c index a967f651..299a8258 100644 --- a/src/math/llrint.c +++ b/src/math/llrint.c @@ -1,9 +1,23 @@ # define TGSOURCE "llrint.c" #include "_tgmath.h" #include <math.h> +#include "limits.h" +#include "fenv.h" long long int TGFN(llrint)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: + return 0; + + case FP_INFINITE: + feraiseexcept(FE_INVALID); + return signbit(x) ? LLONG_MIN : LLONG_MAX; + + default: + break; + } + return x; } diff --git a/src/math/log.c b/src/math/log.c index ee809cf5..abc29e4e 100644 --- a/src/math/log.c +++ b/src/math/log.c @@ -2,11 +2,24 @@ #include <math.h> #include "_tgmath.h" #include "errno.h" +#include "fenv.h" /** natural logarithm **/ TYPE TGFN(log)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: feraiseexcept(FE_DIVBYZERO); return - INFINITY; + case FP_INFINITE: if (!signbit(x)) { return x; } + default: break; + } + + if (x == 1.0) { + return 0.0; + } + if (x < 0) { + feraiseexcept(FE_INVALID); + return NAN; errno = EDOM; /* ARGUMENT(x) is negative */ return TGHUGE; } diff --git a/src/math/log10.c b/src/math/log10.c index cea1d45b..e504df44 100644 --- a/src/math/log10.c +++ b/src/math/log10.c @@ -2,11 +2,24 @@ #include <math.h> #include "_tgmath.h" #include "errno.h" +#include "fenv.h" /** base-10 logarithm **/ TYPE TGFN(log10)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: feraiseexcept(FE_INVALID); return - INFINITY; + case FP_INFINITE: if (!signbit(x)) { return x; } break; + default: break; + } + + if (x == 1.0) { + return 0.0; + } + if (x < 0) { + feraiseexcept(FE_INVALID); + return NAN; errno = EDOM; /* ARGUMENT(x) is negative */ return TGHUGE; } diff --git a/src/math/log1p.c b/src/math/log1p.c index 30e4a8d1..5f3e2684 100644 --- a/src/math/log1p.c +++ b/src/math/log1p.c @@ -1,9 +1,26 @@ # define TGSOURCE "log1p.c" #include "_tgmath.h" #include <math.h> +#include "fenv.h" TYPE TGFN(log1p)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: if (!signbit(x)) { return x; } break; + default: break; + } + + if (x == -1) { + feraiseexcept(FE_DIVBYZERO); + return - INFINITY; + } + + if (x < -1) { + feraiseexcept(FE_INVALID); + return NAN; + } + return x; } diff --git a/src/math/log2.c b/src/math/log2.c index c438d7ce..55569065 100644 --- a/src/math/log2.c +++ b/src/math/log2.c @@ -1,9 +1,25 @@ #define TGSOURCE "log2.c" #include "_tgmath.h" #include <math.h> +#include "fenv.h" TYPE TGFN(log2)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: feraiseexcept(FE_DIVBYZERO); return - INFINITY; + case FP_INFINITE: if (!signbit(x)) { return x; } break; + default: break; + } + + if (x == 1.0) { + return 0.0; + } + + if (x < 0.0) { + feraiseexcept(FE_INVALID); + return NAN; + } + return x; } diff --git a/src/math/logb.c b/src/math/logb.c index d0bbfaca..8dd29371 100644 --- a/src/math/logb.c +++ b/src/math/logb.c @@ -1,9 +1,16 @@ # define TGSOURCE "logb.c" #include "_tgmath.h" #include <math.h> +#include "fenv.h" TYPE TGFN(logb)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: feraiseexcept(FE_DIVBYZERO); return - INFINITY; + case FP_INFINITE: return x; + default: break; + } + return x; } diff --git a/src/math/lrint.c b/src/math/lrint.c index 16eaf7be..3ce1dd66 100644 --- a/src/math/lrint.c +++ b/src/math/lrint.c @@ -1,9 +1,23 @@ # define TGSOURCE "lrint.c" #include "_tgmath.h" #include <math.h> +#include "limits.h" +#include "fenv.h" long int TGFN(lrint)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: + return 0; + + case FP_INFINITE: + feraiseexcept(FE_INVALID); + return signbit(x) ? LONG_MIN : LONG_MAX; + + default: + break; + } + return x; } diff --git a/src/math/modf.c b/src/math/modf.c index fc09fc2d..63db6094 100644 --- a/src/math/modf.c +++ b/src/math/modf.c @@ -6,6 +6,12 @@ /** decompose floating-point numbers **/ TYPE TGFN(modf)(TYPE value, TYPE *iptr) { + switch (fpclassify(value)) { + case FP_INFINITE: *iptr = value; return copysign(0.0, value); + case FP_NAN: *iptr = value; return value; + default: break; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ @@ -13,7 +19,6 @@ TYPE TGFN(modf)(TYPE value, TYPE *iptr) } /* RETURN_SUCCESS(the signed fractional part of ARGUMENT(value)); */ - (void)iptr; return value; } diff --git a/src/math/nan.c b/src/math/nan.c index bd0cbbc4..c0874e5f 100644 --- a/src/math/nan.c +++ b/src/math/nan.c @@ -1,11 +1,22 @@ # define TGSOURCE "nan.c" #include "_tgmath.h" #include <math.h> +#include "string.h" +#include "stdlib.h" + +#define strtodf strtof +#define strtodl strtold TYPE TGFN(nan)(const char *tagp) { - (void)tagp; - return 0.0; + if (tagp) { + char ncharseq[strlen(tagp) + 6]; + strcpy(ncharseq, "NAN("); + strcat(ncharseq, tagp); + strcat(ncharseq, ")"); + return TGFN(strtod)(ncharseq, (char**)NULL); + } + return TGFN(strtod)("NAN", (char**)NULL); } /* diff --git a/src/math/nearbyint.c b/src/math/nearbyint.c index d5d2cc21..ac9417db 100644 --- a/src/math/nearbyint.c +++ b/src/math/nearbyint.c @@ -4,6 +4,12 @@ TYPE TGFN(nearbyint)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + default: break; + } + return x; } diff --git a/src/math/pow.c b/src/math/pow.c index c0d613e5..8a070005 100644 --- a/src/math/pow.c +++ b/src/math/pow.c @@ -2,10 +2,91 @@ #include <math.h> #include "_tgmath.h" #include "errno.h" +#include "fenv.h" /** exponentiation **/ TYPE TGFN(pow)(TYPE x, TYPE y) { + int classx = fpclassify(x); + int classy = fpclassify(y); + + if (classx == FP_ZERO) { + if (classy == FP_INFINITE && signbit(y)) { + feraiseexcept(FE_DIVBYZERO); + return INFINITY; + } + + if (y < 0) { + feraiseexcept(FE_DIVBYZERO); + /* if (y == odd integer) { */ + return signbit(x) ? - INFINITY : INFINITY; + /* } else { */ + return INFINITY; + /* } */ + } + + if (y > 0) { + /* if (y == odd integer) { */ + return x; + /* } else { */ + return 0.0; + /* } */ + } + } + + if (x == -1.0 && classy == FP_INFINITE) { + return 1; + } + + if (x == 1.0) { + return 1.0; + } + + if (classy == FP_ZERO) { + return 1.0; + } + + if (x < 0 /* && y != integer */) { + feraiseexcept(FE_INVALID); + return NAN; + } + + if (classy == FP_INFINITE) { + int sb = signbit(y); + TYPE ab = TGFN(fabs)(x); + if ((sb && ab < 1.0) || (!sb && ab > 1.0)) { + return INFINITY; + } else { + return 0.0; + } + } + + if (classx == FP_INFINITE) { + if (!signbit(x)) { + if (y < 0) { + return 0.0; + } else if (y > 0) { + return INFINITY; + } + } + + if (y < 0) { + /* if (y == odd integer ) { */ + return -0.0; + /* } else { */ + return 0.0; + /* } */ + } + + if (y > 0) { + /* if (y == odd integer ) { */ + return - INFINITY; + /* } else { */ + return INFINITY; + /* } */ + } + } + if (x < 0 /* && !isintegral(y) */) { errno = EDOM; /* ARGUMENT(x) is negative and ARGUMENT(y) is not an integer */ return TGHUGE; diff --git a/src/math/rint.c b/src/math/rint.c index a1868cbb..a97c4709 100644 --- a/src/math/rint.c +++ b/src/math/rint.c @@ -4,6 +4,12 @@ TYPE TGFN(rint)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + default: break; + } + return x; } diff --git a/src/math/round.c b/src/math/round.c index c154bdbf..5eedd578 100644 --- a/src/math/round.c +++ b/src/math/round.c @@ -1,10 +1,33 @@ # define TGSOURCE "round.c" #include "_tgmath.h" #include <math.h> +#include "fenv.h" + +#ifndef __GNUC__ +#pragma STDC FENV_ACCESS ON +#endif TYPE TGFN(round)(TYPE x) { - return x; + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + default: break; + } + + fenv_t save_env; + feholdexcept(&save_env); + + TYPE ret = TGFN(rint)(x); + + if (fetestexcept(FE_INEXACT)) { + fesetround(FE_TOWARDZERO); + ret = TGFN(rint)(TGFN(copysign)(0.5 + TGFN(fabs)(x), x)); + } + + feupdateenv(&save_env); + + return ret; } /* diff --git a/src/math/scalbn.c b/src/math/scalbn.c index bf5b4aae..f01843ce 100644 --- a/src/math/scalbn.c +++ b/src/math/scalbn.c @@ -4,7 +4,17 @@ TYPE TGFN(scalbn)(TYPE x, int n) { - return x - n; + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + default: break; + } + + if (n == 0) { + return x; + } + + return TGFN(ldexp)(x, n); } /* diff --git a/src/math/sin.c b/src/math/sin.c index b753cecb..e93a663b 100644 --- a/src/math/sin.c +++ b/src/math/sin.c @@ -2,6 +2,7 @@ #include <math.h> #include "_tgmath.h" #include "errno.h" +#include "fenv.h" /** sine **/ TYPE TGFN(sin)(TYPE x) @@ -12,6 +13,12 @@ TYPE TGFN(sin)(TYPE x) TYPE sine = 0.0; TYPE power = 1.0; + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: feraiseexcept(FE_INVALID); return NAN; + default: break; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/sinh.c b/src/math/sinh.c index e5942de6..2e8d498c 100644 --- a/src/math/sinh.c +++ b/src/math/sinh.c @@ -6,6 +6,12 @@ /** hyperbolic sine **/ TYPE TGFN(sinh)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + default: break; + } + if (0) { errno = ERANGE; /* The magnitude of ARGUMENT(x) is too large */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/tan.c b/src/math/tan.c index 2e60bd2f..061e6cdd 100644 --- a/src/math/tan.c +++ b/src/math/tan.c @@ -2,10 +2,17 @@ #include <math.h> #include "_tgmath.h" #include "errno.h" +#include "fenv.h" /** tangent **/ TYPE TGFN(tan)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: feraiseexcept(FE_INVALID); return NAN; + default: break; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/tanh.c b/src/math/tanh.c index d64be283..cf5bf028 100644 --- a/src/math/tanh.c +++ b/src/math/tanh.c @@ -6,6 +6,12 @@ /** hyperbolic tangent **/ TYPE TGFN(tanh)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return signbit(x) ? -1.0 : 1.0; + default: break; + } + if (0) { errno = ERANGE; /* The result cannot be represented */ /* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */ diff --git a/src/math/tgamma.c b/src/math/tgamma.c index b80584ac..63e25a5b 100644 --- a/src/math/tgamma.c +++ b/src/math/tgamma.c @@ -1,9 +1,27 @@ # define TGSOURCE "tgamma.c" #include "_tgmath.h" #include <math.h> +#include "fenv.h" TYPE TGFN(tgamma)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: + feraiseexcept(FE_DIVBYZERO); + return signbit(x) ? - INFINITY : INFINITY; + + case FP_INFINITE: + if (signbit(x)) { + feraiseexcept(FE_INVALID); + return NAN; + } else { + return INFINITY; + } + + default: + break; + } + return x; } diff --git a/src/math/trunc.c b/src/math/trunc.c index 4b1081a6..4ee27109 100644 --- a/src/math/trunc.c +++ b/src/math/trunc.c @@ -4,6 +4,12 @@ TYPE TGFN(trunc)(TYPE x) { + switch (fpclassify(x)) { + case FP_ZERO: return x; + case FP_INFINITE: return x; + default: break; + } + return x; } |