Check for Nan

lang/cem/libcc.ansi/math/.distr

@@ -21,3 +21,4 @@ sinh.c
 sqrt.c
 tan.c
 tanh.c
+isnan.c

lang/cem/libcc.ansi/math/LIST

@@ -19,3 +19,4 @@ floor.c
 hugeval.c
 frexp.e
 modf.e
+isnan.c

lang/cem/libcc.ansi/math/Makefile

@@ -1,4 +1,4 @@
 clean:
 	rm -f asin.o atan2.o atan.o ceil.o fabs.o pow.o log10.o \
 		log.o sin.o sinh.o sqrt.o tan.o tanh.o exp.o ldexp.o \
-		fmod.o floor.o hugeval.o frexp.o modf.o OLIST
+		fmod.o floor.o hugeval.o frexp.o modf.o isnan.o OLIST

lang/cem/libcc.ansi/math/asin.c

@@ -32,6 +32,11 @@ asin_acos(double x, int cosfl)
 		 1.0
 	};
 
+	if (__IsNan(x)) {
+		errno = EDOM;
+		return x;
+	}
+
 	if (negative) {
 		x = -x;
 	}

lang/cem/libcc.ansi/math/atan.c

@@ -8,6 +8,7 @@
 
 #include	<float.h>
 #include	<math.h>
+#include	<errno.h>
 #include	"localmath.h"
 
 double
@@ -42,6 +43,10 @@ atan(double x)
 	int     n;
 	double  g;
 
+	if (__IsNan(x)) {
+		errno = EDOM;
+		return x;
+	}
 	if (neg) {
 		x = -x;
 	}

lang/cem/libcc.ansi/math/exp.c

@@ -36,6 +36,10 @@ exp(double x)
 	int	n;
 	int	negative = x < 0;
 
+	if (__IsNan(x)) {
+		errno = EDOM;
+		return x;
+	}
 	if (x < M_LN_MIN_D) {
 		errno = ERANGE;
 		return 0.0;

lang/cem/libcc.ansi/math/ldexp.c

@@ -14,11 +14,19 @@ ldexp(double fl, int exp)
 	int sign = 1;
 	int currexp;
 
+	if (__IsNan(fl)) {
+		errno = EDOM;
+		return fl;
+	}
 	if (fl == 0.0) return 0.0;
 	if (fl<0) {
 		fl = -fl;
 		sign = -1;
 	}
+	if (fl > DBL_MAX) {		/* for infinity */
+		errno = ERANGE;
+		return sign * fl;
+	}
 	fl = frexp(fl,&currexp);
 	exp += currexp;
 	if (exp > 0) {

lang/cem/libcc.ansi/math/log.c

@@ -7,6 +7,7 @@
 /* $Header$ */
 
 #include	<math.h>
+#include	<float.h>
 #include	<errno.h>
 #include	"localmath.h"
 
@@ -32,6 +33,10 @@ log(double x)
 	double	znum, zden, z, w;
 	int	exponent;
 
+	if (__IsNan(x)) {
+		errno = EDOM;
+		return x;
+	}
 	if (x < 0) {
 		errno = EDOM;
 		return -HUGE_VAL;
@@ -41,6 +46,9 @@ log(double x)
 		return -HUGE_VAL;
 	}
 
+	if (x <= DBL_MAX) {
+	}
+	else return x;	/* for infinity and Nan */
 	x = frexp(x, &exponent);
 	if (x > M_1_SQRT2) {
 		znum = (x - 0.5) - 0.5;

lang/cem/libcc.ansi/math/log10.c

@@ -13,6 +13,10 @@
 double
 log10(double x)
 {
+	if (__IsNan(x)) {
+		errno = EDOM;
+		return x;
+	}
 	if (x < 0) {
 		errno = EDOM;
 		return -HUGE_VAL;

lang/cem/libcc.ansi/math/sin.c

@@ -7,6 +7,8 @@
 /* $Header$ */
 
 #include	<math.h>
+#include	<float.h>
+#include	<errno.h>
 #include	"localmath.h"
 
 static double
@@ -30,14 +32,18 @@ sinus(double x, int cos_flag)
 
 	double  xsqr;
 	double  y;
-	int     neg = 0;
+	int     neg = 1;
 
+	if (__IsNan(x)) {
+		errno = EDOM;
+		return x;
+	}
 	if (x < 0) {
 		x = -x;
-		neg = 1;
+		neg = -1;
 	}
 	if (cos_flag) {
-		neg = 0;
+		neg = 1;
 		y = M_PI_2 + x;
 	}
 	else    y = x;
@@ -46,6 +52,8 @@ sinus(double x, int cos_flag)
 
 	y = y * M_1_PI + 0.5;
 
+	if (y >= DBL_MAX/M_PI) return 0.0;
+
 	/*      Use extended precision to calculate reduced argument.
 		Here we used 12 bits of the mantissa for a1.
 		Also split x in integer part x1 and fraction part x2.
@@ -56,7 +64,7 @@ sinus(double x, int cos_flag)
 		double x1, x2;
 
 		modf(y, &y);
-		if (modf(0.5*y, &x1)) neg = !neg;
+		if (modf(0.5*y, &x1)) neg = -neg;
 		if (cos_flag) y -= 0.5;
 		x2 = modf(x, &x1);
 		x = x1 - y * A1;
@@ -67,7 +75,7 @@ sinus(double x, int cos_flag)
 	}
  
 	if (x < 0) {
-		neg = !neg;
+		neg = -neg;
 		x = -x;
 	}
 
@@ -75,7 +83,7 @@ sinus(double x, int cos_flag)
 
 	y = x * x;
 	x += x * y * POLYNOM7(y, r);
-	return neg ? -x : x;
+	return neg==-1 ? -x : x;
 }
 
 double

lang/cem/libcc.ansi/math/sinh.c

@@ -34,6 +34,10 @@ sinh_cosh(double x, int cosh_flag)
 	int	negative = x < 0;
 	double	y = negative ? -x : x;
 
+	if (__IsNan(x)) {
+		errno = EDOM;
+		return x;
+	}
 	if (! cosh_flag && y <= 1.0) {
 		/* ??? check for underflow ??? */
 		y = y * y;

lang/cem/libcc.ansi/math/sqrt.c

@@ -7,6 +7,7 @@
 /* $Header$ */
 
 #include	<math.h>
+#include	<float.h>
 #include	<errno.h>
 
 #define NITER	5
@@ -17,11 +18,17 @@ sqrt(double x)
 	int exponent;
 	double val;
 
+	if (__IsNan(x)) {
+		errno = EDOM;
+		return x;
+	}
 	if (x <= 0) {
 		if (x < 0) errno = EDOM;
 		return 0;
 	}
 
+	if (x > DBL_MAX) return x;	/* for infinity */
+
 	val = frexp(x, &exponent);
 	if (exponent & 1) {
 		exponent--;

lang/cem/libcc.ansi/math/tan.c

@@ -7,6 +7,8 @@
 /* $Header$ */
 
 #include	<math.h>
+#include	<float.h>
+#include	<errno.h>
 #include	"localmath.h"
 
 double
@@ -34,12 +36,18 @@ tan(double x)
 		 0.49819433993786512270e-6
 	};
 
+	if (__IsNan(x)) {
+		errno = EDOM;
+		return x;
+	}
 	if (negative) x = -x;
  
 	/* ??? avoid loss of significance, error if x is too large ??? */
 
 	y = x * M_2_PI + 0.5;
 
+	if (y >= DBL_MAX/M_PI_2) return 0.0;
+
 	/*      Use extended precision to calculate reduced argument.
 		Here we used 12 bits of the mantissa for a1.
 		Also split x in integer part x1 and fraction part x2.

lang/cem/libcc.ansi/math/tanh.c

@@ -8,6 +8,7 @@
 
 #include	<float.h>
 #include	<math.h>
+#include	<errno.h>
 #include	"localmath.h"
 
 double
@@ -31,6 +32,10 @@ tanh(double x)
 	};
 	int 	negative = x < 0;
 
+	if (__IsNan(x)) {
+		errno = EDOM;
+		return x;
+	}
 	if (negative) x = -x;
 
 	if (x >= 0.5*M_LN_MAX_D) {