summaryrefslogtreecommitdiffstats
path: root/scl.c
diff options
context:
space:
mode:
Diffstat (limited to 'scl.c')
-rw-r--r--scl.c341
1 files changed, 277 insertions, 64 deletions
diff --git a/scl.c b/scl.c
index 1b4507d..51d9eab 100644
--- a/scl.c
+++ b/scl.c
@@ -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
}