summaryrefslogtreecommitdiffstats
path: root/scl.c
diff options
context:
space:
mode:
Diffstat (limited to 'scl.c')
-rw-r--r--scl.c115
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