Import Upstream version 5d2upstream/5d2

1 files changed, 277 insertions, 64 deletions

diff --git a/scl.c b/scl.c
index 1b4507d..51d9eab 100644
--- a/scl.c
+++ b/scl.c
@@ -1,18 +1,18 @@
 /* Copyright (C) 1990, 1991, 1992, 1993, 1994, 1995, 1997 Free Software Foundation, Inc.
- * 
+ *
  * This program is free software; you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation; either version 2, or (at your option)
  * any later version.
- * 
+ *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU General Public License
  * along with this software; see the file COPYING.  If not, write to
- * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+ * 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.
@@ -36,7 +36,7 @@
  *
  * If you write modifications of your own for GUILE, it is your choice
  * whether to permit this exception to apply to your modifications.
- * If you do not wish that, delete this exception notice.  
+ * If you do not wish that, delete this exception notice.
  */
 
 /* "scl.c" non-IEEE utility functions and non-integer arithmetic.
@@ -47,6 +47,10 @@
 #ifdef FLOATS
 # include <math.h>
 
+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 char s_makrect[] = "make-rectangular", s_makpolar[] = "make-polar",
 	    s_magnitude[] = "magnitude", s_angle[] = "angle",
 	    s_real_part[] = "real-part", s_imag_part[] = "imag-part",
@@ -74,24 +78,69 @@ static char s_list_tail[] = "list-tail";
 static char s_str2list[] = "string->list";
 static char s_st_copy[] = "string-copy", s_st_fill[] = "string-fill!";
 static char s_vect2list[] = "vector->list", s_ve_fill[] = "vector-fill!";
+static char s_intexpt[] = "integer-expt";
+
 
 /*** NUMBERS -> STRINGS ***/
 #ifdef FLOATS
-int dblprec;
-static double fx[] = {0.0, 5e-1, 5e-2, 5e-3, 5e-4, 5e-5,
-			   5e-6, 5e-7, 5e-8, 5e-9, 5e-10,
-			   5e-11,5e-12,5e-13,5e-14,5e-15,
-			   5e-16,5e-17,5e-18,5e-19,5e-20};
+static int dbl_mant_dig;
+double dbl_prec(x)
+     double x;
+{
+  int expt;
+  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);
+# endif
+  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) */
+static int apx_log10(x)
+     double x;
+{
+  int expt;
+  double frac = frexp(x, &expt);
+  expt -= 1;
+  if (expt >= 0)
+    return (int)(expt * llog2);
+  return -((int)( -expt * llog2));
+}
+
+static double p10[] = {1.0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7};
+static double lpow10(x, n)
+     double x;
+     int n;
+{
+  if (n >= 0) {
+    while (n > 7) {
+      x *= 1e8;
+      n -= 8;
+    }
+    return x*p10[n];
+  }
+  while (n < -7) {
+    x /= 1e8;
+    n += 8;
+  }
+  return x/p10[-n];
+}
+
+/* DBL2STR_FUZZ is a somewhat arbitrary guard against
+   round off error in scaling f and fprec. */
+#define DBL2STR_FUZZ 0.9
+int dblprec;
 static sizet idbl2str(f, a)
      double f;
      char *a;
 {
-  int efmt, dpt, d, i, wp = dblprec;
+  double fprec = dbl_prec(f);
+  int efmt, dpt, d, i, exp;
   sizet ch = 0;
-  int exp = 0;
 
-  if (f==0.0) goto zero;	/*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;}*/
+  if (f==0.0) {exp = 0; goto zero;}	/*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;}*/
   if (f < 0.0) {f = -f;a[ch++]='-';}
   else if (f > 0.0) ;
   else goto funny;
@@ -99,22 +148,34 @@ static sizet idbl2str(f, a)
     if (ch==0) a[ch++]='+';
   funny: a[ch++]='#'; a[ch++]='.'; a[ch++]='#'; return ch;
   }
-# ifdef DBL_MIN_10_EXP		/* Prevent unnormalized values, as from
-			make-uniform-vector, from causing infinite loops. */
-  while (f < 1.0) {f *= 10.0;  if (exp-- < DBL_MIN_10_EXP) goto funny;}
-  while (f > 10.0) {f *= 0.10; if (exp++ > DBL_MAX_10_EXP) goto funny;}
+  exp = apx_log10(f);
+  f = lpow10(f, -exp);
+  fprec = lpow10(fprec, -exp);
+# ifdef DBL_MIN_10_EXP	/* Prevent random unnormalized values, as from
+			   make-uniform-vector, from causing infinite loops,
+			   but try to print gradually underflowing numbers. */
+  while (f < 1.0) {
+    f *= 10.0;
+    fprec *= 10.0;
+    if (exp-- < DBL_MIN_10_EXP - DBL_DIG - 1) goto funny;
+  }
+  while (f > 10.0) {
+    f /= 10.0;
+    fprec /= 10.0;
+    if (exp++ > DBL_MAX_10_EXP) goto funny;}
 # else
-  while (f < 1.0) {f *= 10.0; exp--;}
-  while (f > 10.0) {f /= 10.0; exp++;}
+  while (f < 1.0) {f *= 10.0; fprec *= 10.0; exp--;}
+  while (f > 10.0) {f /= 10.0; fprec /= 10.0; exp++;}
 # endif
-  if (f+fx[wp] >= 10.0) {f = 1.0; exp++;}
+  fprec *= 0.5;
+  if (f+fprec >= 10.0) {f = 1.0; exp++;}
  zero:
 # ifdef ENGNOT
   dpt = (exp+9999)%3;
   exp -= dpt++;
   efmt = 1;
 # else
-  efmt = (exp < -3) || (exp > wp+2);
+  efmt = (exp < -3) || (exp > dblprec+2);
   if (!efmt)
     if (exp < 0) {
       a[ch++] = '0';
@@ -127,18 +188,20 @@ static sizet idbl2str(f, a)
     dpt = 1;
 # endif
 
-  do {
+  for (i = 30; i--;) {
+    /* printf("  f = %.20g, fprec = %.20g, i = %d\n", f, fprec, i); */
     d = f;
     f -= d;
     a[ch++] = d+'0';
-    if (f < fx[wp])  break;
-    if (f+fx[wp] >= 1.0) {
+    if (f < fprec && f < DBL2STR_FUZZ*fprec) break;
+    if ((f + fprec) >= 1.0 && (f + DBL2STR_FUZZ*fprec) >= 1.0) {
       a[ch-1]++;
       break;
     }
     f *= 10.0;
+    fprec *= 10.0;
     if (!(--dpt))  a[ch++] = '.';
-  } while (wp--);
+  }
 
   if (dpt > 0)
 # ifndef ENGNOT
@@ -565,7 +628,9 @@ SCM istr2flo(str, len, radix)
       switch (c = str[i]) {
       case DIGITS:
 	expon = expon*10 + c-'0';
-	if (expon > MAXEXP)  return BOOL_F; /* exponent too large */
+	if (expon > MAXEXP)
+	  if (1==expsgn || expon > (MAXEXP + dblprec + 1))
+	    return BOOL_F; /* exponent too large */
 	break;
       default:
 	goto out4;
@@ -682,37 +747,26 @@ SCM makdbl (x, y)
 {
   SCM z;
   if ((y==0.0) && (x==0.0)) return flo0;
+  DEFER_INTS;
   if (y==0.0) {
 # ifdef SINGLES
-    float fx;
+    float fx = x;        /* David Yeh <theyeh@uclink.berkeley.edu>
+			    changed this so that MSVC works */
 #  ifndef SINGLESONLY
-    if ((-FLTMAX < x) && (x < FLTMAX) && ((fx=x)==x))
+    if ((-FLTMAX < x) && (x < FLTMAX) && ( (double)fx == x) )
 #  endif
       {
 	NEWCELL(z);
-	DEFER_INTS;
 	CAR(z) = tc_flo;
 	FLO(z) = x;
 	ALLOW_INTS;
 	return z;
       }
 # endif				/* def SINGLES */
-    DEFER_INTS;
-# ifdef NUM_HP
-    CDR(z) = (SCM)num_hp_alloc(sizeof(double));
-# else
-    z = must_malloc_cell(1L*sizeof(double), "real");
-# endif
-    CAR(z) = tc_dblr;
+    z = must_malloc_cell(1L*sizeof(double), (SCM)tc_dblr, "real");
   }
   else {
-    DEFER_INTS;
-# ifdef NUM_HP
-    CDR(z) = (SCM)num_hp_alloc(2L*sizeof(double));
-# else
-    z = must_malloc_cell(2L*sizeof(double), "complex");
-# endif
-    CAR(z) = tc_dblc;
+    z = must_malloc_cell(2L*sizeof(double), (SCM)tc_dblc, "complex");
     IMAG(z) = y;
   }
   REAL(z) = x;
@@ -844,23 +898,23 @@ static SCM vector_equal(x, y)
     if FALSEP(equal(VELTS(x)[i], VELTS(y)[i])) return BOOL_F;
   return BOOL_T;
 }
+#ifdef BIGDIG
 SCM bigequal(x, y)
      SCM x, y;
 {
-#ifdef BIGDIG
   if (0==bigcomp(x, y)) return BOOL_T;
-#endif
   return BOOL_F;
 }
+#endif
+#ifdef FLOATS
 SCM floequal(x, y)
      SCM x, y;
 {
-#ifdef FLOATS
   if (REALPART(x) != REALPART(y)) return BOOL_F;
   if (!(CPLXP(x) && (IMAG(x) != IMAG(y)))) return BOOL_T;
-#endif
   return BOOL_F;
 }
+#endif
 SCM equal(x, y)
      SCM x, y;
 {
@@ -887,7 +941,8 @@ SCM equal(x, y)
 	        if (smobs[i].equalp) return (smobs[i].equalp)(x, y);
 		else return BOOL_F;
 	      }
-	case tc7_bvect: case tc7_uvect: case tc7_ivect:
+	case tc7_bvect:
+	case tc7_uvect: case tc7_ivect: case tc7_svect:
 	case tc7_fvect:	case tc7_cvect: case tc7_dvect: {
 	        SCM (*pred)() = smobs[0x0ff & (tc16_array>>8)].equalp;
 	        if (pred) return (*pred)(x, y);
@@ -1615,9 +1670,8 @@ SCM product(x, y)
       if BIGP(y) return mulbig(BDIGITS(x), NUMDIGS(x), BDIGITS(y), NUMDIGS(y),
 			       BIGSIGN(x) ^ BIGSIGN(y));
       ASRTGO(INEXP(y), bady);
-    bigreal: {
-      double bg = big2dbl(x);
-      return makdbl(bg*REALPART(y), CPLXP(y)?bg*IMAG(y):0.0); }
+    bigreal:
+      return bigdblop('*', x, REALPART(y), CPLXP(y) ? IMAG(y) : 0.0);
     }
     ASRTGO(INEXP(x), badx);
 # else
@@ -1736,7 +1790,6 @@ SCM product(x, y)
     return y;
   }
 }
-
 SCM divide(x, y)
      SCM x, y;
 {
@@ -1767,7 +1820,7 @@ SCM divide(x, y)
         if (z < BIGRAD) {
           SCM w = copybig(x, BIGSIGN(x) ? (y>0) : (y<0));
           return divbigdig(BDIGITS(w), NUMDIGS(w), (BIGDIG)z) ?
-	    makdbl(big2dbl(x)/INUM(y), 0.0) : normbig(w);
+	    bigdblop('/', x, INUM(y), 0.0) : normbig(w);
 	}
 #  ifndef DIGSTOOBIG
         z = pseudolong(z);
@@ -1779,25 +1832,23 @@ SCM divide(x, y)
 	  z = divbigbig(BDIGITS(x), NUMDIGS(x), zdigs, DIGSPERLONG,
 			BIGSIGN(x) ? (y>0) : (y<0), 3);}
 #  endif
-        return z ? z : makdbl(big2dbl(x)/INUM(y), 0.0);
+        return z ? z : bigdblop('/', x, INUM(y), 0.0);
       }
       ASRTGO(NIMP(y), bady);
       if BIGP(y) {
 	z = divbigbig(BDIGITS(x), NUMDIGS(x), BDIGITS(y), NUMDIGS(y),
 		      BIGSIGN(x) ^ BIGSIGN(y), 3);
-	return z ? z : makdbl(big2dbl(x)/big2dbl(y), 0.0);
+	return z ? z : inex_divbigbig(x, y);
       }
       ASRTGO(INEXP(y), bady);
-      if REALP(y) return makdbl(big2dbl(x)/REALPART(y), 0.0);
-      a = big2dbl(x);
-      goto complex_div;
+      return bigdblop('/', x, REALPART(y), CPLXP(y) ? IMAG(y) : 0.0);
     }
 # endif
     ASRTGO(INEXP(x), badx);
     if INUMP(y) {d = INUM(y); goto basic_div;}
 # ifdef BIGDIG
     ASRTGO(NIMP(y), bady);
-    if BIGP(y) {d = big2dbl(y); goto basic_div;}
+    if BIGP(y) return bigdblop('\\', y, REALPART(x), CPLXP(x) ? IMAG(x) : 0.0);
     ASRTGO(INEXP(y), bady);
 # else
     ASRTGO(NIMP(y) && INEXP(y), bady);
@@ -1818,7 +1869,7 @@ SCM divide(x, y)
   if NINUMP(y) {
 # ifdef BIGDIG
     ASRTGO(NIMP(y), bady);
-    if BIGP(y) return makdbl(INUM(x)/big2dbl(y), 0.0);
+    if BIGP(y) return bigdblop('\\', y, INUM(x), 0.0);
 #  ifndef RECKLESS
     if (!(INEXP(y)))
     bady: wta(y, (char *)ARG2, s_divide);
@@ -1905,6 +1956,79 @@ SCM divide(x, y)
   }
 }
 
+SCM scm_intexpt(z1, z2)
+     SCM z1, z2;
+{
+  SCM acc = MAKINUM(1L);
+#ifdef BIGDIG
+  if (INUM0==z1 || acc==z1) return z1;
+  else if (MAKINUM(-1L)==z1) return BOOL_F==evenp(z2)?z1:acc;
+#endif
+  ASSERT(INUMP(z2), z2, ARG2, s_intexpt);
+  z2 = INUM(z2);
+  if (z2 < 0) {
+    z2 = -z2;
+    z1 = divide(z1, UNDEFINED);
+  }
+  if INUMP(z1) {
+    long tmp, iacc = 1, iz1 = INUM(z1);
+    while (1) {
+      if (0==z2) {
+	acc = long2num(iacc);
+#ifndef RECKLESS
+	if (FALSEP(z1))
+	  errout: wta(UNDEFINED, (char *)OVFLOW, s_intexpt);
+#endif
+	return acc;
+      }
+      if (1==z2) {
+	tmp = iacc*iz1;
+	if (tmp/iacc != iz1) {
+	overflow:
+	  z1 = long2num(iz1);
+	  acc = long2num(iacc);
+	  ASRTGO(NFALSEP(z1) && NFALSEP(acc), errout);
+	  goto gencase;
+	}
+	acc = long2num(tmp);
+	ASRTGO(NFALSEP(acc), errout);
+	return acc;
+      }
+      if (z2 & 1) {
+	tmp = iacc*iz1;
+	if (tmp/iacc != iz1) goto overflow;
+	iacc = tmp;
+	z2 = z2 - 1;		/* so jumping to gencase works  */
+      }
+      tmp = iz1*iz1;
+      if (tmp/iz1 != iz1) goto overflow;
+      iz1 = tmp;
+      z2 >>= 1;
+    }
+  }
+  ASSERT(NIMP(z1), z1, ARG1, s_intexpt);
+#ifdef FLOATS			/* Move to scl.c ? */
+  if REALP(z1) {
+    double dacc = 1.0, dz1 = REALPART(z1);
+    while(1) {
+      if (0==z2) return makdbl(dacc, 0.0);
+      if (1==z2) return makdbl(dacc*dz1, 0.0);
+      if (z2 & 1) dacc = dacc*dz1;
+      dz1 = dz1*dz1;
+      z2 >>= 1;
+    }
+  }
+#endif
+ gencase:
+  while(1) {
+    if (0==z2) return acc;
+    if (1==z2) return product(acc, z1);
+    if (z2 & 1) acc = product(acc, z1);
+    z1 = product(z1, z1);
+    z2 >>= 1;
+  }
+}
+
 #ifdef FLOATS
 double lasinh(x)
      double x;
@@ -2172,6 +2296,81 @@ double big2dbl(b)
   if (tc16_bigneg==TYP16(b)) return -ans;
   return ans;
 }
+static double big2scaldbl(b, expt)
+     SCM b;
+     int expt;
+{
+  double ans = 0.0;
+  int i = NUMDIGS(b) - 1;
+  BIGDIG *digits = BDIGITS(b);
+  while (i > (expt/BITSPERDIG)) {
+    ans = digits[i] + BIGRAD*ans;
+    i--;
+  }
+  ans = ldexp(ans, BITSPERDIG - expt);
+  /*
+  if (expt = (expt % BITSPERDIG)) {
+    ans = (digits[i] >> expt) +
+      (1L << (BITSPERDIG - expt))*ans;
+  }
+  if ((1L << (BITSPERDIG - expt - 1)) & digits[i])
+    ans += 1;
+  */
+  if (tc16_bigneg==TYP16(b)) return -ans;
+  return ans;
+}
+static SCM bigdblop(op, b, re, im)
+     int op;
+     SCM b;
+     double re, im;
+{
+  double bm = 0.0;
+  int i = 0;
+  if (NUMDIGS(b)*BITSPERDIG < DBL_MAX_EXP) {
+    bm = big2dbl(b);
+  }
+  else {
+    i = INUM(scm_intlength(b));
+    if (i < DBL_MAX_EXP) {
+      i = 0;
+      bm = big2dbl(b);
+    }
+    else {
+      i = i + 1 - DBL_MAX_EXP;
+      bm = big2scaldbl(b, i);
+    }
+  }
+  switch (op) {
+  case '*':
+    return makdbl(ldexp(bm*re, i), 0.0==im ? 0.0 : ldexp(bm*im, i));
+  case '/': {
+    double d = re*re + im*im;
+    return makdbl(ldexp(bm*(re/d), i), ldexp(-bm*(im/d), i));
+  }
+  case '\\':
+    return makdbl(ldexp(re/bm, -i), 0.0==im ? 0.0 : ldexp(im/bm, -i));
+  }
+}
+static SCM inex_divbigbig(a, b)
+     SCM a, b;
+{
+  double r;
+  if ((NUMDIGS(a)*BITSPERDIG < DBL_MAX_EXP) &&
+      (NUMDIGS(b)*BITSPERDIG < DBL_MAX_EXP))
+    r = big2dbl(a) / big2dbl(b);
+  else {
+    int i = INUM(scm_intlength(a));
+    int j = INUM(scm_intlength(b));
+    i = (i > j) ? i : j;
+    if (i < DBL_MAX_EXP)
+      r = big2dbl(a) / big2dbl(b);
+    else {
+      i = i + 1 - DBL_MAX_EXP;
+      r = big2scaldbl(a, i) / big2scaldbl(b, i);
+    }
+  }
+  return makdbl(r, 0.0);
+}
 # endif
 #endif
 
@@ -2333,6 +2532,7 @@ static iproc subr2s[] = {
 	{s_memv, memv},
 	{s_assv, assv},
 #endif
+	{s_intexpt, scm_intexpt},
 	{s_list_tail, list_tail},
 	{s_ve_fill, vector_fill},
 	{s_st_fill, string_fill},
@@ -2386,13 +2586,13 @@ static dblproc cxrs[] = {
 #endif
 
 #ifdef FLOATS
-# ifndef DBL_DIG
+/* # ifndef DBL_DIG -- also needed for ifndef DBL_MANT_DIG */
 static void add1(f, fsum)
      double f, *fsum;
 {
   *fsum = f + 1.0;
 }
-# endif
+/* #endif */
 #endif
 
 void init_scl()
@@ -2413,9 +2613,8 @@ void init_scl()
   FLO(flo0) = 0.0;
 # else
   DEFER_INTS;
-  flo0 = must_malloc_cell(1L*sizeof(double), "real");
+  flo0 = must_malloc_cell(1L*sizeof(double), (SCM)tc_dblr, "real");
   REAL(flo0) = 0.0;
-  CAR(flo0) = tc_dblr;
   ALLOW_INTS;
 # endif
 # ifdef DBL_DIG
@@ -2432,5 +2631,19 @@ void init_scl()
     dblprec = dblprec-1;
   }
 # endif /* DBL_DIG */
+# ifdef DBL_MANT_DIG
+  dbl_mant_dig = DBL_MANT_DIG;
+# else
+  {
+    double fsum = 0.0, eps = 1.0;
+    int i = 0;
+    while (fsum != 1.0) {
+      eps = 0.5 * eps;
+      add1(eps, &fsum);
+      i++;
+    }
+    dbl_mant_dig = i;
+  }
+# endif /* DBL_MANT_DIG */
 #endif
 }