aboutsummaryrefslogtreecommitdiffstats
path: root/scl.c
diff options
context:
space:
mode:
Diffstat (limited to 'scl.c')
-rw-r--r--scl.c898
1 files changed, 588 insertions, 310 deletions
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 <math.h>
# 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
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-04-16 04:36:46 +0000