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