math: improve POSIX compliance

JIRA: RTOS-974

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

math/common.c

@@ -14,6 +14,9 @@
  */
 
 #include <stdint.h>
+#include <errno.h>
+#include <math.h>
+#include <float.h>
 #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;
 }
 

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)) {

math/hyper.c

@@ -14,25 +14,53 @@
  */
 
 #include <math.h>
+#include <float.h>
+#include <errno.h>
 
 
 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));
+	}
 }
 
 

math/power.c

@@ -15,6 +15,7 @@
 
 #include <arch.h> /* Needed for __ieee754_sqrt */
 #include <math.h>
+#include <float.h>
 #include <limits.h>
 #include <errno.h>
 #include <stdlib.h>
@@ -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,17 +48,14 @@ 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)) {
 			errno = EDOM;
 			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;
 }

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;
 }