Merge pull request #70 from alexfru/master

Floating point improvements

include/math.h

@@ -4,16 +4,16 @@
  * specifies the terms and conditions for redistribution.
  */
 
-double  fabs(), floor(), ceil(), fmod(), ldexp();
-double  sqrt(), hypot(), atof();
-double  sin(), cos(), tan(), asin(), acos(), atan(), atan2();
-double  exp(), log(), log10(), pow();
-double  sinh(), cosh(), tanh();
-double  gamma();
-double  j0(), j1(), jn(), y0(), y1(), yn();
+double  fabs(double), floor(double), ceil(double);
+double  sqrt(double), hypot(double, double);
+double  sin(double), cos(double), tan(double);
+double  asin(double), acos(double), atan(double), atan2(double, double);
+double  exp(double), log(double), log10(double), pow(double, double);
+double  sinh(double), cosh(double), tanh(double);
+double  j0(double), j1(double), jn(int, double);
+double  y0(double), y1(double), yn(int, double);
 
-#define HUGE        1.701411733192644270e38
-#define LOGHUGE     39
+#define HUGE_VAL    3.40282347e+38 /* TBD??? use infinity? */
 
 int isnanf(float x);
 int isnan(double x);

include/stdlib.h

@@ -96,12 +96,10 @@ long    random (void);
 char    *setstate (char *);
 void    srandom (unsigned);
 
-#ifndef __SMALLER_C__
 double  atof (const char *);
 double  strtod (const char *, char **);
 char    *ecvt (double, int, int *, int *);
 char    *fcvt (double, int, int *, int *);
 char    *gcvt (double, int, char *);
-#endif
 
 #endif /* _STDLIB_H_ */

src/cmd/smlrc/smlrc.c

@@ -672,7 +672,7 @@ int ParamLevel = 0; // 1+ if parsing params, 0 otherwise
 
 unsigned char SyntaxStack0[SYNTAX_STACK_MAX];
 int SyntaxStack1[SYNTAX_STACK_MAX];
-int SyntaxStackCnt = 8; // number of explicitly initialized elements in SyntaxStack0[]
+int SyntaxStackCnt;
 
 // all code
 
@@ -9234,6 +9234,7 @@ int main(int argc, char** argv)
     tokChar
   }; // SyntaxStackCnt must be initialized to the number of elements in SyntaxStackInit[][]
   memcpy(SyntaxStack0, SyntaxStackInit, sizeof SyntaxStackInit);
+  SyntaxStackCnt = sizeof SyntaxStackInit / sizeof SyntaxStackInit[0];
 
 #ifdef __SMALLER_C__
 #ifdef DETERMINE_VA_LIST

src/libc/gen/atof.c

@@ -4,16 +4,27 @@
  * All rights reserved.  The Berkeley software License Agreement
  * specifies the terms and conditions for redistribution.
  */
+#include <float.h>
 #include <math.h>
 #include <ctype.h>
 
+/*
+ * BIG = 2**(DBL_MANT_DIG+3) defines how many decimal digits
+ * to take into account from the input. It doesn't make sense
+ * to use more digits than log10(2**DBL_MANT_DIG)+1.
+ * BIG is equal 2**27 or 2**56, depending on whether double
+ * is single or double precision.
+ */
+#define BIG (8 * (double)(1L << (DBL_MANT_DIG/2)) * \
+        (double)(1L << (DBL_MANT_DIG/2 + DBL_MANT_DIG%2)))
+
 double
 atof(p)
 register char *p;
 {
 	register int c;
 	double fl, flexp, exp5;
-	double big = 72057594037927936.;  /*2^56*/
+	double big = BIG;
 	int nd;
 	register int eexp, exp, neg, negexp, bexp;
 
@@ -73,7 +84,7 @@ register char *p;
 	}
 
 
-	if ((nd+exp*negexp) < -LOGHUGE){
+	if ((nd+exp*negexp) < DBL_MIN_10_EXP - 2) {
 		fl = 0;
 		exp = 0;
 	}

src/libc/gen/ldexp.c

@@ -1,3 +1,4 @@
+#include <float.h>
 #include <math.h>
 
 double
@@ -5,7 +6,6 @@ ldexp(fr, exp)
         double fr;
         int exp;
 {
-        double	huge	= 1.701411834604692293e38;
 	int neg;
 	int i;
 
@@ -20,13 +20,13 @@ ldexp(fr, exp)
 		i = i-1;
 	}
 	exp = exp+i;
-	if (exp > 127) {
+	if (exp >= DBL_MAX_EXP) {
 		if (neg)
-			return(-huge);
+			return(-HUGE_VAL);
 		else
-			return(huge);
+			return(HUGE_VAL);
         }
-	if (exp < -127)
+	if (exp < DBL_MIN_EXP - 2)
 		return(0);
 	while (exp > 30) {
 		fr = fr*(1L<<30);

src/libc/stdio/doprnt.c

@@ -41,8 +41,10 @@
 #include <float.h>
 #include <math.h>
 
-/* Max number conversion buffer length: a long in base 2, plus NUL byte. */
-#define MAXNBUF	(sizeof(long) * 8 + 1)
+/* Max number conversion buffer length. */
+#define MAXNBUF	\
+	(1/*sign*/ + DBL_MAX_10_EXP+1/*max integral digits*/ + \
+	1/*.*/ + DBL_DIG+1/*max fractional digits*/ + 1/*NUL*/)
 
 static unsigned char *ksprintn (unsigned char *buf, unsigned long v, unsigned char base,
 	int width, unsigned char *lp);
@@ -403,7 +405,7 @@ number:		if (sign && ((long) ul != 0L)) {
 					nbuf [size + 1] = 0;
 				}
 			}
-			if (neg)
+			if (neg || sign)
 				size++;
 			if (! ladjust && width && padding == ' ' &&
 			    (width -= size) > 0)
@@ -411,8 +413,11 @@ number:		if (sign && ((long) ul != 0L)) {
 					PUTC (' ');
 				} while (--width > 0);
 
-			if (neg)
+			if (neg) {
 				PUTC ('-');
+			} else if (sign) {
+				PUTC ('+');
+			}
 
 			if (! ladjust && width && (width -= size) > 0)
 				do {

src/libm/exp.c

@@ -4,6 +4,7 @@
  *
  * The coefficients are #1069 from Hart and Cheney. (22.35D)
  */
+#include <float.h>
 #include <errno.h>
 #include <math.h>
 
@@ -16,7 +17,7 @@ static double	q1	= .3277251518082914423057964422e6;
 static double	q2	= .1749287689093076403844945335e4;
 static double	log2e	= 1.4426950408889634073599247;
 static double	sqrt2	= 1.4142135623730950488016887;
-static double	maxf	= 10000;
+static double	maxf	= DBL_MAX_10_EXP * 2.5/*>ln(10)*/;
 
 double
 exp(arg)
@@ -32,7 +33,7 @@ double arg;
 		return(0.);
 	if(arg > maxf) {
 		errno = ERANGE;
-		return(HUGE);
+		return(HUGE_VAL);
 	}
 	arg *= log2e;
 	ent = floor(arg);

src/libm/j0.c

@@ -176,7 +176,7 @@ y0(arg)
 	errno = 0;
 	if(arg <= 0.){
 		errno = EDOM;
-		return(-HUGE);
+		return(-HUGE_VAL);
 	}
 	if(arg > 8.){
 		asympt(arg);

src/libm/j1.c

@@ -182,7 +182,7 @@ y1(arg)
 	x = arg;
 	if(x <= 0.){
 		errno = EDOM;
-		return(-HUGE);
+		return(-HUGE_VAL);
 	}
 	if(x > 8.){
 		asympt(x);

src/libm/jn.c

@@ -82,7 +82,7 @@ yn(n,x) int n; double x;{
 
 	if (x <= 0) {
 		errno = EDOM;
-		return(-HUGE);
+		return(-HUGE_VAL);
 	}
 	sign = 1;
 	if(n<0){

src/libm/log.c

@@ -31,7 +31,7 @@ double arg;
 
 	if(arg <= 0.) {
 		errno = EDOM;
-		return(-HUGE);
+		return(-HUGE_VAL);
 	}
 	x = frexp(arg,&exp);
 	while(x<0.5) {

src/libm/tan.c

@@ -64,8 +64,8 @@ tan(arg)
 		if(temp == 0.) {
 			errno = ERANGE;
 			if (sign>0)
-				return(HUGE);
-			return(-HUGE);
+				return(HUGE_VAL);
+			return(-HUGE_VAL);
 		}
 		temp = 1./temp;
 	}