Import Upstream version 5e3upstream/5e3

1 files changed, 80 insertions, 35 deletions

diff --git a/scl.c b/scl.c
index 9fe779f..7f17cc6 100644
--- a/scl.c
+++ b/scl.c
@@ -64,7 +64,7 @@ static char s_makrect[] = "make-rectangular", s_makpolar[] = "make-polar",
 	    s_real_part[] = "real-part", s_imag_part[] = "imag-part",
 	    s_in2ex[] = "inexact->exact",s_ex2in[] = "exact->inexact";
 
-static char s_expt[] = "$expt", s_atan2[] = "$atan2";
+static char s_expt[] = "real-expt", s_atan2[] = "$atan2";
 #endif
 static char s_memv[] = "memv", s_assv[] = "assv";
 
@@ -585,11 +585,9 @@ SCM istr2flo(str, len, radix)
   }
   if (i==len) return BOOL_F;	/* bad if lone `+' or `-' */
 
-# ifdef FLOATS
   if (6==len && ('+'==str[0] || '-'==str[0]))
     if (0==strcmp(str_inf0, &str[1]))
       return makdbl(1./0. * ('+'==str[0] ? 1 : -1), 0.0);
-# endif
 
   if (str[i]=='i' || str[i]=='I') { /* handle `+i' and `-i'   */
     if (lead_sgn==0.0) return BOOL_F; /* must have leading sign */
@@ -1029,9 +1027,9 @@ 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_svect:
-  case tc7_fvect:	case tc7_cvect: case tc7_dvect: {
+  case tc7_Vbool: case tc7_VfixN8: case tc7_VfixZ8:
+  case tc7_VfixN16: case tc7_VfixZ16: case tc7_VfixN32: case tc7_VfixZ32:
+  case tc7_VfloR32: case tc7_VfloC32: case tc7_VfloC64: case tc7_VfloR64: {
     SCM (*pred)() = smobs[0x0ff & (tc16_array>>8)].equalp;
     if (pred) return (*pred)(x, y);
     else return BOOL_F;
@@ -1942,11 +1940,15 @@ SCM product(x, y)
     return y;
   }
 }
+  /* Use "Smith's formula" to extend dynamic range */
+  /* David Goldberg
+     What Every Computer Scientist Should Know About Floating-Point Arithmetic
+     http://cch.loria.fr/documentation/IEEE754/ACM/goldberg.pdf */
 SCM divide(x, y)
      SCM x, y;
 {
 #ifdef FLOATS
-  double d, r, i, a;
+  double den, a = 1.0;
   if (NINUMP(x)) {
 # ifndef RECKLESS
     if (!(NIMP(x)))
@@ -1956,10 +1958,14 @@ SCM divide(x, y)
 # ifdef BIGDIG
       if (BIGP(x)) return makdbl(1.0/big2dbl(x), 0.0);
 # endif
+      /* reciprocal */
       ASRTGO(INEXP(x), badx);
       if (REALP(x)) return makdbl(1.0/REALPART(x), 0.0);
-      r = REAL(x);  i = IMAG(x);  d = r*r+i*i;
-      return makdbl(r/d, -i/d);
+      {
+	y = x;
+	a = 1.0;
+	goto real_over_complex;
+      }
     }
 # ifdef BIGDIG
     if (BIGP(x)) {
@@ -1999,7 +2005,7 @@ SCM divide(x, y)
     }
 # endif
     ASRTGO(INEXP(x), badx);
-    if (INUMP(y)) {d = INUM(y); goto basic_div;}
+    if (INUMP(y)) {den = INUM(y); goto basic_div;}
 # ifdef BIGDIG
     ASRTGO(NIMP(y), bady);
     if (BIGP(y)) return bigdblop('\\', y, REALPART(x), CPLXP(x) ? IMAG(x) : 0.0);
@@ -2008,13 +2014,28 @@ SCM divide(x, y)
     ASRTGO(NIMP(y) && INEXP(y), bady);
 # endif
     if (REALP(y)) {
-      d = REALPART(y);
-    basic_div: return makdbl(REALPART(x)/d, CPLXP(x)?IMAG(x)/d:0.0);
+      den = REALPART(y);
+    basic_div: return makdbl(REALPART(x)/den, CPLXP(x)?IMAG(x)/den:0.0);
     }
     a = REALPART(x);
-    if (REALP(x)) goto complex_div;
-    r = REAL(y);  i = IMAG(y);  d = r*r+i*i;
-    return makdbl((a*r+IMAG(x)*i)/d, (IMAG(x)*r-a*i)/d);
+    if (REALP(x)) goto real_over_complex;
+    /* Both x and y are complex */
+    /* Use "Smith's formula" to extend dynamic range */
+    {
+      double b = IMAG(x);
+      double c = REALPART(y);
+      double d = IMAG(y);
+      if ((d > 0 ? d : -d) < (c > 0 ? c : -c)) {
+	double r = d/c;
+	double i = c + d*r;
+	return makdbl((a + b*r)/i, (b - a*r)/i);
+      }
+      {
+	double r = c/d;
+	double i = d + c*r;
+	return makdbl((b + a*r)/i, (-a + b*r)/i);
+      }
+    }
   }
   if (UNBNDP(y)) {
     if ((MAKINUM(1L)==x) || (MAKINUM(-1L)==x)) return x;
@@ -2036,9 +2057,23 @@ SCM divide(x, y)
 # endif
     if (REALP(y)) return makdbl(INUM(x)/REALPART(y), 0.0);
     a = INUM(x);
-  complex_div:
-    r = REAL(y);  i = IMAG(y);  d = r*r+i*i;
-    return makdbl((a*r)/d, (-a*i)/d);
+  real_over_complex:
+    /* Both x and y are complex */
+    /* Use "Smith's formula" to extend dynamic range */
+    {
+      double c = REALPART(y);
+      double d = IMAG(y);
+      if ((d > 0 ? d : -d) < (c > 0 ? c : -c)) {
+	double r = d/c;
+	double i = c + d*r;
+	return makdbl((a)/i, (- a*r)/i);
+      }
+      {
+	double r = c/d;
+	double i = d + c*r;
+	return makdbl((a*r)/i, (-a)/i);
+      }
+    }
   }
 #else
 # ifdef BIGDIG
@@ -2370,7 +2405,17 @@ SCM scm_magnitude(z)
   if (CPLXP(z))
     {
       double i = IMAG(z), r = REAL(z);
-      return makdbl(sqrt(i*i+r*r), 0.0);
+      if (i < 0) i = -i;
+      if (r < 0) r = -r;
+      if (i < r) {
+	double q = i / r;
+	return makdbl(r * sqrt(1 + q * q), 0.0);
+      }
+      if (0.0==i) return i;
+      {
+	double q = r / i;
+	return makdbl(i * sqrt(1 + q * q), 0.0);
+      }
     }
   return makdbl(fabs(REALPART(z)), 0.0);
 }
@@ -2848,23 +2893,23 @@ static dblproc cxrs[] = {
 	{"ceiling", ceil},
 	{"truncate", scm_truncate},
 	{"round", scm_round},
-	{"$sqrt", sqrt},
 	{"$abs", fabs},
-	{"$exp", exp},
-	{"$log", log},
-	{"$log10", log10},
-	{"$sin", sin},
-	{"$cos", cos},
-	{"$tan", tan},
-	{"$asin", asin},
-	{"$acos", acos},
-	{"$atan", atan},
-	{"$sinh", sinh},
-	{"$cosh", cosh},
-	{"$tanh", tanh},
-	{"$asinh", asinh},
-	{"$acosh", acosh},
-	{"$atanh", atanh},
+	{"real-sqrt", sqrt},
+	{"real-exp", exp},
+	{"real-ln", log},
+	{"real-log10", log10},
+	{"real-sin", sin},
+	{"real-cos", cos},
+	{"real-tan", tan},
+	{"real-asin", asin},
+	{"real-acos", acos},
+	{"real-atan", atan},
+	{"real-sinh", sinh},
+	{"real-cosh", cosh},
+	{"real-tanh", tanh},
+	{"real-asinh", asinh},
+	{"real-acosh", acosh},
+	{"real-atanh", atanh},
 	{0, 0}};
 #endif