From fd5f104f287427fee885583bc398c137674e6af0 Mon Sep 17 00:00:00 2001 From: Bryan Newbold Date: Fri, 3 Mar 2017 00:56:40 -0800 Subject: New upstream version 5f2 --- scl.c | 898 +++++++++++++++++++++++++++++++++++++++++++----------------------- 1 file changed, 588 insertions(+), 310 deletions(-) (limited to 'scl.c') diff --git a/scl.c b/scl.c index 2fb2d17..d7fbf0d 100644 --- a/scl.c +++ b/scl.c @@ -1,5 +1,5 @@ /* "scl.c" non-IEEE utility functions and non-integer arithmetic. - * Copyright (C) 1990, 1991, 1992, 1993, 1994, 1995, 1997, 2005, 2006 Free Software Foundation, Inc. + * Copyright (C) 1990, 1991, 1992, 1993, 1994, 1995, 1997, 2005, 2006, 2013, 2014 Free Software Foundation, Inc. * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as @@ -25,15 +25,21 @@ # 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 SCM inex_divintbig P((SCM a, SCM b)); +# ifndef BIGDIG static int apx_log10 P((double x)); static double lpow10 P((double x, int n)); +# endif static sizet idbl2str P((double f, char *a)); +static sizet pdbl2str P((double f, char *a, sizet ch)); static sizet iflo2str P((SCM flt, char *str)); static void safe_add_1 P((double f, double *fsum)); static long scm_twos_power P((SCM n)); +static double mantexp2dbl P((SCM manstr, int expo)); +static double pmantexp2dbl P((SCM bmant, int expo)); +static SCM ex2in P((SCM z)); +static SCM ilog P((unsigned long m, SCM b, SCM k, unsigned long *n)); static char s_makrect[] = "make-rectangular", s_makpolar[] = "make-polar", s_magnitude[] = "magnitude", s_angle[] = "angle", @@ -63,29 +69,171 @@ 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"; -static char str_inf0[] = "inf.0"; +static char str_inf0[] = "inf.0", str_nan0[] = "nan.0", str_0[] = "0.0"; +static char s_intexpt[] = "integer-expt", s_cintlog[] = "ceiling-integer-log"; +static char s_dfloat_parts[] = "double-float-parts"; +#define s_intlog (&s_cintlog[8]) /*** NUMBERS -> STRINGS ***/ #ifdef FLOATS 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; + +int inf2str(f, a) + double f; + char *a; +{ + sizet ch = 0; + if (f < 0.0) a[ch++] = '-'; + else if (f > 0.0) a[ch++] = '+'; + else { + a[ch++] = '+'; + strcpy(&a[ch], str_nan0); + return ch+sizeof(str_nan0)-1; + } + strcpy(&a[ch], str_inf0); + return ch+sizeof(str_inf0)-1; +} + +/* idbl2str() is the top level double-to-string conversion */ +static sizet idbl2str(f, a) + double f; + char *a; +{ + sizet ch = 0; + if (f==0.0) {strcpy(a, str_0); return sizeof(str_0)-1;} + if (f==2*f) return inf2str(f, a); + if (f < 0.0) {f = -f;a[ch++]='-';} + else if (f > 0.0) ; + else return inf2str(f, a); + return pdbl2str(f, a, ch); +} + +static double llog2 = .30102999566398114; /* log10(2) */ + +/* There are also temporary strings used in number->string conversion. */ +void strrecy(str) + SCM str; +{ + if (IMP(str) || !STRINGP(str)) return; + DEFER_INTS; + must_free(CHARS(str), (sizet)LENGTH(str)); + CAR(str) = MAKINUM(1); + CDR(str) = INUM0; + ALLOW_INTS; +} + +# ifdef BIGDIG + +/* The useful extent of bignums used in floating-point conversions is */ +/* limited. Recycle them explicitly, rather than waiting for GC. */ + +void bigrecy(bgnm) + SCM bgnm; +{ + if (IMP(bgnm) || !BIGP(bgnm)) return; + DEFER_INTS; + must_free(CHARS(bgnm), (sizet)NUMDIGS(bgnm)*sizeof(BIGDIG)); + CAR(bgnm) = INUM0; + CDR(bgnm) = INUM0; + ALLOW_INTS; +} + +/* can convert to string accurately with bignums */ +/* f > 0 */ +/* DBL_MIN_EXP = -1021 */ +/* dbl_mant_dig = 53 */ +static sizet pdbl2str(f, a, ch) + double f; + char *a; + sizet ch; +{ + SCM mant, quo, num; + sizet dp = ch; + int e2, point, ndig = dbl_mant_dig; + /* code from scm_dfloat_parts() */ + double dman = frexp(f, &e2); +# ifdef DBL_MIN_EXP + if (e2 < DBL_MIN_EXP) + ndig -= DBL_MIN_EXP - e2; +# endif + e2 -= ndig; + mant = dbl2big(ldexp(dman, ndig)); + point = ceil(e2*llog2); + /* if (scm_verbose > 1) printf("mantissa = %g -> #x%s; e2 = %d -> %d; point = %d; ndig = %d -> %d\n", dman, CHARS(number2string(mant, MAKINUM(16))), e2+ndig, e2, point, dbl_mant_dig, ndig); */ + if (e2 >= 0) { + /* try first with starved precision */ + { + num = scm_ash(mant, MAKINUM(e2 - point)); + bigrecy(mant); + quo = scm_round_quotient(num, VELTS(pows5)[(long) point]); + if (pmantexp2dbl(quo, point) != f) { + bigrecy(quo); quo = num; + num = scm_ash(quo, MAKINUM(1L)); + bigrecy(quo); + quo = scm_round_quotient(num, VELTS(pows5)[(long) --point]); + } + } + } else { /* e2 <= 0 */ + /* try first with starved precision */ + { + SCM den = scm_ash(MAKINUM(1L), MAKINUM(point - e2)); + num = product(mant, VELTS(pows5)[- (long) point]); + bigrecy(mant); + quo = scm_round_quotient(num, den); + if (pmantexp2dbl(quo, point) != f) { + bigrecy(quo); quo = num; + point--; + num = product(quo, MAKINUM(10)); + if (mant != MAKINUM(1)) bigrecy(quo); + quo = scm_round_quotient(num, den); + } + bigrecy(den); + } + } + bigrecy(num); + a[ch++] = '.'; + /* if (sizeof(UBIGLONG)>=sizeof(double)) /\* Is ulong larger than mantissa? *\/ */ + /* ch += iulong2str(num2ulong(quo, (char *)ARG1, s_number2string), 10, &a[ch]); */ + /* else */ { + SCM str = number2string(quo, MAKINUM(10)); + int len = LENGTH(str), i = 0; + bigrecy(quo); + point += len - 1; + while (i < len) a[ch++] = CHARS(str)[i++]; + strrecy(str); + } + a[dp] = a[dp+1]; a[++dp] = '.'; +# ifdef ENGNOT + while ((dp+1 < ch) && (point+9999)%3) { a[dp] = a[dp+1]; a[++dp] = '.'; point--; } +# endif /* ENGNOT */ + while ('0'==a[--ch]); ch++; + if (point != 0) { + a[ch++] = 'e'; + return ch + ilong2str(point, 10, &a[ch]); + } else return ch; +} +# else /* ~BIGDIG */ + +/* DBL2STR_FUZZ is a somewhat arbitrary guard against + round off error in scaling f and fprec. */ +# define DBL2STR_FUZZ 0.9 +int dblprec; +/* static double dbl_eps; */ double dbl_prec(x) double x; { int expt; double frac = frexp(x, &expt); -# ifdef DBL_MIN_EXP +# 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 +# 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; { @@ -116,51 +264,17 @@ static double lpow10(x, n) return x/p10[-n]; } -int inf2str(f, a) - double f; - char *a; -{ - sizet ch = 0; - if (f < 0.0) a[ch++] = '-'; - else if (f > 0.0) a[ch++] = '+'; - else { - a[ch++] = '0'; a[ch++] = '/'; a[ch++] = '0'; - return ch; - } - while (str_inf0[ch - 1]) { - a[ch] = str_inf0[ch - 1]; - ch++; - } -/* # ifdef COMPACT_INFINITY_NOTATION */ -/* else a[ch++] = '0'; */ -/* # else */ -/* a[ch++] = (f != f) ? '0' : '1'; */ -/* # endif */ -/* a[ch++] = '/'; a[ch++] = '0'; */ - return ch; -} - -/* 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) +/* f is finite and positive */ +static sizet pdbl2str(f, a, ch) double f; char *a; + sizet ch; { double fprec = dbl_prec(f); - int efmt, dpt, d, i, exp; - sizet ch = 0; - - if (f==0.0) {exp = 0; goto zero;} /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;}*/ - if (f==2*f) return inf2str(f, a); - if (f < 0.0) {f = -f;a[ch++]='-';} - else if (f > 0.0) ; - else return inf2str(f, a); - exp = apx_log10(f); + int efmt, dpt, d, i, exp = apx_log10(f); f = lpow10(f, -exp); fprec = lpow10(fprec, -exp); -# ifdef DBL_MIN_10_EXP /* Prevent random unnormalized values, as from +# 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) { @@ -173,18 +287,18 @@ static sizet idbl2str(f, a) fprec /= 10.0; if (exp++ > DBL_MAX_10_EXP) return inf2str(f, a); } -# else +# else while (f < 1.0) {f *= 10.0; fprec *= 10.0; exp--;} while (f > 10.0) {f /= 10.0; fprec /= 10.0; exp++;} -# endif +# endif fprec *= 0.5; if (f+fprec >= 10.0) {f = 1.0; exp++;} - zero: -# ifdef ENGNOT + /* zero: */ +# ifdef ENGNOT dpt = (exp+9999)%3; exp -= dpt++; efmt = 1; -# else +# else efmt = (exp < -3) || (exp > dblprec+2); if (!efmt) if (exp < 0) { @@ -196,7 +310,7 @@ static sizet idbl2str(f, a) dpt = exp+1; else dpt = 1; -# endif +# endif for (i = 30; i--;) { /* printf(" f = %.20g, fprec = %.20g, i = %d\n", f, fprec, i); */ @@ -214,7 +328,7 @@ static sizet idbl2str(f, a) } if (dpt > 0) -# ifndef ENGNOT +# ifndef ENGNOT if ((dpt > 4) && (exp > 6)) { d = (a[0]=='-'?2:1); for (i = ch++; i > d; i--) @@ -222,7 +336,7 @@ static sizet idbl2str(f, a) a[d] = '.'; efmt = 1; } else -# endif +# endif { while (--dpt) a[ch++] = '0'; a[ch++] = '.'; @@ -242,6 +356,7 @@ static sizet idbl2str(f, a) } return ch; } +# endif /* ~BIGDIG */ static sizet iflo2str(flt, str) SCM flt; @@ -263,7 +378,7 @@ static sizet iflo2str(flt, str) } #endif /* FLOATS */ -sizet iuint2str(num, rad, p) +sizet iulong2str(num, rad, p) unsigned long num; int rad; char *p; @@ -281,16 +396,16 @@ sizet iuint2str(num, rad, p) } return j; } -sizet iint2str(num, rad, p) +sizet ilong2str(num, rad, p) long num; int rad; char *p; { if ((num < 0) && !(rad < 0)) { *p++ = '-'; - return 1 + iuint2str((unsigned long) -num, rad, p); + return 1 + iulong2str((unsigned long) -num, rad, p); } - return iuint2str((unsigned long) num, rad < 0 ? -rad : rad, p); + return iulong2str((unsigned long) num, rad < 0 ? -rad : rad, p); } #ifdef BIGDIG static SCM big2str(b, radix) @@ -329,6 +444,7 @@ static SCM big2str(b, radix) for (i = j;j < LENGTH(ss);j++) s[ch+j-i] = s[j]; /* jeh */ resizuve(ss, (SCM)MAKINUM(ch+LENGTH(ss)-i)); /* jeh */ } + bigrecy(t); return ss; } #endif @@ -364,7 +480,7 @@ SCM number2string(x, radix) #endif { char num_buf[INTBUFLEN]; - return makfromstr(num_buf, iint2str(INUM(x), (int)INUM(radix), num_buf)); + return makfromstr(num_buf, ilong2str(INUM(x), (int)INUM(radix), num_buf)); } } /* These print routines are stubbed here so that repl.c doesn't need @@ -406,17 +522,17 @@ int bigprint(exp, port, writing) SCM istr2int(str, len, radix) char *str; long len; - register long radix; + register int radix; { - sizet j; register sizet k, blen = 1; - sizet i = 0; + sizet i = 0, j; int c; SCM res; register BIGDIG *ds; - register unsigned long t2; + register UBIGLONG t2; if (0 >= len) return BOOL_F; /* zero length */ + /* Estimate number of digits; will trim during finish */ if (10==radix) j = 1+(84*len)/(BITSPERDIG*25); else j = (8 < radix) ? 1+(4*len)/BITSPERDIG : 1+(3*len)/BITSPERDIG; switch (str[0]) { /* leading sign */ @@ -425,6 +541,7 @@ SCM istr2int(str, len, radix) } res = mkbig(j, '-'==str[0]); ds = BDIGITS(res); + /* clear allocated digits */ for (k = j;k--;) ds[k] = 0; do { switch (c = str[i++]) { @@ -441,9 +558,9 @@ SCM istr2int(str, len, radix) k = 0; t2 = c; moretodo: - while(k < blen) { -/* printf("k = %d, blen = %d, t2 = %ld, ds[k] = %d\n", k, blen, t2, ds[k]);*/ - t2 += ds[k]*radix; + /* Add t2 into temp bignum */ + while (k < blen) { + t2 += ((UBIGLONG)ds[k])*radix; ds[k++] = BIGLO(t2); t2 = BIGDN(t2); } @@ -463,7 +580,7 @@ SCM istr2int(str, len, radix) SCM istr2int(str, len, radix) register char *str; long len; - register long radix; + register int radix; { register long n = 0, ln; register int c; @@ -514,31 +631,114 @@ SCM istr2int(str, len, radix) #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; + +/* [Jaffer] */ +/* In Clinger's "How to Read Floating-Point Numbers Accurately" */ +/* the key issue is that successive rounding does not have the */ +/* same effect as one rounding operation. With bignums it is */ +/* a simple matter to accumulate digits without rounding. */ +/* The approach here is to compute the binary value without rounding, */ +/* then round explicitly when scaling the bigint to fit into the */ +/* mantissa. */ + +/* manstr is the mantissa with the decimal point removed and point + (the exponent) adjusted appropriately */ +double mantexp2dbl(manstr, point) + SCM manstr; + int point; +{ + SCM bmant = istr2int(CHARS(manstr), LENGTH(manstr), 10L); + double val = pmantexp2dbl(bmant, point); + bigrecy(bmant); + return val; +} +double pmantexp2dbl(bmant, point) + SCM bmant; + int point; +{ + double ans; + if (BOOL_F == bmant) return REAL(scm_narn); + if (point >= 0) { + SCM quo, num = product(bmant, VELTS(pows5)[(long) point]); + int bex = INUM(scm_intlength(num)) - dbl_mant_dig; + if (bex > 0) { + SCM den = scm_ash(MAKINUM(1L), MAKINUM(bex)); + quo = scm_round_quotient(num, den); + bigrecy(den); + } + else + quo = scm_ash(num, MAKINUM(- bex)); + /* quo may not be a bignum */ + if (INUMP(quo)) ans = ldexp((double)(INUM(quo)), bex + point); + else { + sizet i = NUMDIGS(quo); + sizet j = i - (dbl_mant_dig + BITSPERDIG - 1)/BITSPERDIG; + BIGDIG *digits = BDIGITS(quo); + ans = 0.0; + while (i-- > j) ans = digits[i] + ldexp(ans, BITSPERDIG); + bex += j * BITSPERDIG; + ans = ldexp(ans, bex + point); + } + if (num != quo) bigrecy(quo); + if (bmant != MAKINUM(1L)) bigrecy(num); + return ans; + } else { + int maxpow = LENGTH(pows5) - 1; + SCM scl = (-point <= maxpow) ? + VELTS(pows5)[(long) -point] : + product(VELTS(pows5)[(long)maxpow], VELTS(pows5)[(long)-point-maxpow]); + int bex = /* bex < 0 */ + INUM(scm_intlength(bmant)) - INUM(scm_intlength(scl)) - dbl_mant_dig; + SCM num = scm_ash(bmant, MAKINUM(-bex)); + SCM quo = scm_round_quotient(num, scl); + if (INUM(scm_intlength(quo)) > dbl_mant_dig) { + bex++; /* too many bits of quotient */ + quo = scm_round_quotient(num, scm_ash(scl, MAKINUM(1L))); + } + if (-point > maxpow) bigrecy(scl); + bigrecy(num); + ans = ldexp(int2dbl(quo), bex + point); + bigrecy(quo); + return ans; } +} + +# else /* def BIGDIG */ + +double mantexp2dbl(manstr, point) + SCM manstr; + int point; +{ + register int c, i = 0; + double res = 0.0; + char *str = CHARS(manstr); + int len = LENGTH(manstr); + do { /* check initial digits */ + switch (c = str[i]) { + case DIGITS: + c = c - '0'; + res = res * 10 + c; + break; + case 'D': case 'E': case 'F': + case 'd': case 'e': case 'f': + default: + goto out1; + } + } while (++i < len); + out1: + if (point >= 0) + while (point--) res *= 10.0; + else + while (point++) { +# ifdef _UNICOS + res *= 0.1; +# else + res /= 10.0; # endif - d = INUM(n); - if (0==d) return 0; - out: - do { - d4 = 15 & d; - c += twostab[d4]; - d >>= 4; - } while (0==d4); - return c; + } + return res; } + # endif /* def BIGDIG */ SCM istr2flo(str, len, radix) @@ -546,37 +746,38 @@ SCM istr2flo(str, len, radix) register long len; register long radix; { - register int c, i = 0; - double lead_sgn = 0.0; + register int c, i = 0, j = 0; + int lead_sgn = 0; double res = 0.0, tmp = 0.0; int flg = 0; int point = 0; - SCM second; + int shrtp = 0; + SCM second, manstr; + char *mant; if (i >= len) return BOOL_F; /* zero length */ switch (*str) { /* leading sign */ - case '-': lead_sgn = -1.0; i++; break; - case '+': lead_sgn = 1.0; i++; break; + case '-': lead_sgn = -1; i++; break; + case '+': lead_sgn = 1; i++; break; } if (i==len) return BOOL_F; /* bad if lone `+' or `-' */ - if (6==len && ('+'==str[0] || '-'==str[0])) - if (0==strcmp(str_inf0, &str[1])) + if (6==len && ('+'==str[0] || '-'==str[0])) { + if (0==strcasecmp(str_inf0, &str[1])) return makdbl(1./0. * ('+'==str[0] ? 1 : -1), 0.0); - + else if (0==strcasecmp(str_nan0, &str[1])) + return scm_narn; + } if (str[i]=='i' || str[i]=='I') { /* handle `+i' and `-i' */ - if (lead_sgn==0.0) return BOOL_F; /* must have leading sign */ + if (lead_sgn==0) return BOOL_F; /* must have leading sign */ if (++i < len) return BOOL_F; /* `i' not last character */ - return makdbl(0.0, lead_sgn); - } - /* # ifdef COMPACT_INFINITY_NOTATION */ - if (0.0 != lead_sgn && str[i]=='/') { - res = 1; - flg = 1; - goto out1; + return makdbl(0.0, (double)lead_sgn); } - /* # endif */ + + manstr = makstr(len); + mant = CHARS(manstr); + do { /* check initial digits */ switch (c = str[i]) { case DIGITS: @@ -593,6 +794,7 @@ SCM istr2flo(str, len, radix) c = c-'a'+10; accum1: if (c >= radix) return BOOL_F; /* bad digit for radix */ + mant[j++] = str[i]; res = res * radix + c; flg = 1; /* res is valid */ break; @@ -607,16 +809,15 @@ SCM istr2flo(str, len, radix) /* By here, must have seen a digit, or must have next char be a `.' with radix==10 */ - if (!flg) - if (!(str[i]=='.' && radix==10)) - return BOOL_F; + if ((!flg) && (!(str[i]=='.' && radix==10))) return BOOL_F; while (str[i]=='#') { /* optional sharps */ res *= radix; + mant[j++] = '0'; if (++i==len) goto done; } - if (str[i]=='/') { + if (str[i]=='/' && i+1 < len) { while (++i < len) { switch (c = str[i]) { case DIGITS: @@ -640,6 +841,7 @@ SCM istr2flo(str, len, radix) if (i < len) while (str[i]=='#') { /* optional sharps */ tmp *= radix; + mant[j++] = '0'; if (++i==len) break; } res /= tmp; @@ -654,6 +856,7 @@ SCM istr2flo(str, len, radix) point--; res = res*10.0 + c-'0'; flg = 1; + mant[j++] = str[i]; break; default: goto out3; @@ -668,60 +871,56 @@ SCM istr2flo(str, len, radix) } switch (str[i]) { /* exponent */ + case 'f': case 'F': + case 's': case 'S': + shrtp = !0; case 'd': case 'D': case 'e': case 'E': - case 'f': case 'F': - case 'l': case 'L': - case 's': case 'S': { - int expsgn = 1, expon = 0; - if (radix != 10) return BOOL_F; /* only in radix 10 */ - if (++i==len) return BOOL_F; /* bad exponent */ - switch (str[i]) { - case '-': expsgn=(-1); - case '+': if (++i==len) return BOOL_F; /* bad exponent */ - } - if (str[i] < '0' || str[i] > '9') return BOOL_F; /* bad exponent */ - do { - switch (c = str[i]) { - case DIGITS: - expon = expon*10 + c-'0'; - /* if (expon > MAXEXP) */ - /* if (1==expsgn || expon > (MAXEXP + dblprec + 1)) */ - /* return BOOL_F; /\* exponent too large *\/ */ - break; - default: - goto out4; + case 'l': case 'L': { + int expsgn = 1, expon = 0; + if (radix != 10) return BOOL_F; /* only in radix 10 */ + if (++i==len) return BOOL_F; /* bad exponent */ + switch (str[i]) { + case '-': expsgn=(-1); + case '+': if (++i==len) return BOOL_F; /* bad exponent */ } - } while (++i < len); - out4: - point += expsgn*expon; - } + if (str[i] < '0' || str[i] > '9') return BOOL_F; /* bad exponent */ + do { + switch (c = str[i]) { + case DIGITS: + expon = expon*10 + c-'0'; + /* if (expon > MAXEXP) */ + /* if (1==expsgn || expon > (MAXEXP + dblprec + 1)) */ + /* return BOOL_F; /\* exponent too large *\/ */ + break; + default: + goto out4; + } + } while (++i < len); + out4: + point += expsgn*expon; + } } adjust: - if (point >= 0) - while (point--) res *= 10.0; - else -# ifdef _UNICOS - while (point++) res *= 0.1; -# else - while (point++) res /= 10.0; -# endif + mant[j] = 0; + manstr = resizuve(manstr, MAKINUM(j)); + if (radix == 10) res = mantexp2dbl(manstr, point); done: /* at this point, we have a legitimate floating point result */ - if (lead_sgn==-1.0) res = -res; - if (i==len) return makdbl(res, 0.0); + if (lead_sgn==-1) res = -res; + if (i==len) return shrtp ? makflo(res) : makdbl(res, 0.0); if (str[i]=='i' || str[i]=='I') { /* pure imaginary number */ - if (lead_sgn==0.0) return BOOL_F; /* must have leading sign */ + if (lead_sgn==0) return BOOL_F; /* must have leading sign */ if (++i < len) return BOOL_F; /* `i' not last character */ return makdbl(0.0, res); } switch (str[i++]) { - case '-': lead_sgn = -1.0; break; - case '+': lead_sgn = 1.0; break; + case '-': lead_sgn = -1; break; + case '+': lead_sgn = 1; break; case '@': { /* polar input for complex number */ /* get a `real' for angle */ second = istr2flo(&str[i], (long)(len-i), radix); @@ -737,7 +936,7 @@ SCM istr2flo(str, len, radix) /* at this point, last char must be `i' */ if (str[len-1] != 'i' && str[len-1] != 'I') return BOOL_F; /* handles `x+i' and `x-i' */ - if (i==(len-1)) return makdbl(res, lead_sgn); + if (i==(len-1)) return makdbl(res, (double)lead_sgn); /* get a `ureal' for complex part */ second = istr2flo(&str[i], (long)((len-i)-1), radix); if (IMP(second)) return BOOL_F; @@ -1003,9 +1202,10 @@ SCM equal(x, y) if (smobs[i].equalp) return (smobs[i].equalp)(x, y); else return BOOL_F; } - 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: { + case tc7_VfixN8: case tc7_VfixZ8: case tc7_VfixN16: case tc7_VfixZ16: + case tc7_VfixN32: case tc7_VfixZ32: case tc7_VfixN64: case tc7_VfixZ64: + case tc7_VfloR32: case tc7_VfloC32: case tc7_VfloC64: case tc7_VfloR64: + case tc7_Vbool: { SCM (*pred)() = smobs[0x0ff & (tc16_array>>8)].equalp; if (pred) return (*pred)(x, y); else return BOOL_F; @@ -1035,21 +1235,50 @@ SCM scm_complex_p(obj) } # 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; +} + int scm_bigdblcomp(b, d) SCM b; double d; { - sizet dlen, blen; + int dlen = 99; + sizet blen = 0; 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 (!(d==2*d && d != 0.0)) { + 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); + double bd = int2dbl(b); if (bd > d) return -1; if (bd < d) return 1; return 0; @@ -1357,7 +1586,7 @@ SCM negativep(x) } static char s_exactprob[] = "not representable as inexact"; -SCM lmax(x, y) +SCM scm_max(x, y) SCM x, y; { #ifdef FLOATS @@ -1382,7 +1611,7 @@ SCM lmax(x, y) ASRTGO(REALP(y), bady); big_dbl: if (-1 != scm_bigdblcomp(x, REALPART(y))) return y; - z = big2dbl(x); + z = int2dbl(x); ASRTER(0==scm_bigdblcomp(x, z), x, s_exactprob, s_max); return makdbl(z, 0.0); } @@ -1441,7 +1670,7 @@ SCM lmax(x, y) return ((long)x < (long)y) ? y : x; } -SCM lmin(x, y) +SCM scm_min(x, y) SCM x, y; { #ifdef FLOATS @@ -1466,7 +1695,7 @@ SCM lmin(x, y) ASRTGO(REALP(y), bady); big_dbl: if (1 != scm_bigdblcomp(x, REALPART(y))) return y; - z = big2dbl(x); + z = int2dbl(x); ASRTER(0==scm_bigdblcomp(x, z), x, s_exactprob, s_min); return makdbl(z, 0.0); } @@ -1549,7 +1778,7 @@ SCM sum(x, y) return addbig(BDIGITS(x), NUMDIGS(x), BIGSIGN(x), y, 0); } ASRTGO(INEXP(y), bady); - bigreal: return makdbl(big2dbl(x)+REALPART(y), CPLXP(y)?IMAG(y):0.0); + bigreal: return makdbl(int2dbl(x)+REALPART(y), CPLXP(y)?IMAG(y):0.0); } ASRTGO(INEXP(x), badx); # else @@ -1647,7 +1876,7 @@ SCM difference(x, y) if (NINUMP(x)) { # ifndef RECKLESS if (!(NIMP(x))) - badx: wta(x, (char *)ARG1, s_difference); + badx: wta(x, (char *)ARG1, s_difference); # endif if (UNBNDP(y)) { # ifdef BIGDIG @@ -1666,12 +1895,12 @@ SCM difference(x, y) if (BIGP(x)) { if (BIGP(y)) return (NUMDIGS(x) < NUMDIGS(y)) ? addbig(BDIGITS(x), NUMDIGS(x), BIGSIGN(x), y, 0x0100) : - addbig(BDIGITS(y), NUMDIGS(y), BIGSIGN(y) ^ 0x0100, x, 0); + addbig(BDIGITS(y), NUMDIGS(y), BIGSIGN(y) ^ 0x0100, x, 0); ASRTGO(INEXP(y), bady); - return makdbl(big2dbl(x)-REALPART(y), CPLXP(y)?-IMAG(y):0.0); + return makdbl(int2dbl(x)-REALPART(y), CPLXP(y)?-IMAG(y):0.0); } ASRTGO(INEXP(x), badx); - if (BIGP(y)) return makdbl(REALPART(x)-big2dbl(y), CPLXP(x)?IMAG(x):0.0); + if (BIGP(y)) return makdbl(REALPART(x)-int2dbl(y), CPLXP(x)?IMAG(x):0.0); ASRTGO(INEXP(y), bady); # else ASRTGO(INEXP(x), badx); @@ -1701,12 +1930,12 @@ SCM difference(x, y) } # ifndef RECKLESS if (!(INEXP(y))) - bady: wta(y, (char *)ARG2, s_difference); + bady: wta(y, (char *)ARG2, s_difference); # endif # else # ifndef RECKLESS if (!(NIMP(y) && INEXP(y))) - bady: wta(y, (char *)ARG2, s_difference); + bady: wta(y, (char *)ARG2, s_difference); # endif # endif return makdbl(INUM(x)-REALPART(y), CPLXP(y)?-IMAG(y):0.0); @@ -1739,7 +1968,7 @@ SCM difference(x, y) if (NINUMP(y)) { # ifndef RECKLESS if (!(NIMP(y) && BIGP(y))) - bady: wta(y, (char *)ARG2, s_difference); + bady: wta(y, (char *)ARG2, s_difference); # endif { # ifndef DIGSTOOBIG @@ -1932,7 +2161,7 @@ SCM divide(x, y) # endif if (UNBNDP(y)) { # ifdef BIGDIG - if (BIGP(x)) return makdbl(1.0/big2dbl(x), 0.0); + if (BIGP(x)) return inex_divintbig(MAKINUM(1L), x); # endif /* reciprocal */ ASRTGO(INEXP(x), badx); @@ -1953,19 +2182,23 @@ SCM divide(x, y) if (z < 0) z = -z; if (z < BIGRAD) { SCM w = copybig(x, BIGSIGN(x) ? (y>0) : (y<0)); - return divbigdig(BDIGITS(w), NUMDIGS(w), (BIGDIG)z) ? - bigdblop('/', x, INUM(y), 0.0) : normbig(w); + int sts = divbigdig(BDIGITS(w), NUMDIGS(w), (BIGDIG)z); + if (sts) { + bigrecy(w); + return bigdblop('/', x, INUM(y), 0.0); + } + else return normbig(w); } # ifndef DIGSTOOBIG z = pseudolong(z); z = divbigbig(BDIGITS(x), NUMDIGS(x), (BIGDIG *)&z, DIGSPERLONG, - BIGSIGN(x) ? (y>0) : (y<0), 3); + BIGSIGN(x) ? (y>0) : (y<0), 4); # else { BIGDIG zdigs[DIGSPERLONG]; longdigs(z, zdigs); z = divbigbig(BDIGITS(x), NUMDIGS(x), zdigs, DIGSPERLONG, - BIGSIGN(x) ? (y>0) : (y<0), 3); + BIGSIGN(x) ? (y>0) : (y<0), 4); } # endif return z ? z : bigdblop('/', x, INUM(y), 0.0); @@ -1973,8 +2206,8 @@ SCM divide(x, y) 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 : inex_divbigbig(x, y); + BIGSIGN(x) ^ BIGSIGN(y), 4); + return z ? z : inex_divintbig(x, y); } ASRTGO(INEXP(y), bady); return bigdblop('/', x, REALPART(y), CPLXP(y) ? IMAG(y) : 0.0); @@ -2020,7 +2253,7 @@ SCM divide(x, y) if (NINUMP(y)) { # ifdef BIGDIG ASRTGO(NIMP(y), bady); - if (BIGP(y)) return bigdblop('\\', y, INUM(x), 0.0); + if (BIGP(y)) return inex_divintbig(x, y); /* bigdblop('\\', y, INUM(x), 0.0); */ # ifndef RECKLESS if (!(INEXP(y))) bady: wta(y, (char *)ARG2, s_divide); @@ -2064,25 +2297,29 @@ SCM divide(x, y) if (z < 0) z = -z; if (z < BIGRAD) { SCM w = copybig(x, BIGSIGN(x) ? (y>0) : (y<0)); - if (divbigdig(BDIGITS(w), NUMDIGS(w), (BIGDIG)z)) goto ov; + int sts = divbigdig(BDIGITS(w), NUMDIGS(w), (BIGDIG)z); + if (sts) { + bigrecy(w); + goto ov; + } return w; } # ifndef DIGSTOOBIG z = pseudolong(z); z = divbigbig(BDIGITS(x), NUMDIGS(x), (BIGDIG *)&z, DIGSPERLONG, - BIGSIGN(x) ? (y>0) : (y<0), 3); + BIGSIGN(x) ? (y>0) : (y<0), 4); # else { BIGDIG zdigs[DIGSPERLONG]; longdigs(z, zdigs); z = divbigbig(BDIGITS(x), NUMDIGS(x), zdigs, DIGSPERLONG, - BIGSIGN(x) ? (y>0) : (y<0), 3); + BIGSIGN(x) ? (y>0) : (y<0), 4); } # endif } else { ASRTGO(NIMP(y) && BIGP(y), bady); z = divbigbig(BDIGITS(x), NUMDIGS(x), BDIGITS(y), NUMDIGS(y), - BIGSIGN(x) ^ BIGSIGN(y), 3); + BIGSIGN(x) ^ BIGSIGN(y), 4); } if (!z) goto ov; return z; @@ -2123,11 +2360,50 @@ SCM divide(x, y) } } +SCM ilog(m, b, k, n) + unsigned long m; + SCM b, k; + unsigned long *n; +{ + /* printf("call ilog %ld ", m); scm_iprin1(b, cur_outp, 0); printf(" "); scm_iprin1(k, cur_outp, 0); printf("\n"); */ + if (BOOL_T==lessp(k, b)) return k; + *n += m; + { + SCM q = ilog(2*m, product(b, b), lquotient(k, b), n); + /* printf("return ilog "); scm_iprin1(q, cur_outp, 0); printf("\n"); */ + if (BOOL_T==lessp(q, b)) return q; + *n += m; + return lquotient(q, b); + } +} + +SCM scm_intlog(base, k) + SCM base, k; +{ + unsigned long n = 1; + ASRTER(INUMP(base) || (NIMP(base) && BIGP(base)), base, ARG1, s_intlog); + ASRTER(BOOL_T==lessp(MAKINUM(1L), base), base, OUTOFRANGE, s_intlog); + ASRTER((INUMP(k) && k > 0) || (NIMP(k) && TYP16(k)==tc16_bigpos), k, ARG2, s_intlog); + if (BOOL_T==lessp(k, base)) return INUM0; + ilog(1, base, lquotient(k, base), &n); + return MAKINUM(n); +} + +#ifdef INUMS_ONLY +# define eqv eqp +#endif +SCM scm_cintlog(base, k) + SCM base, k; +{ + SCM il = scm_intlog(base, k); + return (BOOL_T==eqv(k, scm_intexpt(base, il))) ? il : sum(MAKINUM(1L), il); +} + SCM scm_intexpt(z1, z2) SCM z1, z2; { SCM acc = MAKINUM(1L); - int recip = 0; + long iz2; #ifdef FLOATS double dacc, dz1; #endif @@ -2135,25 +2411,20 @@ SCM scm_intexpt(z1, z2) ASRTER(INUMP(z2), z2, ARG2, s_intexpt); if (acc==z1) return z1; if (MAKINUM(-1L)==z1) return BOOL_F==evenp(z2)?z1:acc; - z2 = INUM(z2); - if (z2 < 0) { - z2 = -z2; - recip = 1; /* z1 = divide(z1, UNDEFINED); */ + iz2 = INUM(z2); + if (iz2 < 0L) { + iz2 = -iz2; + z1 = divide(z1, UNDEFINED); } if (INUMP(z1)) { - long tmp, iacc = 1, iz1 = INUM(z1); -#ifdef FLOATS - if (recip) { dz1 = iz1; goto flocase; } -#endif + long tmp, iacc = 1L, iz1 = INUM(z1); while (1) { - if (0==z2) { + if (0L==iz2) { acc = long2num(iacc); break; } - if (0==iz1) - if (0==recip) return z1; - else goto overflow; - if (1==z2) { + if (0L==iz1) return z1; + if (1L==iz2) { tmp = iacc*iz1; if (tmp/iacc != iz1) { overflow: @@ -2165,20 +2436,20 @@ SCM scm_intexpt(z1, z2) acc = long2num(tmp); break; } - if (z2 & 1) { + if (iz2 & 1) { tmp = iacc*iz1; if (tmp/iacc != iz1) goto overflow; iacc = tmp; - z2 = z2 - 1; /* so jumping to gencase works */ + iz2 = iz2 - 1L; /* so jumping to gencase works */ } tmp = iz1*iz1; if (tmp/iz1 != iz1) goto overflow; iz1 = tmp; - z2 >>= 1; + iz2 >>= 1; } #ifndef RECKLESS if (FALSEP(acc)) - errout: wta(UNDEFINED, (char *)OVFLOW, s_intexpt); + errout: wta(UNDEFINED, (char *)OVFLOW, s_intexpt); #endif goto ret; } @@ -2186,27 +2457,50 @@ SCM scm_intexpt(z1, z2) #ifdef FLOATS if (REALP(z1)) { dz1 = REALPART(z1); - flocase: dacc = 1.0; while(1) { - if (0==z2) break; - if (1==z2) {dacc = dacc*dz1; break;} - if (z2 & 1) dacc = dacc*dz1; + if (0L==iz2) break; + if (1L==iz2) {dacc = dacc*dz1; break;} + if (iz2 & 1) dacc = dacc*dz1; dz1 = dz1*dz1; - z2 >>= 1; + iz2 >>= 1; } - return makdbl(recip ? 1.0/dacc : dacc, 0.0); + return makdbl(dacc, 0.0); } #endif gencase: - while(1) { - 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; + { + SCM tz1 = z1; + while(1) { + SCM tmp; + if (0L==iz2) break; + if (1L==iz2) { + tmp = acc; + acc = product(tmp, tz1); +#ifdef BIGDIG + bigrecy(tmp); +#endif + break; + } + if (iz2 & 1) { + tmp = acc; + acc = product(tmp, tz1); +#ifdef BIGDIG + bigrecy(tmp); +#endif + } + tmp = tz1; + tz1 = product(tmp, tmp); +#ifdef BIGDIG + if (tmp != z1 && tmp != acc) bigrecy(tmp); +#endif + iz2 >>= 1; + } +#ifdef BIGDIG + if (tz1 != acc && tz1 != z1) bigrecy(tz1); +#endif } - ret: return recip ? divide(acc, UNDEFINED) : acc; + ret: return acc; } #ifdef FLOATS @@ -2258,7 +2552,7 @@ void two_doubles(z1, z2, sstring, xy) else { # ifdef BIGDIG ASRTGO(NIMP(z1), badz1); - if (BIGP(z1)) xy->x = big2dbl(z1); + if (BIGP(z1)) xy->x = int2dbl(z1); else { # ifndef RECKLESS if (!(REALP(z1))) @@ -2274,7 +2568,7 @@ void two_doubles(z1, z2, sstring, xy) else { # ifdef BIGDIG ASRTGO(NIMP(z2), badz2); - if (BIGP(z2)) xy->y = big2dbl(z2); + if (BIGP(z2)) xy->y = int2dbl(z2); else { # ifndef RECKLESS if (!(REALP(z2))) @@ -2425,7 +2719,7 @@ SCM ex2in(z) ASRTGO(NIMP(z), badz); if (INEXP(z)) return z; # ifdef BIGDIG - if (BIGP(z)) return makdbl(big2dbl(z), 0.0); + if (BIGP(z)) return makdbl(int2dbl(z), 0.0); # endif badz: wta(z, (char *)ARG1, s_ex2in); } @@ -2508,50 +2802,50 @@ SCM dbl2big(d) BIGDIG *digits; SCM ans; double u = (d < 0)?-d:d; + ASRTER(u == u && u != u/2, makdbl(d, 0.0), OVFLOW, "dbl2big"); while (0 != floor(u)) {u /= BIGRAD;i++;} ans = mkbig(i, d < 0); digits = BDIGITS(ans); while (i--) { - u *= BIGRAD; + u = ldexp(u, BITSPERDIG); c = floor(u); u -= c; digits[i] = c; } - ASRTER(0==u, INUM0, OVFLOW, "dbl2big"); - return ans; + ASRTER(0==u, makdbl(d, 0.0), OVFLOW, "dbl2big"); + return normbig(ans); } -double big2dbl(b) +/* This turns out to need explicit rounding for bignums */ +double int2dbl(b) SCM b; { - double ans = 0.0; - sizet i = NUMDIGS(b); - BIGDIG *digits = BDIGITS(b); - while (i--) ans = digits[i] + BIGRAD*ans; - 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; + if (INUMP(b)) return (double)(INUM(b)); + { + SCM num = scm_iabs(b), quo = num; + int bex = INUM(scm_intlength(num)) - dbl_mant_dig; + double ans = 0.0; + if (bex > 0) { + SCM den = scm_ash(MAKINUM(1L), MAKINUM(bex)); + quo = scm_round_quotient(num, den); + bigrecy(den); + } + /* quo may not be a bignum */ + if (INUMP(quo)) + ans = ldexp((double)(INUM(quo)), bex); + else { + sizet i = NUMDIGS(quo); + sizet j = i - (dbl_mant_dig + BITSPERDIG - 1)/BITSPERDIG; + BIGDIG *digits = BDIGITS(quo); + if (j < 0) j = 0; + while (i-- > j) ans = digits[i] + ldexp(ans, BITSPERDIG); + bex += j * BITSPERDIG; + if (bex > 0) ans = ldexp(ans, bex); + } + if (num != b) bigrecy(num); + if (quo != b) bigrecy(quo); + if (tc16_bigneg==TYP16(b)) return -ans; + return ans; + } } static SCM bigdblop(op, b, re, im) int op; @@ -2561,54 +2855,65 @@ static SCM bigdblop(op, b, re, im) double bm = 0.0; int i = 0; if (NUMDIGS(b)*BITSPERDIG < DBL_MAX_EXP) { - bm = big2dbl(b); + bm = int2dbl(b); } else { i = INUM(scm_intlength(b)); if (i < DBL_MAX_EXP) { i = 0; - bm = big2dbl(b); + bm = int2dbl(b); } else { i = i + 1 - DBL_MAX_EXP; - bm = big2scaldbl(b, i); + bm = ldexp(int2dbl(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 '/': + if (0.0==im) return makdbl(bm/re, 0.0); + { + 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)); default: return UNSPECIFIED; } } -static SCM inex_divbigbig(a, b) +/* now able to return unnomalized doubles. */ +static SCM inex_divintbig(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); + { + int sgn = (((INUMP(a) ? (INUM(a) < 0):BIGSIGN(a))==0) ^ + (BIGSIGN(b)==0)) ? -1 : 1; + SCM ma = scm_abs(a); + SCM mb = scm_abs(b); + int la = INUM(scm_intlength(ma)); + int lb = INUM(scm_intlength(mb)); + if (la <= DBL_MAX_EXP && lb <= DBL_MAX_EXP) { + r = int2dbl(a) / int2dbl(b); + } + else if (la > DBL_MAX_EXP && lb > DBL_MAX_EXP) { + int k = (la > lb ? la : lb) - DBL_MAX_EXP; + r = sgn * + int2dbl(scm_ash(ma, MAKINUM(-k))) / + int2dbl(scm_ash(mb, MAKINUM(-k))); + } else if (la > lb) { + int k = la - DBL_MAX_EXP; + r = sgn * ldexp(int2dbl(scm_ash(ma, MAKINUM(-k))) / int2dbl(mb), k); + } else { + int k = lb - DBL_MAX_EXP; + r = sgn * ldexp(int2dbl(ma) / int2dbl(scm_ash(mb, MAKINUM(-k))), -k); } } return makdbl(r, 0.0); } -static char s_dfloat_parts[] = "double-float-parts"; SCM scm_dfloat_parts(f) SCM f; { @@ -2629,42 +2934,10 @@ SCM scm_make_dfloat(mant, expt) double dmant = num2dbl(mant, (char *)ARG1, s_make_dfloat); int e = INUM(expt); ASRTER(INUMP(expt), expt, ARG2, s_make_dfloat); - ASRTER((dmant < 0 ? -dmant : dmant)<=max_dbl_int, mant, + ASRTER((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 @@ -2812,8 +3085,8 @@ static iproc subr1s[] = { static iproc asubrs[] = { {s_difference, difference}, {s_divide, divide}, - {s_max, lmax}, - {s_min, lmin}, + {s_max, scm_max}, + {s_min, scm_min}, {s_sum, sum}, {s_product, product}, {0, 0}}; @@ -2826,7 +3099,6 @@ static iproc subr2s[] = { {s_expt, expt}, # ifdef BIGDIG {s_make_dfloat, scm_make_dfloat}, - {s_next_dfloat, scm_next_dfloat}, # endif #endif #ifdef INUMS_ONLY @@ -2837,6 +3109,8 @@ static iproc subr2s[] = { {s_assv, assv}, #endif {s_intexpt, scm_intexpt}, + {s_intlog, scm_intlog}, + {s_cintlog, scm_cintlog}, {s_list_tail, list_tail}, {s_ve_fill, vector_fill}, {s_st_fill, string_fill}, @@ -2926,9 +3200,10 @@ void init_scl() IMAG(scm_narn) = 0.0/0.0; ALLOW_INTS; # endif -# ifdef DBL_DIG +# ifndef BIGDIG +# ifdef DBL_DIG dblprec = (DBL_DIG > 20) ? 20 : DBL_DIG; -# else +# else { /* determine floating point precision */ double f = 0.1; volatile double fsum = 1.0+f; @@ -2939,7 +3214,8 @@ void init_scl() } dblprec = dblprec-1; } -# endif /* DBL_DIG */ +# endif /* DBL_DIG */ +# endif /* !BIGDIG */ # ifdef DBL_MANT_DIG dbl_mant_dig = DBL_MANT_DIG; # else @@ -2957,7 +3233,9 @@ void init_scl() # 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); + /* dbl_eps = ldexp(1.0, - dbl_mant_dig); */ sysintern("double-float-mantissa-length", MAKINUM(dbl_mant_dig)); + sysintern("bigdbl:powers-of-5", + pows5 = make_vector(MAKINUM(326), MAKINUM(1))); #endif } -- cgit v1.2.3