diff options
Diffstat (limited to 'scl.c')
| -rw-r--r-- | scl.c | 115 |
1 files changed, 80 insertions, 35 deletions
@@ -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 |