From deda2c0fd8689349fea2a900199a76ff7ecb319e Mon Sep 17 00:00:00 2001 From: Bryan Newbold Date: Mon, 20 Feb 2017 00:05:26 -0800 Subject: Import Upstream version 5d6 --- scl.c | 251 +++++++++++++++++++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 212 insertions(+), 39 deletions(-) (limited to 'scl.c') 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 +# ifndef PLAN9 +# include +# 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 } -- cgit v1.2.3