From aee847b4182c38b77f06a4059b6a2646152d4391 Mon Sep 17 00:00:00 2001 From: Adam Debek Date: Fri, 10 Jan 2025 15:45:45 +0100 Subject: [PATCH] math: improve POSIX compliance JIRA: RTOS-974 --- include/math.h | 4 ++++ math/common.c | 8 ++++++++ math/exp.c | 37 +++++++++++++++++++++++++++++++++++-- math/hyper.c | 32 ++++++++++++++++++++++++++++++-- math/power.c | 26 +++++++++++++++++--------- math/trig.c | 24 +++++++++++++++++------- 6 files changed, 111 insertions(+), 20 deletions(-) diff --git a/include/math.h b/include/math.h index 85d50ba6..c80ea4ad 100644 --- a/include/math.h +++ b/include/math.h @@ -25,6 +25,10 @@ extern "C" { #endif +#define MATH_ERRNO 1 +#define MATH_ERREXCEPT 2 +#define math_errhandling (MATH_ERRNO) + #define FP_NAN 0 #define FP_INFINITE 1 #define FP_ZERO 2 diff --git a/math/common.c b/math/common.c index 5600536c..97c772ed 100644 --- a/math/common.c +++ b/math/common.c @@ -14,6 +14,9 @@ */ #include +#include +#include +#include #include "common.h" @@ -133,6 +136,11 @@ double quickPow(double x, int e) } } + /* ia32 use 80-bit extended precision registers so isinf() may not be true */ + if ((isinf(res) != 0) || (fabs(res) > DBL_MAX)) { + errno = ERANGE; + } + return res; } diff --git a/math/exp.c b/math/exp.c index 7cc39ccb..f9f0c5d4 100644 --- a/math/exp.c +++ b/math/exp.c @@ -23,6 +23,14 @@ double frexp(double x, int *exp) { conv_t *conv = (conv_t *)&x; + if (isnan(x) != 0) { + return NAN; + } + + if (isinf(x) != 0) { + return x; + } + *exp = 0; if ((x == 0.0) || (x == -0.0)) { @@ -45,6 +53,10 @@ double ldexp(double x, int exp) conv_t *conv = (conv_t *)&x; int exponent = 0; + if (isnan(x) != 0) { + return NAN; + } + if ((x == 0.0) || (x == -0.0)) { return x; } @@ -79,7 +91,10 @@ double log(double x) conv_t *conv = (conv_t *)&tmp; int exp = 0, i; - if (x < 0.0) { + if (isnan(x) != 0) { + return NAN; + } + else if (x < 0.0) { errno = EDOM; return NAN; } @@ -132,6 +147,10 @@ double modf(double x, double *intpart) int exp = conv->i.exponent - 1023; uint64_t m, mask = 0xfffffffffffffLL; + if (isnan(x) != 0) { + return NAN; + } + if (exp > 52) { *intpart = x; return (conv->i.sign ? -0.0 : 0.0); @@ -176,7 +195,13 @@ double exp(double x) double res, resi, powx, e, factorial; int i; - if (x > 710.0) { + if (isnan(x) != 0) { + return NAN; + } + + /* Values of x greater than 709.79 will cause overflow, returning INFINITY */ + if (x > 709.79) { + errno = ERANGE; return HUGE_VAL; } @@ -208,6 +233,10 @@ double ceil(double x) { double ipart, fpart; + if (isnan(x) != 0) { + return NAN; + } + fpart = modf(x, &ipart); if ((x > 0.0) && ((fpart + x) != x)) { @@ -222,6 +251,10 @@ double floor(double x) { double ipart, fpart; + if (isnan(x) != 0) { + return NAN; + } + fpart = modf(x, &ipart); if ((x < 0.0) && ((fpart + x) != x)) { diff --git a/math/hyper.c b/math/hyper.c index dc4d410f..c78a7ea8 100644 --- a/math/hyper.c +++ b/math/hyper.c @@ -14,25 +14,53 @@ */ #include +#include +#include double cosh(double x) { + double y; + + if (isnan(x) != 0) { + return NAN; + } + if ((x == INFINITY) || (x == -INFINITY)) { return INFINITY; } - return ((exp(x) + exp(-x)) / 2.0); + /* Make sure exp(x) is not infinity */ + if (x < 709.78) { + return ((exp(x) + exp(-x)) / 2.0); + } + else { + y = cosh(x / 2.0); + return ((2.0 * y * y) - 1.0); + } } double sinh(double x) { + double y; + + if (isnan(x) != 0) { + return NAN; + } + if ((x == INFINITY) || (x == -INFINITY)) { return x; } - return ((exp(x) - exp(-x)) / 2.0); + /* Ensure exp(x) does not overflow. Values of x less than 709.78 are guaranteed to be safe. */ + if (x < 709.78) { + return ((exp(x) - exp(-x)) / 2.0); + } + else { + y = sinh(x / 3.0); + return ((3.0 * y) + (4.0 * y * y * y)); + } } diff --git a/math/power.c b/math/power.c index 143a165a..9330b585 100644 --- a/math/power.c +++ b/math/power.c @@ -15,6 +15,7 @@ #include /* Needed for __ieee754_sqrt */ #include +#include #include #include #include @@ -27,11 +28,16 @@ double pow(double base, double exponent) double res = 1.0; int s = 1; + if ((base == 1.0) || (exponent == 0.0) || (exponent == -0.0)) { + return 1.0; + } + + if ((isnan(base) != 0) || (isnan(exponent) != 0)) { + return NAN; + } + if ((base == 0.0) || (base == -0.0)) { - if ((exponent == 0.0) || (exponent == -0.0)) { - return 1.0; - } - else if (exponent < 0.0) { + if (exponent < 0.0) { if (base == 0.0) { return INFINITY; } @@ -42,9 +48,6 @@ double pow(double base, double exponent) return 0.0; } - else if ((exponent == 0.0) || (exponent == -0.0)) { - return 1.0; - } if (base < 0.0) { if (!isInteger(exponent)) { @@ -52,7 +55,7 @@ double pow(double base, double exponent) return NAN; } - s = (fmod(exponent, 2) > 0.0) ? -1 : 1; + s = (int)((fmod(exponent, 2) > 0.0) ? -1 : 1); base = -base; } @@ -61,7 +64,12 @@ double pow(double base, double exponent) } exponent *= log(base); - res = s * exp(exponent); + res = (double)s * exp(exponent); + + /* ia32 use 80-bit extended precision registers so isinf() may not be true */ + if ((isinf(res) != 0) || (fabs(res) > DBL_MAX)) { + errno = ERANGE; + } return res; } diff --git a/math/trig.c b/math/trig.c index 376fadf9..eece1dca 100644 --- a/math/trig.c +++ b/math/trig.c @@ -109,14 +109,21 @@ double sin(double x) double tan(double x) { - double c = cos(x); + double c; - if ((c > 0.0) || (c < 0.0)) { - return (sin(x) / c); + if (isinf(x) != 0) { + errno = EDOM; + return NAN; } - errno = EDOM; - return NAN; + c = cos(x); + + if (c != 0.0) { + return (sin(x) / c); + } + else { + return 0.0; + } } @@ -197,6 +204,10 @@ double atan(double x) double atan2(double y, double x) { + if ((isnan(y) != 0) || (isnan(x) != 0)) { + return NAN; + } + if (x > 0.0) { return atan(y / x); } @@ -215,6 +226,5 @@ double atan2(double y, double x) return -M_PI_2; } - errno = EDOM; - return NAN; + return 0.0; }