math: improve POSIX compliance

JIRA: RTOS-974

math/common.c

@@ -14,6 +14,9 @@
  */
 
 #include <stdint.h>
+#include <errno.h>
+#include <math.h>
+#include <float.h>
 #include "common.h"
 
 
@@ -127,6 +130,10 @@ double quickPow(double x, int e)
 		}
 	}
 
+	if (res < DBL_MIN || res > DBL_MAX) {
+		errno = ERANGE;
+	}
+
 	return res;
 }
 

math/exp.c

@@ -22,9 +22,13 @@
 double frexp(double x, int* exp)
 {
 	conv_t *conv = (conv_t *)&x;
-
 	*exp = 0;
 
+	if (isnan(x) != 0) {
+		errno = EDOM;
+		return NAN;
+	}
+
 	if (x == 0.0 || x == -0.0)
 		return x;
 
@@ -43,6 +47,11 @@ double ldexp(double x, int exp)
 	conv_t *conv = (conv_t *)&x;
 	int exponent = 0;
 
+	if (isnan(x) != 0) {
+		errno = EDOM;
+		return NAN;
+	}
+
 	if (x == 0.0 || x == -0.0)
 		return x;
 
@@ -75,7 +84,7 @@ double log(double x)
 	conv_t *conv = (conv_t *)&tmp;
 	int exp = 0, i;
 
-	if (x < 0) {
+	if (isnan(x) != 0 || x < 0) {
 		errno = EDOM;
 		return NAN;
 	}
@@ -127,6 +136,11 @@ double modf(double x, double* intpart)
 	int exp = conv->i.exponent - 1023;
 	uint64_t m, mask = 0xfffffffffffffLL;
 
+	if (isnan(x) != 0) {
+		errno = EDOM;
+		return NAN;
+	}
+
 	if (exp > 52) {
 		*intpart = x;
 		return conv->i.sign? -0.0 : 0.0;
@@ -170,8 +184,16 @@ double exp(double x)
 	double res, resi, powx, e, factorial;
 	int i;
 
-	if (x > 710.0)
+	if (isnan(x) != 0) {
+		errno = EDOM;
+		return NAN;
+	}
+
+	/* 710.47 is intended - bare in mind cosh */
+	if (x > 710.47) {
+		errno = ERANGE;
 		return HUGE_VAL;
+	}
 
 	/* Get floor of exponent */
 	x = modf(x, &e);
@@ -200,6 +222,11 @@ double ceil(double x)
 {
 	double ipart, fpart;
 
+	if (isnan(x) != 0) {
+		errno = EDOM;
+		return NAN;
+	}
+
 	fpart = modf(x, &ipart);
 
 	if (x > 0 && fpart + x != x)
@@ -213,6 +240,11 @@ double floor(double x)
 {
 	double ipart, fpart;
 
+	if (isnan(x) != 0) {
+		errno = EDOM;
+		return NAN;
+	}
+
 	fpart = modf(x, &ipart);
 
 	if (x < 0 && fpart + x != x)

math/hyper.c

@@ -14,10 +14,16 @@
  */
 
 #include <math.h>
+#include <errno.h>
 
 
 double cosh(double x)
 {
+	if (isnan(x) != 0) {
+		errno = EDOM;
+		return NAN;
+	}
+
 	if (x == INFINITY || x == -INFINITY)
 		return INFINITY;
 
@@ -27,6 +33,11 @@ double cosh(double x)
 
 double sinh(double x)
 {
+	if (isnan(x) != 0) {
+		errno = EDOM;
+		return NAN;
+	}
+
 	if (x == INFINITY || x == -INFINITY)
 		return x;
 

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,17 @@ 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) {
+		errno = EDOM;
+		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;
 			else if (base == -0.0)
@@ -40,9 +47,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)) {
@@ -60,6 +64,10 @@ double pow(double base, double exponent)
 	exponent *= log(base);
 	res = s * exp(exponent);
 
+	if (res < DBL_MIN || res > DBL_MAX) {
+		errno = ERANGE;
+	}
+
 	return res;
 }
 

math/trig.c

@@ -180,6 +180,11 @@ double atan(double x)
 
 double atan2(double y, double x)
 {
+	if (isnan(y) != 0 || isnan(x) != 0) {
+		errno = EDOM;
+		return NAN;
+	}
+
 	if (x > 0) {
 		return atan(y / x);
 	}