Import Upstream version 5d6upstream/5d6

1 files changed, 212 insertions, 39 deletions

diff --git a/scl.c b/scl.c
index 51d9eab..57d020e 100644
--- a/scl.c
+++ b/scl.c
@@ -15,26 +15,26 @@
  * the Free Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111, USA.
  *
  * As a special exception, the Free Software Foundation gives permission
- * for additional uses of the text contained in its release of GUILE.
+ * for additional uses of the text contained in its release of SCM.
  *
- * The exception is that, if you link the GUILE library with other files
+ * The exception is that, if you link the SCM library with other files
  * to produce an executable, this does not by itself cause the
  * resulting executable to be covered by the GNU General Public License.
  * Your use of that executable is in no way restricted on account of
- * linking the GUILE library code into it.
+ * linking the SCM library code into it.
  *
  * This exception does not however invalidate any other reasons why
  * the executable file might be covered by the GNU General Public License.
  *
  * This exception applies only to the code released by the
- * Free Software Foundation under the name GUILE.  If you copy
+ * Free Software Foundation under the name SCM.  If you copy
  * code from other Free Software Foundation releases into a copy of
- * GUILE, as the General Public License permits, the exception does
+ * SCM, as the General Public License permits, the exception does
  * not apply to the code that you add in this way.  To avoid misleading
  * anyone as to the status of such modified files, you must delete
  * this exception notice from them.
  *
- * If you write modifications of your own for GUILE, it is your choice
+ * If you write modifications of your own for SCM, it is your choice
  * whether to permit this exception to apply to your modifications.
  * If you do not wish that, delete this exception notice.
  */
@@ -45,11 +45,19 @@
 #include "scm.h"
 
 #ifdef FLOATS
-# include <math.h>
+# ifndef PLAN9
+#  include <math.h>
+# endif
 
 static double big2scaldbl P((SCM b, int expt));
 static SCM bigdblop P((int op, SCM b, double re, double im));
 static SCM inex_divbigbig P((SCM a, SCM b));
+static int apx_log10 P((double x));
+static double lpow10 P((double x, int n));
+static sizet idbl2str P((double f, char *a));
+static sizet iflo2str P((SCM flt, char *str));
+static void add1 P((double f, double *fsum));
+static long scm_twos_power P((SCM n));
 
 static char s_makrect[] = "make-rectangular", s_makpolar[] = "make-polar",
 	    s_magnitude[] = "magnitude", s_angle[] = "angle",
@@ -83,7 +91,13 @@ static char s_intexpt[] = "integer-expt";
 
 /*** NUMBERS -> STRINGS ***/
 #ifdef FLOATS
-static int dbl_mant_dig;
+# ifndef DBL_MANT_DIG
+#  define DBL_MANT_DIG dbl_mant_dig
+# endif
+static int dbl_mant_dig = 0;
+static double max_dbl_int;       /* Integers less than or equal to max_dbl_int
+				    are representable exactly as doubles. */
+static double dbl_eps;
 double dbl_prec(x)
      double x;
 {
@@ -91,10 +105,10 @@ double dbl_prec(x)
   double frac = frexp(x, &expt);
 # ifdef DBL_MIN_EXP
   if (0.0==x || expt < DBL_MIN_EXP) /* gradual underflow */
-    return ldexp(1.0, -dbl_mant_dig) * ldexp(1.0, DBL_MIN_EXP);
+    return ldexp(1.0, -DBL_MANT_DIG) * ldexp(1.0, DBL_MIN_EXP);
 # endif
-  if (1.0==frac) return ldexp(1.0, expt - dbl_mant_dig + 1);
-  return ldexp(1.0, expt - dbl_mant_dig);
+  if (1.0==frac) return ldexp(1.0, expt - DBL_MANT_DIG + 1);
+  return ldexp(1.0, expt - DBL_MANT_DIG);
 }
 
 static double llog2 = 0.3010299956639812; /* log10(2) */
@@ -130,7 +144,7 @@ static double lpow10(x, n)
 
 /* DBL2STR_FUZZ is a somewhat arbitrary guard against
    round off error in scaling f and fprec. */
-#define DBL2STR_FUZZ 0.9
+# define DBL2STR_FUZZ 0.9
 int dblprec;
 static sizet idbl2str(f, a)
      double f;
@@ -479,14 +493,23 @@ SCM istr2int(str, len, radix)
       if (c >= radix) return BOOL_F; /* bad digit for radix */
       ln = n;
       n = n * radix - c;
-      /* Negation is a workaround for HP700 cc bug */
-      if (n > ln || (-n > -MOST_NEGATIVE_FIXNUM)) goto ovfl;
+      if (n > ln
+# ifdef hpux
+	  || (-n > -MOST_NEGATIVE_FIXNUM) /* workaround for HP700 cc bug */
+# endif
+	  ) goto ovfl;
       break;
     default:
       return BOOL_F;		/* not a digit */
     }
   } while (i < len);
-  if (!lead_neg) if ((n = -n) > MOST_POSITIVE_FIXNUM) goto ovfl;
+  if (lead_neg) {
+    if (n < MOST_NEGATIVE_FIXNUM) goto ovfl;
+  }
+  else {
+    if (n < -MOST_POSITIVE_FIXNUM) goto ovfl;
+    n = -n;
+  }
   return MAKINUM(n);
  ovfl:				/* overflow scheme integer */
   return BOOL_F;
@@ -494,6 +517,34 @@ SCM istr2int(str, len, radix)
 #endif
 
 #ifdef FLOATS
+# ifdef BIGDIG
+static char twostab[] = {4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0};
+static long scm_twos_power(n)
+     SCM n;
+{
+  long d, c = 0;
+  int d4;
+#  ifdef BIGDIG
+  if NINUMP(n) {
+    BIGDIG *ds;
+    int i = 0;
+    ds = BDIGITS(n);
+    while (0==(d = ds[i++])) c += BITSPERDIG;
+    goto out;
+  }
+#  endif
+  d = INUM(n);
+  if (0==d) return 0;
+ out:
+  do {
+    d4 = 15 & d;
+    c += twostab[d4];
+    d >>= 4;
+  } while (0==d4);
+  return c;
+}
+# endif /* def BIGDIG */
+
 SCM istr2flo(str, len, radix)
      register char *str;
      register long len;
@@ -965,6 +1016,29 @@ SCM numberp(x)
   return BOOL_F;
 }
 #ifdef FLOATS
+# ifdef BIGDIG
+int scm_bigdblcomp(b, d)
+     SCM b;
+     double d;
+{
+  sizet dlen, blen;
+  int dneg = d < 0 ? 1 : 0;
+  int bneg = BIGSIGN(b) ? 1 : 0;
+  if (bneg < dneg) return -1;
+  if (bneg > dneg) return 1;
+  frexp(d, &dlen);
+  blen = INUM(scm_intlength(b));
+  if (blen > dlen) return dneg ? 1 : -1;
+  if (blen < dlen) return dneg ? -1 : 1;
+  if ((blen <= DBL_MANT_DIG) || (blen - scm_twos_power(b)) <= DBL_MANT_DIG) {
+    double bd = big2dbl(b);
+    if (bd > d) return -1;
+    if (bd < d) return 1;
+    return 0;
+  }
+  return bigcomp(b, dbl2big(d));
+}
+# endif
 SCM realp(x)
      SCM x;
 {
@@ -1018,7 +1092,8 @@ SCM eqp(x, y)
       if BIGP(y) return (0==bigcomp(x, y)) ? BOOL_T : BOOL_F;
       ASRTGO(INEXP(y), bady);
     bigreal:
-      return (REALP(y) && (big2dbl(x)==REALPART(y))) ? BOOL_T : BOOL_F;
+      return (REALP(y) && (0==scm_bigdblcomp(x, REALPART(y)))) ?
+	BOOL_T : BOOL_F;
     }
     ASRTGO(INEXP(x), badx);
 # else
@@ -1091,7 +1166,7 @@ SCM lessp(x, y)
       ASRTGO(NIMP(y), bady);
       if BIGP(y) return (1==bigcomp(x, y)) ? BOOL_T : BOOL_F;
       ASRTGO(REALP(y), bady);
-      return (big2dbl(x) < REALPART(y)) ? BOOL_T : BOOL_F;
+      return (1==scm_bigdblcomp(x, REALPART(y))) ? BOOL_T : BOOL_F;
     }
     ASRTGO(REALP(x), badx);
 # else
@@ -1100,7 +1175,7 @@ SCM lessp(x, y)
     if INUMP(y) return (REALPART(x) < ((double)INUM(y))) ? BOOL_T : BOOL_F;
 # ifdef BIGDIG
     ASRTGO(NIMP(y), bady);
-    if BIGP(y) return (REALPART(x) < big2dbl(y)) ? BOOL_T : BOOL_F;
+    if BIGP(y) return (-1==scm_bigdblcomp(y, REALPART(x))) ? BOOL_T : BOOL_F;
     ASRTGO(REALP(y), bady);
 # else
     ASRTGO(NIMP(y) && REALP(y), bady);
@@ -1248,10 +1323,12 @@ SCM negativep(x)
   return (x < INUM0) ? BOOL_T : BOOL_F;
 }
 
+static char s_exactprob[] = "not representable as inexact";
 SCM lmax(x, y)
      SCM x, y;
 {
 #ifdef FLOATS
+  SCM t;
   double z;
 #endif
   if UNBNDP(y) {
@@ -1270,8 +1347,11 @@ SCM lmax(x, y)
       ASRTGO(NIMP(y), bady);
       if BIGP(y) return (1==bigcomp(x, y)) ? y : x;
       ASRTGO(REALP(y), bady);
+    big_dbl:
+      if (-1 != scm_bigdblcomp(x, REALPART(y))) return y;
       z = big2dbl(x);
-      return (z < REALPART(y)) ? y : makdbl(z, 0.0);
+      ASSERT(0==scm_bigdblcomp(x, z), x, s_exactprob, s_max);
+      return makdbl(z, 0.0);
     }
     ASRTGO(REALP(x), badx);
 # else
@@ -1280,7 +1360,9 @@ SCM lmax(x, y)
     if INUMP(y) return (REALPART(x) < (z = INUM(y))) ? makdbl(z, 0.0) : x;
 # ifdef BIGDIG
     ASRTGO(NIMP(y), bady);
-    if BIGP(y) return (REALPART(x) < (z = big2dbl(y))) ? makdbl(z, 0.0) : x;
+    if BIGP(y) {
+      t = y; y = x; x = t; goto big_dbl;
+    }
     ASRTGO(REALP(y), bady);
 # else
     ASRTGO(NIMP(y) && REALP(y), bady);
@@ -1330,6 +1412,7 @@ SCM lmin(x, y)
      SCM x, y;
 {
 #ifdef FLOATS
+  SCM t;
   double z;
 #endif
   if UNBNDP(y) {
@@ -1348,8 +1431,11 @@ SCM lmin(x, y)
       ASRTGO(NIMP(y), bady);
       if BIGP(y) return (-1==bigcomp(x, y)) ? y : x;
       ASRTGO(REALP(y), bady);
+    big_dbl:
+      if (1 != scm_bigdblcomp(x, REALPART(y))) return y;
       z = big2dbl(x);
-      return (z > REALPART(y)) ? y : makdbl(z, 0.0);
+      ASSERT(0==scm_bigdblcomp(x, z), x, s_exactprob, s_min);
+      return makdbl(z, 0.0);
     }
     ASRTGO(REALP(x), badx);
 # else
@@ -1358,7 +1444,9 @@ SCM lmin(x, y)
     if INUMP(y) return (REALPART(x) > (z = INUM(y))) ? makdbl(z, 0.0) : x;
 # ifdef BIGDIG
     ASRTGO(NIMP(y), bady);
-    if BIGP(y) return (REALPART(x) > (z = big2dbl(y))) ? makdbl(z, 0.0) : x;
+    if BIGP(y) {
+      t = y; y = x; x = t; goto big_dbl;
+    }
     ASRTGO(REALP(y), bady);
 # else
     ASRTGO(NIMP(y) && REALP(y), bady);
@@ -1960,6 +2048,10 @@ SCM scm_intexpt(z1, z2)
      SCM z1, z2;
 {
   SCM acc = MAKINUM(1L);
+  int recip = 0;
+#ifdef FLOATS
+  double dacc, dz1;
+#endif
 #ifdef BIGDIG
   if (INUM0==z1 || acc==z1) return z1;
   else if (MAKINUM(-1L)==z1) return BOOL_F==evenp(z2)?z1:acc;
@@ -1968,18 +2060,17 @@ SCM scm_intexpt(z1, z2)
   z2 = INUM(z2);
   if (z2 < 0) {
     z2 = -z2;
-    z1 = divide(z1, UNDEFINED);
+    recip = 1; /*    z1 = divide(z1, UNDEFINED); */
   }
   if INUMP(z1) {
     long tmp, iacc = 1, iz1 = INUM(z1);
+#ifdef FLOATS
+    if (recip) { dz1 = iz1; goto flocase; }
+#endif
     while (1) {
       if (0==z2) {
 	acc = long2num(iacc);
-#ifndef RECKLESS
-	if (FALSEP(z1))
-	  errout: wta(UNDEFINED, (char *)OVFLOW, s_intexpt);
-#endif
-	return acc;
+	break;
       }
       if (1==z2) {
 	tmp = iacc*iz1;
@@ -1991,8 +2082,7 @@ SCM scm_intexpt(z1, z2)
 	  goto gencase;
 	}
 	acc = long2num(tmp);
-	ASRTGO(NFALSEP(acc), errout);
-	return acc;
+	break;
       }
       if (z2 & 1) {
 	tmp = iacc*iz1;
@@ -2005,28 +2095,37 @@ SCM scm_intexpt(z1, z2)
       iz1 = tmp;
       z2 >>= 1;
     }
+#ifndef RECKLESS
+    if (FALSEP(acc))
+      errout: wta(UNDEFINED, (char *)OVFLOW, s_intexpt);
+#endif
+    goto ret;
   }
   ASSERT(NIMP(z1), z1, ARG1, s_intexpt);
-#ifdef FLOATS			/* Move to scl.c ? */
+#ifdef FLOATS
   if REALP(z1) {
-    double dacc = 1.0, dz1 = REALPART(z1);
+    dz1 = REALPART(z1);
+  flocase:
+    dacc = 1.0;
     while(1) {
-      if (0==z2) return makdbl(dacc, 0.0);
-      if (1==z2) return makdbl(dacc*dz1, 0.0);
+      if (0==z2) break;
+      if (1==z2) {dacc = dacc*dz1; break;}
       if (z2 & 1) dacc = dacc*dz1;
       dz1 = dz1*dz1;
       z2 >>= 1;
     }
+    return makdbl(recip ? 1.0/dacc : dacc, 0.0);
   }
 #endif
  gencase:
   while(1) {
-    if (0==z2) return acc;
-    if (1==z2) return product(acc, z1);
+    if (0==z2) break;
+    if (1==z2) {acc = product(acc, z1); break;}
     if (z2 & 1) acc = product(acc, z1);
     z1 = product(z1, z1);
     z2 >>= 1;
   }
+ ret: return recip ? divide(acc, UNDEFINED) : acc;
 }
 
 #ifdef FLOATS
@@ -2242,8 +2341,11 @@ SCM in2ex(z)
 # ifdef BIGDIG
   {
     double u = floor(REALPART(z)+0.5);
-    if ((u <= MOST_POSITIVE_FIXNUM) && (-u <= -MOST_NEGATIVE_FIXNUM)) {
-      /* Negation is a workaround for HP700 cc bug */
+    if ((u <= MOST_POSITIVE_FIXNUM)
+#  ifdef hpux
+	&& (-u <= -MOST_NEGATIVE_FIXNUM) /* workaround for HP700 cc bug */
+#  endif
+	) {
       SCM ans = MAKINUM((long)u);
       if (INUM(ans)==(long)u) return ans;
     }
@@ -2349,6 +2451,8 @@ static SCM bigdblop(op, b, re, im)
   }
   case '\\':
     return makdbl(ldexp(re/bm, -i), 0.0==im ? 0.0 : ldexp(im/bm, -i));
+  default:
+    return UNSPECIFIED;
   }
 }
 static SCM inex_divbigbig(a, b)
@@ -2371,6 +2475,64 @@ static SCM inex_divbigbig(a, b)
   }
   return makdbl(r, 0.0);
 }
+
+static char s_dfloat_parts[] = "double-float-parts";
+SCM scm_dfloat_parts(f)
+     SCM f;
+{
+  int expt, ndig = DBL_MANT_DIG;
+  double mant = frexp(num2dbl(f, (char *)ARG1, s_dfloat_parts), &expt);
+#  ifdef DBL_MIN_EXP
+  if (expt < DBL_MIN_EXP)
+    ndig -= DBL_MIN_EXP - expt;
+#  endif
+  mant *= ldexp(1.0, ndig);
+  expt -= ndig;
+  return scm_values(dbl2big(mant), MAKINUM(expt), EOL, s_dfloat_parts);
+}
+static char s_make_dfloat[] = "make-double-float";
+SCM scm_make_dfloat(mant, expt)
+     SCM mant, expt;
+{
+  double dmant = num2dbl(mant, (char *)ARG1, s_make_dfloat);
+  int e = INUM(expt);
+  ASSERT(INUMP(expt), expt, ARG2, s_make_dfloat);
+  ASSERT((dmant < 0 ? -dmant : dmant)<=max_dbl_int, mant,
+	 OUTOFRANGE, s_make_dfloat);
+  return makdbl(ldexp(dmant, e), 0.0);
+}
+static char s_next_dfloat[] = "next-double-float";
+SCM scm_next_dfloat(f1, f2)
+     SCM f1, f2;
+{
+  int e, neg = 0;
+  double d1 = num2dbl(f1, (char *)ARG1, s_next_dfloat);
+  double dif = num2dbl(f2, (char *)ARG2, s_next_dfloat) - d1;
+  double d = frexp(d1, &e), eps = dbl_eps;
+  if (d1 < 0) {neg = 1; dif = -dif; d = -d;}
+  if (dif > 0) {
+#  ifdef DBL_MIN_EXP
+    if (e < DBL_MIN_EXP)
+      eps = ldexp(eps, DBL_MIN_EXP - e);
+    else if (0.0==d)
+      eps = ldexp(1.0, DBL_MIN_EXP - DBL_MANT_DIG);
+#  endif
+    d = ldexp(d + eps, e);
+  }
+  else if (dif < 0) {
+#  ifdef DBL_MIN_EXP
+    if (e < DBL_MIN_EXP)
+      eps = ldexp(eps, DBL_MIN_EXP - e);
+    else if (0.0==d)
+      eps = ldexp(-1.0, DBL_MIN_EXP - DBL_MANT_DIG);
+#  endif
+    if (0.5==d) eps *= 0.5;
+    d = ldexp(d - eps, e);
+  }
+  else if (0.0==dif)
+    return f1;
+  return makdbl(neg ? -d : d, 0.0);
+}
 # endif
 #endif
 
@@ -2490,6 +2652,9 @@ static iproc subr1s[] = {
 	{s_angle, angle},
 	{s_in2ex, in2ex},
 	{s_ex2in, ex2in},
+# ifdef BIGDIG
+	{s_dfloat_parts, scm_dfloat_parts},
+# endif
 #else
 	{"real?", numberp},
 	{"rational?", numberp},
@@ -2524,6 +2689,10 @@ static iproc subr2s[] = {
 	{s_makpolar, makpolar},
 	{s_atan2, latan2},
 	{s_expt, expt},
+# ifdef BIGDIG
+	{s_make_dfloat, scm_make_dfloat},
+	{s_next_dfloat, scm_next_dfloat},
+# endif
 #endif
 #ifdef INUMS_ONLY
 	{s_memv, memq},
@@ -2634,7 +2803,7 @@ void init_scl()
 # ifdef DBL_MANT_DIG
   dbl_mant_dig = DBL_MANT_DIG;
 # else
-  {
+  if (!DBL_MANT_DIG) {		/* means we #defined it. */
     double fsum = 0.0, eps = 1.0;
     int i = 0;
     while (fsum != 1.0) {
@@ -2645,5 +2814,9 @@ void init_scl()
     dbl_mant_dig = i;
   }
 # endif /* DBL_MANT_DIG */
+  max_dbl_int = pow(2.0, dbl_mant_dig - 1.0);
+  max_dbl_int = max_dbl_int + (max_dbl_int - 1.0);
+  dbl_eps = ldexp(1.0, -dbl_mant_dig);
+  sysintern("double-float-mantissa-length", MAKINUM(DBL_MANT_DIG));
 #endif
 }