summaryrefslogtreecommitdiffstats
path: root/scl.c
diff options
context:
space:
mode:
Diffstat (limited to 'scl.c')
-rw-r--r--scl.c251
1 files changed, 212 insertions, 39 deletions
diff --git a/scl.c b/scl.c
index 51d9eab..57d020e 100644
--- a/scl.c
+++ b/scl.c
@@ -15,26 +15,26 @@
* the Free Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111, USA.
*
* As a special exception, the Free Software Foundation gives permission
- * for additional uses of the text contained in its release of GUILE.
+ * for additional uses of the text contained in its release of SCM.
*
- * The exception is that, if you link the GUILE library with other files
+ * The exception is that, if you link the SCM library with other files
* to produce an executable, this does not by itself cause the
* resulting executable to be covered by the GNU General Public License.
* Your use of that executable is in no way restricted on account of
- * linking the GUILE library code into it.
+ * linking the SCM library code into it.
*
* This exception does not however invalidate any other reasons why
* the executable file might be covered by the GNU General Public License.
*
* This exception applies only to the code released by the
- * Free Software Foundation under the name GUILE. If you copy
+ * Free Software Foundation under the name SCM. If you copy
* code from other Free Software Foundation releases into a copy of
- * GUILE, as the General Public License permits, the exception does
+ * SCM, as the General Public License permits, the exception does
* not apply to the code that you add in this way. To avoid misleading
* anyone as to the status of such modified files, you must delete
* this exception notice from them.
*
- * If you write modifications of your own for GUILE, it is your choice
+ * If you write modifications of your own for SCM, it is your choice
* whether to permit this exception to apply to your modifications.
* If you do not wish that, delete this exception notice.
*/
@@ -45,11 +45,19 @@
#include "scm.h"
#ifdef FLOATS
-# include <math.h>
+# ifndef PLAN9
+# 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 int apx_log10 P((double x));
+static double lpow10 P((double x, int n));
+static sizet idbl2str P((double f, char *a));
+static sizet iflo2str P((SCM flt, char *str));
+static void add1 P((double f, double *fsum));
+static long scm_twos_power P((SCM n));
static char s_makrect[] = "make-rectangular", s_makpolar[] = "make-polar",
s_magnitude[] = "magnitude", s_angle[] = "angle",
@@ -83,7 +91,13 @@ static char s_intexpt[] = "integer-expt";
/*** NUMBERS -> STRINGS ***/
#ifdef FLOATS
-static int dbl_mant_dig;
+# ifndef DBL_MANT_DIG
+# define DBL_MANT_DIG dbl_mant_dig
+# endif
+static int dbl_mant_dig = 0;
+static double max_dbl_int; /* Integers less than or equal to max_dbl_int
+ are representable exactly as doubles. */
+static double dbl_eps;
double dbl_prec(x)
double x;
{
@@ -91,10 +105,10 @@ double dbl_prec(x)
double frac = frexp(x, &expt);
# ifdef DBL_MIN_EXP
if (0.0==x || expt < DBL_MIN_EXP) /* gradual underflow */
- return ldexp(1.0, -dbl_mant_dig) * ldexp(1.0, DBL_MIN_EXP);
+ return ldexp(1.0, -DBL_MANT_DIG) * ldexp(1.0, DBL_MIN_EXP);
# endif
- if (1.0==frac) return ldexp(1.0, expt - dbl_mant_dig + 1);
- return ldexp(1.0, expt - dbl_mant_dig);
+ if (1.0==frac) return ldexp(1.0, expt - DBL_MANT_DIG + 1);
+ return ldexp(1.0, expt - DBL_MANT_DIG);
}
static double llog2 = 0.3010299956639812; /* log10(2) */
@@ -130,7 +144,7 @@ static double lpow10(x, n)
/* DBL2STR_FUZZ is a somewhat arbitrary guard against
round off error in scaling f and fprec. */
-#define DBL2STR_FUZZ 0.9
+# define DBL2STR_FUZZ 0.9
int dblprec;
static sizet idbl2str(f, a)
double f;
@@ -479,14 +493,23 @@ SCM istr2int(str, len, radix)
if (c >= radix) return BOOL_F; /* bad digit for radix */
ln = n;
n = n * radix - c;
- /* Negation is a workaround for HP700 cc bug */
- if (n > ln || (-n > -MOST_NEGATIVE_FIXNUM)) goto ovfl;
+ if (n > ln
+# ifdef hpux
+ || (-n > -MOST_NEGATIVE_FIXNUM) /* workaround for HP700 cc bug */
+# endif
+ ) goto ovfl;
break;
default:
return BOOL_F; /* not a digit */
}
} while (i < len);
- if (!lead_neg) if ((n = -n) > MOST_POSITIVE_FIXNUM) goto ovfl;
+ if (lead_neg) {
+ if (n < MOST_NEGATIVE_FIXNUM) goto ovfl;
+ }
+ else {
+ if (n < -MOST_POSITIVE_FIXNUM) goto ovfl;
+ n = -n;
+ }
return MAKINUM(n);
ovfl: /* overflow scheme integer */
return BOOL_F;
@@ -494,6 +517,34 @@ SCM istr2int(str, len, radix)
#endif
#ifdef FLOATS
+# ifdef BIGDIG
+static char twostab[] = {4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0};
+static long scm_twos_power(n)
+ SCM n;
+{
+ long d, c = 0;
+ int d4;
+# ifdef BIGDIG
+ if NINUMP(n) {
+ BIGDIG *ds;
+ int i = 0;
+ ds = BDIGITS(n);
+ while (0==(d = ds[i++])) c += BITSPERDIG;
+ goto out;
+ }
+# endif
+ d = INUM(n);
+ if (0==d) return 0;
+ out:
+ do {
+ d4 = 15 & d;
+ c += twostab[d4];
+ d >>= 4;
+ } while (0==d4);
+ return c;
+}
+# endif /* def BIGDIG */
+
SCM istr2flo(str, len, radix)
register char *str;
register long len;
@@ -965,6 +1016,29 @@ SCM numberp(x)
return BOOL_F;
}
#ifdef FLOATS
+# ifdef BIGDIG
+int scm_bigdblcomp(b, d)
+ SCM b;
+ double d;
+{
+ sizet dlen, blen;
+ int dneg = d < 0 ? 1 : 0;
+ int bneg = BIGSIGN(b) ? 1 : 0;
+ if (bneg < dneg) return -1;
+ if (bneg > dneg) return 1;
+ frexp(d, &dlen);
+ blen = INUM(scm_intlength(b));
+ if (blen > dlen) return dneg ? 1 : -1;
+ if (blen < dlen) return dneg ? -1 : 1;
+ if ((blen <= DBL_MANT_DIG) || (blen - scm_twos_power(b)) <= DBL_MANT_DIG) {
+ double bd = big2dbl(b);
+ if (bd > d) return -1;
+ if (bd < d) return 1;
+ return 0;
+ }
+ return bigcomp(b, dbl2big(d));
+}
+# endif
SCM realp(x)
SCM x;
{
@@ -1018,7 +1092,8 @@ SCM eqp(x, y)
if BIGP(y) return (0==bigcomp(x, y)) ? BOOL_T : BOOL_F;
ASRTGO(INEXP(y), bady);
bigreal:
- return (REALP(y) && (big2dbl(x)==REALPART(y))) ? BOOL_T : BOOL_F;
+ return (REALP(y) && (0==scm_bigdblcomp(x, REALPART(y)))) ?
+ BOOL_T : BOOL_F;
}
ASRTGO(INEXP(x), badx);
# else
@@ -1091,7 +1166,7 @@ SCM lessp(x, y)
ASRTGO(NIMP(y), bady);
if BIGP(y) return (1==bigcomp(x, y)) ? BOOL_T : BOOL_F;
ASRTGO(REALP(y), bady);
- return (big2dbl(x) < REALPART(y)) ? BOOL_T : BOOL_F;
+ return (1==scm_bigdblcomp(x, REALPART(y))) ? BOOL_T : BOOL_F;
}
ASRTGO(REALP(x), badx);
# else
@@ -1100,7 +1175,7 @@ SCM lessp(x, y)
if INUMP(y) return (REALPART(x) < ((double)INUM(y))) ? BOOL_T : BOOL_F;
# ifdef BIGDIG
ASRTGO(NIMP(y), bady);
- if BIGP(y) return (REALPART(x) < big2dbl(y)) ? BOOL_T : BOOL_F;
+ if BIGP(y) return (-1==scm_bigdblcomp(y, REALPART(x))) ? BOOL_T : BOOL_F;
ASRTGO(REALP(y), bady);
# else
ASRTGO(NIMP(y) && REALP(y), bady);
@@ -1248,10 +1323,12 @@ SCM negativep(x)
return (x < INUM0) ? BOOL_T : BOOL_F;
}
+static char s_exactprob[] = "not representable as inexact";
SCM lmax(x, y)
SCM x, y;
{
#ifdef FLOATS
+ SCM t;
double z;
#endif
if UNBNDP(y) {
@@ -1270,8 +1347,11 @@ SCM lmax(x, y)
ASRTGO(NIMP(y), bady);
if BIGP(y) return (1==bigcomp(x, y)) ? y : x;
ASRTGO(REALP(y), bady);
+ big_dbl:
+ if (-1 != scm_bigdblcomp(x, REALPART(y))) return y;
z = big2dbl(x);
- return (z < REALPART(y)) ? y : makdbl(z, 0.0);
+ ASSERT(0==scm_bigdblcomp(x, z), x, s_exactprob, s_max);
+ return makdbl(z, 0.0);
}
ASRTGO(REALP(x), badx);
# else
@@ -1280,7 +1360,9 @@ SCM lmax(x, y)
if INUMP(y) return (REALPART(x) < (z = INUM(y))) ? makdbl(z, 0.0) : x;
# ifdef BIGDIG
ASRTGO(NIMP(y), bady);
- if BIGP(y) return (REALPART(x) < (z = big2dbl(y))) ? makdbl(z, 0.0) : x;
+ if BIGP(y) {
+ t = y; y = x; x = t; goto big_dbl;
+ }
ASRTGO(REALP(y), bady);
# else
ASRTGO(NIMP(y) && REALP(y), bady);
@@ -1330,6 +1412,7 @@ SCM lmin(x, y)
SCM x, y;
{
#ifdef FLOATS
+ SCM t;
double z;
#endif
if UNBNDP(y) {
@@ -1348,8 +1431,11 @@ SCM lmin(x, y)
ASRTGO(NIMP(y), bady);
if BIGP(y) return (-1==bigcomp(x, y)) ? y : x;
ASRTGO(REALP(y), bady);
+ big_dbl:
+ if (1 != scm_bigdblcomp(x, REALPART(y))) return y;
z = big2dbl(x);
- return (z > REALPART(y)) ? y : makdbl(z, 0.0);
+ ASSERT(0==scm_bigdblcomp(x, z), x, s_exactprob, s_min);
+ return makdbl(z, 0.0);
}
ASRTGO(REALP(x), badx);
# else
@@ -1358,7 +1444,9 @@ SCM lmin(x, y)
if INUMP(y) return (REALPART(x) > (z = INUM(y))) ? makdbl(z, 0.0) : x;
# ifdef BIGDIG
ASRTGO(NIMP(y), bady);
- if BIGP(y) return (REALPART(x) > (z = big2dbl(y))) ? makdbl(z, 0.0) : x;
+ if BIGP(y) {
+ t = y; y = x; x = t; goto big_dbl;
+ }
ASRTGO(REALP(y), bady);
# else
ASRTGO(NIMP(y) && REALP(y), bady);
@@ -1960,6 +2048,10 @@ SCM scm_intexpt(z1, z2)
SCM z1, z2;
{
SCM acc = MAKINUM(1L);
+ int recip = 0;
+#ifdef FLOATS
+ double dacc, dz1;
+#endif
#ifdef BIGDIG
if (INUM0==z1 || acc==z1) return z1;
else if (MAKINUM(-1L)==z1) return BOOL_F==evenp(z2)?z1:acc;
@@ -1968,18 +2060,17 @@ SCM scm_intexpt(z1, z2)
z2 = INUM(z2);
if (z2 < 0) {
z2 = -z2;
- z1 = divide(z1, UNDEFINED);
+ recip = 1; /* z1 = divide(z1, UNDEFINED); */
}
if INUMP(z1) {
long tmp, iacc = 1, iz1 = INUM(z1);
+#ifdef FLOATS
+ if (recip) { dz1 = iz1; goto flocase; }
+#endif
while (1) {
if (0==z2) {
acc = long2num(iacc);
-#ifndef RECKLESS
- if (FALSEP(z1))
- errout: wta(UNDEFINED, (char *)OVFLOW, s_intexpt);
-#endif
- return acc;
+ break;
}
if (1==z2) {
tmp = iacc*iz1;
@@ -1991,8 +2082,7 @@ SCM scm_intexpt(z1, z2)
goto gencase;
}
acc = long2num(tmp);
- ASRTGO(NFALSEP(acc), errout);
- return acc;
+ break;
}
if (z2 & 1) {
tmp = iacc*iz1;
@@ -2005,28 +2095,37 @@ SCM scm_intexpt(z1, z2)
iz1 = tmp;
z2 >>= 1;
}
+#ifndef RECKLESS
+ if (FALSEP(acc))
+ errout: wta(UNDEFINED, (char *)OVFLOW, s_intexpt);
+#endif
+ goto ret;
}
ASSERT(NIMP(z1), z1, ARG1, s_intexpt);
-#ifdef FLOATS /* Move to scl.c ? */
+#ifdef FLOATS
if REALP(z1) {
- double dacc = 1.0, dz1 = REALPART(z1);
+ dz1 = REALPART(z1);
+ flocase:
+ dacc = 1.0;
while(1) {
- if (0==z2) return makdbl(dacc, 0.0);
- if (1==z2) return makdbl(dacc*dz1, 0.0);
+ if (0==z2) break;
+ if (1==z2) {dacc = dacc*dz1; break;}
if (z2 & 1) dacc = dacc*dz1;
dz1 = dz1*dz1;
z2 >>= 1;
}
+ return makdbl(recip ? 1.0/dacc : dacc, 0.0);
}
#endif
gencase:
while(1) {
- if (0==z2) return acc;
- if (1==z2) return product(acc, z1);
+ if (0==z2) break;
+ if (1==z2) {acc = product(acc, z1); break;}
if (z2 & 1) acc = product(acc, z1);
z1 = product(z1, z1);
z2 >>= 1;
}
+ ret: return recip ? divide(acc, UNDEFINED) : acc;
}
#ifdef FLOATS
@@ -2242,8 +2341,11 @@ SCM in2ex(z)
# ifdef BIGDIG
{
double u = floor(REALPART(z)+0.5);
- if ((u <= MOST_POSITIVE_FIXNUM) && (-u <= -MOST_NEGATIVE_FIXNUM)) {
- /* Negation is a workaround for HP700 cc bug */
+ if ((u <= MOST_POSITIVE_FIXNUM)
+# ifdef hpux
+ && (-u <= -MOST_NEGATIVE_FIXNUM) /* workaround for HP700 cc bug */
+# endif
+ ) {
SCM ans = MAKINUM((long)u);
if (INUM(ans)==(long)u) return ans;
}
@@ -2349,6 +2451,8 @@ static SCM bigdblop(op, b, re, im)
}
case '\\':
return makdbl(ldexp(re/bm, -i), 0.0==im ? 0.0 : ldexp(im/bm, -i));
+ default:
+ return UNSPECIFIED;
}
}
static SCM inex_divbigbig(a, b)
@@ -2371,6 +2475,64 @@ static SCM inex_divbigbig(a, b)
}
return makdbl(r, 0.0);
}
+
+static char s_dfloat_parts[] = "double-float-parts";
+SCM scm_dfloat_parts(f)
+ SCM f;
+{
+ int expt, ndig = DBL_MANT_DIG;
+ double mant = frexp(num2dbl(f, (char *)ARG1, s_dfloat_parts), &expt);
+# ifdef DBL_MIN_EXP
+ if (expt < DBL_MIN_EXP)
+ ndig -= DBL_MIN_EXP - expt;
+# endif
+ mant *= ldexp(1.0, ndig);
+ expt -= ndig;
+ return scm_values(dbl2big(mant), MAKINUM(expt), EOL, s_dfloat_parts);
+}
+static char s_make_dfloat[] = "make-double-float";
+SCM scm_make_dfloat(mant, expt)
+ SCM mant, expt;
+{
+ double dmant = num2dbl(mant, (char *)ARG1, s_make_dfloat);
+ int e = INUM(expt);
+ ASSERT(INUMP(expt), expt, ARG2, s_make_dfloat);
+ ASSERT((dmant < 0 ? -dmant : dmant)<=max_dbl_int, mant,
+ OUTOFRANGE, s_make_dfloat);
+ return makdbl(ldexp(dmant, e), 0.0);
+}
+static char s_next_dfloat[] = "next-double-float";
+SCM scm_next_dfloat(f1, f2)
+ SCM f1, f2;
+{
+ int e, neg = 0;
+ double d1 = num2dbl(f1, (char *)ARG1, s_next_dfloat);
+ double dif = num2dbl(f2, (char *)ARG2, s_next_dfloat) - d1;
+ double d = frexp(d1, &e), eps = dbl_eps;
+ if (d1 < 0) {neg = 1; dif = -dif; d = -d;}
+ if (dif > 0) {
+# ifdef DBL_MIN_EXP
+ if (e < DBL_MIN_EXP)
+ eps = ldexp(eps, DBL_MIN_EXP - e);
+ else if (0.0==d)
+ eps = ldexp(1.0, DBL_MIN_EXP - DBL_MANT_DIG);
+# endif
+ d = ldexp(d + eps, e);
+ }
+ else if (dif < 0) {
+# ifdef DBL_MIN_EXP
+ if (e < DBL_MIN_EXP)
+ eps = ldexp(eps, DBL_MIN_EXP - e);
+ else if (0.0==d)
+ eps = ldexp(-1.0, DBL_MIN_EXP - DBL_MANT_DIG);
+# endif
+ if (0.5==d) eps *= 0.5;
+ d = ldexp(d - eps, e);
+ }
+ else if (0.0==dif)
+ return f1;
+ return makdbl(neg ? -d : d, 0.0);
+}
# endif
#endif
@@ -2490,6 +2652,9 @@ static iproc subr1s[] = {
{s_angle, angle},
{s_in2ex, in2ex},
{s_ex2in, ex2in},
+# ifdef BIGDIG
+ {s_dfloat_parts, scm_dfloat_parts},
+# endif
#else
{"real?", numberp},
{"rational?", numberp},
@@ -2524,6 +2689,10 @@ static iproc subr2s[] = {
{s_makpolar, makpolar},
{s_atan2, latan2},
{s_expt, expt},
+# ifdef BIGDIG
+ {s_make_dfloat, scm_make_dfloat},
+ {s_next_dfloat, scm_next_dfloat},
+# endif
#endif
#ifdef INUMS_ONLY
{s_memv, memq},
@@ -2634,7 +2803,7 @@ void init_scl()
# ifdef DBL_MANT_DIG
dbl_mant_dig = DBL_MANT_DIG;
# else
- {
+ if (!DBL_MANT_DIG) { /* means we #defined it. */
double fsum = 0.0, eps = 1.0;
int i = 0;
while (fsum != 1.0) {
@@ -2645,5 +2814,9 @@ void init_scl()
dbl_mant_dig = i;
}
# endif /* DBL_MANT_DIG */
+ max_dbl_int = pow(2.0, dbl_mant_dig - 1.0);
+ max_dbl_int = max_dbl_int + (max_dbl_int - 1.0);
+ dbl_eps = ldexp(1.0, -dbl_mant_dig);
+ sysintern("double-float-mantissa-length", MAKINUM(DBL_MANT_DIG));
#endif
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-04-28 07:37:07 +0000