Browse Source

math: fix exp10 not to raise invalid exception on NaN

This was not caught earlier because gcc incorrectly generates quiet
relational operators that never raise exceptions.
Szabolcs Nagy 10 years ago
parent
commit
18daae3135
3 changed files with 13 additions and 4 deletions
  1. 4 1
      src/math/exp10.c
  2. 4 1
      src/math/exp10f.c
  3. 5 2
      src/math/exp10l.c

+ 4 - 1
src/math/exp10.c

@@ -1,5 +1,6 @@
 #define _GNU_SOURCE
 #define _GNU_SOURCE
 #include <math.h>
 #include <math.h>
+#include <stdint.h>
 #include "libc.h"
 #include "libc.h"
 
 
 double exp10(double x)
 double exp10(double x)
@@ -11,7 +12,9 @@ double exp10(double x)
 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15
 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15
 	};
 	};
 	double n, y = modf(x, &n);
 	double n, y = modf(x, &n);
-	if (fabs(n) < 16) {
+	union {double f; uint64_t i;} u = {n};
+	/* fabs(n) < 16 without raising invalid on nan */
+	if ((u.i>>52 & 0x7ff) < 0x3ff+4) {
 		if (!y) return p10[(int)n+15];
 		if (!y) return p10[(int)n+15];
 		y = exp2(3.32192809488736234787031942948939 * y);
 		y = exp2(3.32192809488736234787031942948939 * y);
 		return y * p10[(int)n+15];
 		return y * p10[(int)n+15];

+ 4 - 1
src/math/exp10f.c

@@ -1,5 +1,6 @@
 #define _GNU_SOURCE
 #define _GNU_SOURCE
 #include <math.h>
 #include <math.h>
+#include <stdint.h>
 #include "libc.h"
 #include "libc.h"
 
 
 float exp10f(float x)
 float exp10f(float x)
@@ -9,7 +10,9 @@ float exp10f(float x)
 		1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7
 		1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7
 	};
 	};
 	float n, y = modff(x, &n);
 	float n, y = modff(x, &n);
-	if (fabsf(n) < 8) {
+	union {float f; uint32_t i;} u = {n};
+	/* fabsf(n) < 8 without raising invalid on nan */
+	if ((u.i>>23 & 0xff) < 0x7f+3) {
 		if (!y) return p10[(int)n+7];
 		if (!y) return p10[(int)n+7];
 		y = exp2f(3.32192809488736234787031942948939f * y);
 		y = exp2f(3.32192809488736234787031942948939f * y);
 		return y * p10[(int)n+7];
 		return y * p10[(int)n+7];

+ 5 - 2
src/math/exp10l.c

@@ -2,13 +2,14 @@
 #include <float.h>
 #include <float.h>
 #include <math.h>
 #include <math.h>
 #include "libc.h"
 #include "libc.h"
+#include "libm.h"
 
 
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
 long double exp10l(long double x)
 long double exp10l(long double x)
 {
 {
 	return exp10(x);
 	return exp10(x);
 }
 }
-#else
+#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
 long double exp10l(long double x)
 long double exp10l(long double x)
 {
 {
 	static const long double p10[] = {
 	static const long double p10[] = {
@@ -18,7 +19,9 @@ long double exp10l(long double x)
 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15
 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15
 	};
 	};
 	long double n, y = modfl(x, &n);
 	long double n, y = modfl(x, &n);
-	if (fabsl(n) < 16) {
+	union ldshape u = {n};
+	/* fabsl(n) < 16 without raising invalid on nan */
+	if ((u.i.se & 0x7fff) < 0x3fff+4) {
 		if (!y) return p10[(int)n+15];
 		if (!y) return p10[(int)n+15];
 		y = exp2l(3.32192809488736234787031942948939L * y);
 		y = exp2l(3.32192809488736234787031942948939L * y);
 		return y * p10[(int)n+15];
 		return y * p10[(int)n+15];