aboutsummaryrefslogtreecommitdiffstats
path: root/subr.c
diff options
context:
space:
mode:
Diffstat (limited to 'subr.c')
-rwxr-xr-x[-rw-r--r--]subr.c471
1 files changed, 332 insertions, 139 deletions
diff --git a/subr.c b/subr.c
index 95eed21..57e035a 100644..100755
--- a/subr.c
+++ b/subr.c
@@ -1,5 +1,5 @@
/* "subr.c" integer and other Scheme procedures
- * Copyright (C) 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997 Free Software Foundation, Inc.
+ * Copyright (C) 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 2013 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
@@ -36,9 +36,10 @@ static char s_symbol2string[] = "symbol->string",
extern char s_inexactp[];
#define s_exactp (s_inexactp+2)
static char s_oddp[] = "odd?", s_evenp[] = "even?";
-static char s_quotient[] = "quotient",
+static char s_rquotient[] = "round-quotient",
s_remainder[] = "remainder", s_modulo[] = "modulo";
static char s_gcd[] = "gcd";
+#define s_quotient (s_rquotient+6)
static char s_ci_eq[] = "char-ci=?",
s_ch_lessp[] = "char<?", s_ch_leqp[] = "char<=?",
@@ -306,6 +307,104 @@ SCM evenp(n)
#endif
return (4 & (int)n) ? BOOL_F : BOOL_T;
}
+
+SCM scm_round_quotient(num, den)
+ SCM num, den;
+{
+ register long quo, rem;
+ /* if (scm_verbose > 1) */
+ /* printf("%s / %s\n", */
+ /* CHARS(number2string(num, MAKINUM(10))), */
+ /* CHARS(number2string(den, MAKINUM(10)))); */
+#ifdef BIGDIG
+ if (NINUMP(num)) {
+ long w;
+ ASRTER(NIMP(num) && BIGP(num), num, ARG1, s_rquotient);
+ if (NINUMP(den)) {
+ ASRTGO(NIMP(den) && BIGP(den), badden);
+ return divbigbig(BDIGITS(num), NUMDIGS(num), BDIGITS(den), NUMDIGS(den),
+ BIGSIGN(num) ^ BIGSIGN(den), 3);
+ }
+ if (!(quo = INUM(den))) goto ov;
+ if (1==quo) return num;
+ /* divbigdig() hasn't been extended to perform rounding */
+ /* if (quo < 0) quo = -quo; */
+ /* if (quo < BIGRAD) { */
+ /* w = copybig(num, BIGSIGN(num) ? (den>0) : (den<0)); */
+ /* divbigdig(BDIGITS(w), NUMDIGS(w), (BIGDIG)quo); */
+ /* return normbig(w); */
+ /* } */
+# ifndef DIGSTOOBIG
+ w = pseudolong(quo);
+ return divbigbig(BDIGITS(num), NUMDIGS(num), (BIGDIG *)&w, DIGSPERLONG,
+ BIGSIGN(num) ? (den>0) : (den<0), 3);
+# else
+ { BIGDIG quodigs[DIGSPERLONG];
+ longdigs(quo, quodigs);
+ return divbigbig(BDIGITS(num), NUMDIGS(num), quodigs, DIGSPERLONG,
+ BIGSIGN(num) ? (den>0) : (den<0), 3);
+ }
+# endif
+ }
+ if (NINUMP(den)) {
+# ifndef RECKLESS
+ if (!(NIMP(den) && BIGP(den)))
+ badden: wta(den, (char *)ARG2, s_rquotient);
+# endif
+ if (NUMDIGS(den) > DIGSPERLONG ||
+ (NUMDIGS(den)==DIGSPERLONG && BDIGITS(den)[DIGSPERLONG-1] >= BIGRAD/2))
+ return INUM0;
+ quo = num2long(den, (char *)ARG2, s_rquotient);
+ rem = INUM(num)%quo;
+ if (labs(2*rem) > labs(quo))
+ return MAKINUM(((INUM(num) < 0)==(quo < 0)) ? 1 : -1);
+ else return INUM0;
+ }
+#else
+ ASRTER(INUMP(num), num, ARG1, s_rquotient);
+ ASRTER(INUMP(den), den, ARG2, s_rquotient);
+#endif
+ if ((quo = INUM(den))==0)
+ ov: wta(den, (char *)OVFLOW, s_rquotient);
+ quo = INUM(num)/quo;
+ {
+# if (__TURBOC__==1)
+ rem = ((den<0) ? -INUM(num) : INUM(num))%INUM(den);
+# else
+ rem = INUM(num)%INUM(den);
+# endif
+#ifdef BADIVSGNS
+ if (rem==0) ;
+ else if (rem < 0) {
+ if (num < 0) ;
+ else quo--;
+ } else if (num < 0) quo++;
+#endif
+ if ((1 & quo)
+ ? labs(2*rem) >= labs(INUM(den))
+ : labs(2*rem) > labs(INUM(den)))
+ quo = quo + (((INUM(num) < 0)==(INUM(den) < 0)) ? 1 : -1);
+ }
+ if (!FIXABLE(quo))
+#ifdef BIGDIG
+ return long2big(quo);
+#else
+ wta(num, (char *)OVFLOW, s_rquotient);
+#endif
+ return MAKINUM(quo);
+}
+
+/* SCM scm_round_quotient(num, den) */
+/* SCM num, den; */
+/* { */
+/* SCM quo = lquotient(num, den); */
+/* SCM rem = lremainder(num, den); */
+/* if (BOOL_T==((BOOL_T==evenp(quo) ? greaterp : greqp) */
+/* (scm_ash(scm_iabs(rem), MAKINUM(1L)), scm_iabs(den)))) */
+/* quo = sum(quo, MAKINUM(negativep(num)==negativep(den) ? 1L : -1L)); */
+/* return quo; */
+/* } */
+
SCM lquotient(x, y)
SCM x, y;
{
@@ -319,8 +418,7 @@ SCM lquotient(x, y)
return divbigbig(BDIGITS(x), NUMDIGS(x), BDIGITS(y), NUMDIGS(y),
BIGSIGN(x) ^ BIGSIGN(y), 2);
}
- z = INUM(y);
- ASRTGO(z, ov);
+ if (!(z = INUM(y))) goto ov;
if (1==z) return x;
if (z < 0) z = -z;
if (z < BIGRAD) {
@@ -458,9 +556,9 @@ SCM lgcd(x, y)
SCM x, y;
{
register long u, v, k, t;
- tailrec:
if (UNBNDP(y)) return UNBNDP(x) ? INUM0 : x;
#ifdef BIGDIG
+ tailrec:
if (NINUMP(x)) {
big_gcd:
ASRTER(NIMP(x) && BIGP(x), x, ARG1, s_gcd);
@@ -567,7 +665,7 @@ SCM scm_big_ior P((BIGDIG *x, sizet nx, int xsgn, SCM bigy));
SCM scm_big_and P((BIGDIG *x, sizet nx, int xsgn, SCM bigy, int zsgn));
SCM scm_big_xor P((BIGDIG *x, sizet nx, int xsgn, SCM bigy));
SCM scm_big_test P((BIGDIG *x, sizet nx, int xsgn, SCM bigy));
-SCM scm_big_ash P((SCM x, long cnt));
+SCM scm_big_ash P((SCM x, int cnt));
SCM scm_copy_big_dec(b, sign)
SCM b;
@@ -770,8 +868,6 @@ static SCM scm_copy_big_2scomp(x, blen, sign)
sizet i;
if (INUMP(x)) {
long lx = INUM(x);
- if (nres < (LONG_BIT + BITSPERDIG - 1)/BITSPERDIG)
- nres = (LONG_BIT + BITSPERDIG - 1)/BITSPERDIG;
res = mkbig(nres, sign);
rds = BDIGITS(res);
if (lx < 0) {
@@ -848,15 +944,13 @@ static void scm_2scomp1(b)
SCM scm_big_ash(x, cnt)
SCM x;
- long cnt;
+ int cnt;
{
SCM res;
- BIGDIG *resds;
- unsigned long d;
- int sign, ishf;
- long i, fshf, blen, n;
+ BIGDIG *resds, d;
+ int sign, i, ishf, fshf, blen, n;
if (INUMP(x)) {
- blen = LONG_BIT;
+ blen = INUM(scm_intlength(x));
sign = INUM(x) < 0 ? 0x0100 : 0;
}
else {
@@ -871,12 +965,13 @@ SCM scm_big_ash(x, cnt)
resds = BDIGITS(res);
n = NUMDIGS(res) - ishf - 1;
for (i = 0; i < n; i++) {
- d = (resds[i + ishf]>>fshf) |
- ((resds[i + ishf + 1])<<(BITSPERDIG - fshf) & (BIGRAD - 1));
+ d = (resds[i + ishf]>>fshf);
+ if (fshf)
+ d |= ((resds[i + ishf + 1])<<(BITSPERDIG - fshf) & (BIGRAD - 1));
resds[i] = d;
}
d = (resds[i + ishf]>>fshf);
- if (sign) d |= ((BIGRAD - 1)<<(BITSPERDIG - fshf) & (BIGRAD - 1));
+ if (sign && fshf) d |= ((BIGRAD - 1)<<(BITSPERDIG - fshf) & (BIGRAD - 1));
resds[i] = d;
n = NUMDIGS(res);
d = sign ? BIGRAD - 1 : 0;
@@ -888,16 +983,18 @@ SCM scm_big_ash(x, cnt)
fshf = cnt % BITSPERDIG;
res = scm_copy_big_2scomp(x, blen + cnt, sign);
resds = BDIGITS(res);
- for (i = NUMDIGS(res) - 1; i > ishf; i--) {
- d = (((resds[i - ishf])<<fshf) & (BIGRAD - 1)) |
- ((resds[i - ishf - 1])>>(BITSPERDIG - fshf));
- resds[i] = d;
- }
- d = (((resds[i - ishf])<<fshf) & (BIGRAD - 1));
+ /* if (scm_verbose>1){for (i=NUMDIGS(res); i--;) printf(" %08x",resds[i]); printf("\n");} */
+ for (i = NUMDIGS(res) - 1; i > ishf; i--)
+ if (fshf) {
+ d = (((resds[i - ishf])<<fshf) & (BIGRAD - 1))
+ | ((resds[i - ishf - 1])>>(BITSPERDIG - fshf));
+ resds[i] = d;
+ } else resds[i] = resds[i - ishf];
+ d = fshf ? (((resds[i - ishf])<<fshf) & (BIGRAD - 1)) : resds[i - ishf];
resds[i] = d;
- for (i--; i >= 0; i--)
- resds[i] = 0;
+ for (i--; i >= 0; i--) resds[i] = 0;
}
+ /* if (scm_verbose>1){for (i=NUMDIGS(res); i--;) printf(" %08x",resds[i]); printf("\n");} */
if (sign) scm_2scomp1(res);
return normbig(res);
}
@@ -1170,17 +1267,18 @@ SCM scm_lognot(n)
SCM scm_ash(n, cnt)
SCM n, cnt;
{
- SCM res = INUM(n);
+ SCM res;
+ long ni = INUM(n);
+ int icnt = INUM(cnt);
ASRTER(INUMP(cnt), cnt, ARG2, s_ash);
- cnt = INUM(cnt);
if (INUMP(n)) {
- if (cnt < 0) {
- if (-cnt >= LONG_BIT) return INUM0;
- return MAKINUM(SRS(res, -cnt));
+ if (icnt < 0) {
+ if (-icnt >= LONG_BIT) return INUM0;
+ return MAKINUM(SRS(ni, -icnt));
}
- if (cnt >= LONG_BIT) goto ovflow;
- res = MAKINUM(res<<cnt);
- if (INUM(res)>>cnt != INUM(n))
+ if (icnt >= LONG_BIT) goto ovflow;
+ res = MAKINUM(ni<<icnt);
+ if (INUM(res)>>icnt != INUM(n))
goto ovflow;
else
return res;
@@ -1188,8 +1286,8 @@ SCM scm_ash(n, cnt)
#ifdef BIGDIG
ASRTER(NIMP(n) && BIGP(n), n, ARG1, s_ash);
ovflow:
- if (0==cnt) return n;
- return scm_big_ash(n, cnt);
+ if (0==icnt) return n;
+ return scm_big_ash(n, icnt);
#else
ovflow:
wta(n, INUMP(n) ? (char *)OVFLOW : (char *)ARG1, s_ash);
@@ -1201,10 +1299,11 @@ SCM scm_bitfield(n, start, end)
SCM n, start, end;
{
int sign;
+ int istart = INUM(start);
+ int iend = INUM(end);
ASRTER(INUMP(start), start, ARG2, s_bitfield);
ASRTER(INUMP(end), end, ARG3, s_bitfield);
- start = INUM(start); end = INUM(end);
- ASRTER(end >= start, MAKINUM(end), OUTOFRANGE, s_bitfield);
+ ASRTER(iend >= istart, MAKINUM(iend), OUTOFRANGE, s_bitfield);
#ifdef BIGDIG
if (NINUMP(n)) {
BIGDIG *ds;
@@ -1212,29 +1311,29 @@ SCM scm_bitfield(n, start, end)
ASRTER(NIMP(n) && BIGP(n), n, ARG1, s_bitfield);
sign = BIGSIGN(n);
big:
- if (sign) n = scm_copy_big_2scomp(n, (sizet)end, 0);
- n = scm_big_ash(n, -start);
+ if (sign) n = scm_copy_big_2scomp(n, (sizet)iend, 0);
+ n = scm_big_ash(n, -istart);
if (INUMP(n)) {
- if (end - start >= LONG_BIT - 2) return n;
- return MAKINUM(INUM(n) & ((1L<<(end - start)) - 1));
+ if (iend - istart >= LONG_BIT - 2) return n;
+ return MAKINUM(INUM(n) & ((1L<<(iend - istart)) - 1));
}
nd = NUMDIGS(n);
ds = BDIGITS(n);
- i = (end - start) / BITSPERDIG;
+ i = (iend - istart) / BITSPERDIG;
if (i >= nd) return n;
- ds[i] &= ((1 << ((end - start) % BITSPERDIG)) - 1);
+ ds[i] &= ((1 << ((iend - istart) % BITSPERDIG)) - 1);
for (++i; i < nd; i++) ds[i] = 0;
return normbig(n);
}
- if (end >= LONG_BIT - 2) {
+ if (iend >= LONG_BIT - 2) {
sign = INUM(n) < 0;
goto big;
}
#else
ASRTER(INUMP(n), n, ARG1, s_bitfield);
- ASRTER(end < LONG_BIT - 2, MAKINUM(end), OUTOFRANGE, s_bitfield);
+ ASRTER(iend < LONG_BIT - 2, MAKINUM(iend), OUTOFRANGE, s_bitfield);
#endif
- return MAKINUM((INUM(n)>>start) & ((1L<<(end - start)) - 1));
+ return MAKINUM((INUM(n)>>istart) & ((1L<<(iend - istart)) - 1));
}
SCM scm_bitif(mask, n0, n1)
@@ -1298,8 +1397,12 @@ SCM scm_bitwise_bit_count(n)
if (NINUMP(n)) {
sizet i; BIGDIG *ds, d;
ASRTER(NIMP(n) && BIGP(n), n, ARG1, s_bitwise_bit_count);
- if (BIGSIGN(n))
- return scm_lognot(scm_bitwise_bit_count(difference(MAKINUM(-1L), n)));
+ if (BIGSIGN(n)) {
+ SCM df = difference(MAKINUM(-1L), n);
+ SCM bc = scm_bitwise_bit_count(df);
+ bigrecy(df);
+ return scm_lognot(bc);
+ }
ds = BDIGITS(n);
for (i = NUMDIGS(n); i--; )
for (d = ds[i]; d; d >>= 4) c += logtab[15 & d];
@@ -1325,7 +1428,12 @@ SCM scm_logcount(n)
#ifdef BIGDIG
if (NINUMP(n)) {
ASRTER(NIMP(n) && BIGP(n), n, ARG1, s_logcount);
- if (BIGSIGN(n)) return scm_bitwise_bit_count(difference(MAKINUM(-1L), n));
+ if (BIGSIGN(n)) {
+ SCM df = difference(MAKINUM(-1L), n);
+ SCM bc = scm_bitwise_bit_count(df);
+ bigrecy(df);
+ return bc;
+ }
return scm_bitwise_bit_count(n);
}
#else
@@ -1347,7 +1455,12 @@ SCM scm_intlength(n)
if (NINUMP(n)) {
BIGDIG *ds, d;
ASRTER(NIMP(n) && BIGP(n), n, ARG1, s_intlength);
- if (BIGSIGN(n)) return scm_intlength(difference(MAKINUM(-1L), n));
+ if (BIGSIGN(n)) {
+ SCM df = difference(MAKINUM(-1L), n);
+ SCM si = scm_intlength(df);
+ bigrecy(df);
+ return si;
+ }
ds = BDIGITS(n);
d = ds[c = NUMDIGS(n)-1];
for (c *= BITSPERDIG; d; d >>= 4) {c += 4; l = ilentab[15 & d];}
@@ -1761,6 +1874,7 @@ SCM mkbig(nlen, sign)
ALLOW_INTS;
return v;
}
+/* big2inum() frees bignum b when it returns an INUM */
SCM big2inum(b, l)
SCM b;
sizet l;
@@ -1769,9 +1883,14 @@ SCM big2inum(b, l)
BIGDIG *tmp = BDIGITS(b);
while (l--) num = BIGUP(num) + tmp[l];
if (TYP16(b)==tc16_bigpos) {
- if (POSFIXABLE(num)) return MAKINUM(num);
+ if (POSFIXABLE(num)) {
+ bigrecy(b);
+ return MAKINUM(num);
+ }}
+ else if (UNEGFIXABLE(num)) {
+ bigrecy(b);
+ return MAKINUM(-(long)num);
}
- else if (UNEGFIXABLE(num)) return MAKINUM(-(long)num);
return b;
}
char s_adjbig[] = "adjbig";
@@ -1824,7 +1943,7 @@ SCM long2big(n)
if (n < 0) n = -n;
while (i < DIGSPERLONG) {
digits[i++] = BIGLO(n);
- n = BIGDN((unsigned long)n);
+ n = BIGDN(n);
}
return ans;
}
@@ -1853,7 +1972,7 @@ int bigcomp(x, y)
if (ysign > xsign) return -1;
if ((ylen = NUMDIGS(y)) > (xlen = NUMDIGS(x))) return (xsign) ? -1 : 1;
if (ylen < xlen) return (xsign) ? 1 : -1;
- while(xlen-- && (BDIGITS(y)[xlen]==BDIGITS(x)[xlen]));
+ while (xlen-- && (BDIGITS(y)[xlen]==BDIGITS(x)[xlen]));
if (-1==xlen) return 0;
return (BDIGITS(y)[xlen] > BDIGITS(x)[xlen]) ?
(xsign ? -1 : 1) : (xsign ? 1 : -1);
@@ -1890,7 +2009,7 @@ SCM addbig(x, nx, xsgn, bigy, sgny)
sizet nx; /* Assumes nx <= NUMDIGS(bigy) */
int xsgn, sgny; /* Assumes xsgn and sgny equal either 0 or 0x0100 */
{
- long num = 0;
+ SBIGLONG num = 0;
sizet i = 0, ny = NUMDIGS(bigy);
SCM z = copybig(bigy, BIGSIGN(bigy) ^ sgny);
BIGDIG *zds = BDIGITS(z);
@@ -1938,7 +2057,7 @@ SCM mulbig(x, nx, y, ny, sgn)
int sgn;
{
sizet i = 0, j = nx + ny;
- unsigned long n = 0;
+ UBIGLONG n = 0;
SCM z = mkbig(j, sgn);
BIGDIG *zds = BDIGITS(z);
while (j--) zds[j] = 0;
@@ -1946,7 +2065,7 @@ SCM mulbig(x, nx, y, ny, sgn)
j = 0;
if (x[i]) {
do {
- n += zds[i + j] + ((unsigned long) x[i] * y[j]);
+ n += zds[i + j] + ((UBIGLONG) x[i] * y[j]);
zds[i + j++] = BIGLO(n);
n = BIGDN(n);
} while (j < ny);
@@ -1955,12 +2074,12 @@ SCM mulbig(x, nx, y, ny, sgn)
} while (++i < nx);
return normbig(z);
}
-unsigned int divbigdig(ds, h, div)
+UBIGLONG divbigdig(ds, h, div)
BIGDIG *ds;
sizet h;
BIGDIG div;
{
- register unsigned long t2 = 0;
+ register UBIGLONG t2 = 0L;
while(h--) {
t2 = BIGUP(t2) + ds[h];
ds[h] = t2 / div;
@@ -1975,7 +2094,7 @@ SCM divbigint(x, z, sgn, mode)
{
if (z < 0) z = -z;
if (z < BIGRAD) {
- register unsigned long t2 = 0;
+ register UBIGLONG t2 = 0;
register BIGDIG *ds = BDIGITS(x);
sizet nd = NUMDIGS(x);
while(nd--) t2 = (BIGUP(t2) + ds[nd]) % z;
@@ -1984,7 +2103,7 @@ SCM divbigint(x, z, sgn, mode)
}
{
# ifndef DIGSTOOBIG
- unsigned long t2 = pseudolong(z);
+ UBIGLONG t2 = pseudolong(z);
return divbigbig(BDIGITS(x), NUMDIGS(x), (BIGDIG *)&t2,
DIGSPERLONG, sgn, mode);
# else
@@ -1994,108 +2113,180 @@ SCM divbigint(x, z, sgn, mode)
# endif
}
}
-SCM divbigbig(x, nx, y, ny, sgn, modes)
+
+static SCM scm_copy_big_ash1 P((BIGDIG *xds, sizet xlen, BIGDIG dscl));
+/* Make a copy of 2*xds and divide by dscl if dscl > 0 */
+SCM scm_copy_big_ash1 (xds, xlen, dscl)
+ BIGDIG *xds;
+ sizet xlen;
+ BIGDIG dscl;
+{
+ sizet rlen = xlen + 1, i;
+ SCM dencell = mkbig(rlen, 0);
+ BIGDIG *dends = BDIGITS(dencell);
+ dends[xlen] = xds[xlen-1]>>(BITSPERDIG - 1);
+ for (i = xlen - 1; i > 0; i--)
+ dends[i] = (((xds[i])<<1) & (BIGRAD - 1))
+ | ((xds[i-1])>>(BITSPERDIG - 1));
+ dends[0] = (((xds[0])<<1) & (BIGRAD - 1));
+ while(rlen && !dends[rlen-1]) rlen--;
+ if (dscl) {
+ divbigdig(dends, rlen, dscl);
+ while(rlen && !dends[rlen-1]) rlen--;
+ }
+ SETNUMDIGS(dencell, rlen, TYP16(dencell));
+ return dencell;
+}
+
+SCM divbigbig(x, xlen, y, ylen, sgn, mode)
BIGDIG *x, *y;
- sizet nx, ny;
- int sgn, modes;
- /* modes description
+ sizet xlen, ylen;
+ int sgn, mode;
+ /* mode description
0 remainder
1 modulo
2 quotient
- 3 quotient but returns 0 if division is not exact. */
-{
- sizet i = 0, j = 0;
- long num = 0;
- unsigned long t2 = 0;
- SCM z, newy;
- BIGDIG d = 0, qhat, *zds, *yds;
- /* algorithm requires nx >= ny */
- if (nx < ny)
- switch (modes) {
+ 3 quotient with round-toward-even
+ 4 quotient but returns NULL if division is not exact. */
+{
+ int roundup = 0; /* used for round-quotient */
+ sizet i = 0, j = 0; /* loop indexes */
+ SBIGLONG dds = 0; /* double-digit signed */
+ UBIGLONG ddu = 0; /* double-digit unsigned */
+ SCM quocell, dencell;
+ sizet rlen;
+ BIGDIG *quods, /* quotient digits */
+ *dends, /* scaled denominator digits */
+ dscl = 0, /* unscale quotient from scaled divisor */
+ qhat;
+ while(!y[ylen-1]) ylen--; /* in case y came in as a psuedolong */
+ if (xlen < ylen)
+ switch (mode) {
case 0: /* remainder -- just return x */
- z = mkbig(nx, sgn); zds = BDIGITS(z);
- do {zds[i] = x[i];} while (++i < nx);
- return z;
+ quocell = mkbig(xlen, sgn); quods = BDIGITS(quocell);
+ do {quods[i] = x[i];} while (++i < xlen);
+ return quocell;
case 1: /* modulo -- return y-x */
- z = mkbig(ny, sgn); zds = BDIGITS(z);
+ quocell = mkbig(ylen, sgn); quods = BDIGITS(quocell);
do {
- num += (long) y[i] - x[i];
- if (num < 0) {zds[i] = num + BIGRAD; num = -1;}
- else {zds[i] = num; num = 0;}
- } while (++i < nx);
- while (i < ny) {
- num += y[i];
- if (num < 0) {zds[i++] = num + BIGRAD; num = -1;}
- else {zds[i++] = num; num = 0;}
+ dds += (long) y[i] - x[i];
+ if (dds < 0) {quods[i] = dds + BIGRAD; dds = -1;}
+ else {quods[i] = dds; dds = 0;}
+ } while (++i < xlen);
+ while (i < ylen) {
+ dds += y[i];
+ if (dds < 0) {quods[i++] = dds + BIGRAD; dds = -1;}
+ else {quods[i++] = dds; dds = 0;}
}
goto doadj;
case 2: return INUM0; /* quotient is zero */
- case 3: return 0; /* the division is not exact */
+ case 3: /* round-toward-even */
+ /* Use dencell and dends variables to double the numerator */
+ dencell = scm_copy_big_ash1(x, xlen, dscl);
+ dends = BDIGITS(dencell);
+ rlen = NUMDIGS(dencell);
+ if (rlen < ylen) return INUM0;;
+ if (rlen > ylen) goto retone;
+ i = rlen;
+ while (i-- && (y[i]==dends[i]));
+ if (-1==i || (y[i] > dends[i])) return INUM0;
+ retone:
+ return MAKINUM(sgn ? -1 : 1);
+ case 4: return 0; /* the division is not exact */
}
-
- z = mkbig(nx==ny ? nx+2 : nx+1, sgn); zds = BDIGITS(z);
- if (nx==ny) zds[nx+1] = 0;
- while(!y[ny-1]) ny--; /* in case y came in as a psuedolong */
- if (y[ny-1] < (BIGRAD>>1)) { /* normalize operands */
- d = BIGRAD/(y[ny-1]+1);
- newy = mkbig(ny, 0); yds = BDIGITS(newy);
- while(j < ny)
- {t2 += (unsigned long) y[j]*d; yds[j++] = BIGLO(t2); t2 = BIGDN(t2);}
- y = yds; j = 0; t2 = 0;
- while(j < nx)
- {t2 += (unsigned long) x[j]*d; zds[j++] = BIGLO(t2); t2 = BIGDN(t2);}
- zds[j] = t2;
+ /* main algorithm requires xlen >= ylen */
+ quocell = mkbig(xlen==ylen ? xlen+2 : xlen+1, sgn); quods = BDIGITS(quocell);
+ if (xlen==ylen) quods[xlen+1] = 0;
+ if (y[ylen-1] < (BIGRAD>>1)) { /* normalize operands */
+ dscl = BIGRAD/(y[ylen-1]+1);
+ dencell = mkbig(ylen, 0); dends = BDIGITS(dencell);
+ while(j < ylen) {
+ ddu += (UBIGLONG) y[j]*dscl;
+ dends[j++] = BIGLO(ddu); ddu = BIGDN(ddu);
+ }
+ j = 0; ddu = 0; /* y = dends; */
+ while(j < xlen) {
+ ddu += (UBIGLONG) x[j]*dscl;
+ quods[j++] = BIGLO(ddu); ddu = BIGDN(ddu);
+ }
+ quods[j] = ddu;
+ } else {
+ dends = y;
+ quods[j = xlen] = 0;
+ while (j--) quods[j] = x[j];
}
- else {zds[j = nx] = 0; while (j--) zds[j] = x[j];}
- j = nx==ny ? nx+1 : nx; /* dividend needs more digits than divisor */
+ j = xlen==ylen ? xlen+1 : xlen; /* dividend needs more digits than divisor */
do { /* loop over digits of quotient */
- if (zds[j]==y[ny-1]) qhat = BIGRAD-1;
- else qhat = (BIGUP(zds[j]) + zds[j-1])/y[ny-1];
+ if (quods[j]==dends[ylen-1]) qhat = BIGRAD-1;
+ else qhat = (BIGUP(quods[j]) + quods[j-1])/dends[ylen-1];
if (!qhat) continue;
- i = 0; num = 0; t2 = 0;
+ i = 0; dds = 0; ddu = 0;
do { /* multiply and subtract */
- t2 += (unsigned long) y[i] * qhat;
- num += zds[j - ny + i] - BIGLO(t2);
- if (num < 0) {zds[j - ny + i] = num + BIGRAD; num = -1;}
- else {zds[j - ny + i] = num; num = 0;}
- t2 = BIGDN(t2);
- } while (++i < ny);
- num += zds[j - ny + i] - t2; /* borrow from high digit; don't update */
- while (num) { /* "add back" required */
- i = 0; num = 0; qhat--;
+ ddu += (UBIGLONG) dends[i] * qhat;
+ dds += quods[j - ylen + i] - BIGLO(ddu);
+ if (dds < 0) {quods[j - ylen + i] = dds + BIGRAD; dds = -1;}
+ else {quods[j - ylen + i] = dds; dds = 0;}
+ ddu = BIGDN(ddu);
+ } while (++i < ylen);
+ dds += quods[j - ylen + i] - ddu; /* borrow from high digit; don't update */
+ while (dds) { /* "add back" required */
+ i = 0; dds = 0; qhat--;
do {
- num += (long) zds[j - ny + i] + y[i];
- zds[j - ny + i] = BIGLO(num);
- num = BIGDN(num);
- } while (++i < ny);
- num--;
+ dds += (long) quods[j - ylen + i] + dends[i];
+ quods[j - ylen + i] = BIGLO(dds);
+ dds = BIGDN(dds);
+ } while (++i < ylen);
+ dds--;
+ }
+ if (mode >= 2) quods[j] = qhat; /* returning quotient */
+ } while (--j >= ylen);
+ switch (mode) {
+ case 4: /* check that remainder==0 */
+ for (j = ylen;j && !quods[j-1];--j) ; if (j) return 0;
+ case 3: /* round toward even */
+ /* Reuse dencell and dends variables to double the remainder */
+ dencell = scm_copy_big_ash1(quods, ylen, dscl);
+ dends = BDIGITS(dencell);
+ rlen = NUMDIGS(dencell);
+ if (rlen > ylen) roundup = 1;
+ else if (rlen < ylen) ;
+ else {
+ i = rlen;
+ while (i-- && (y[i]==dends[i]));
+ if (-1==i) {
+ if (0==roundup && quods[ylen] & 1) roundup = 1;
+ } else if (y[i] < dends[i]) roundup = 1;
+ }
+ case 2: /* move quotient down in quocell */
+ j = (xlen==ylen ? xlen+2 : xlen+1) - ylen;
+ for (i = 0;i < j;i++) quods[i] = quods[i+ylen];
+ ylen = i;
+ if (roundup) {
+ i = 0; dds = 1;
+ while (i < ylen) {
+ dds += quods[i];
+ quods[i++] = BIGLO(dds);
+ dds = BIGDN(dds);
+ if (!dds) break;
+ }
}
- if (modes & 2) zds[j] = qhat;
- } while (--j >= ny);
- switch (modes) {
- case 3: /* check that remainder==0 */
- for (j = ny;j && !zds[j-1];--j) ; if (j) return 0;
- case 2: /* move quotient down in z */
- j = (nx==ny ? nx+2 : nx+1) - ny;
- for (i = 0;i < j;i++) zds[i] = zds[i+ny];
- ny = i;
break;
case 1: /* subtract for modulo */
- i = 0; num = 0; j = 0;
- do {num += y[i] - zds[i];
- j = j | zds[i];
- if (num < 0) {zds[i] = num + BIGRAD; num = -1;}
- else {zds[i] = num; num = 0;}
- } while (++i < ny);
+ i = 0; dds = 0; j = 0;
+ do {dds += dends[i] - quods[i];
+ j = j | quods[i];
+ if (dds < 0) {quods[i] = dds + BIGRAD; dds = -1;}
+ else {quods[i] = dds; dds = 0;}
+ } while (++i < ylen);
if (!j) return INUM0;
case 0: /* just normalize remainder */
- if (d) divbigdig(zds, ny, d);
+ if (dscl) divbigdig(quods, ylen, dscl);
}
doadj:
- for (j = ny;j && !zds[j-1];--j) ;
+ for (j = ylen;j && !quods[j-1];--j) ;
if (j * BITSPERDIG <= sizeof(SCM)*CHAR_BIT)
- if (INUMP(z = big2inum(z, j))) return z;
- return adjbig(z, j);
+ if (INUMP(quocell = big2inum(quocell, j))) return quocell;
+ return adjbig(quocell, j);
}
#endif
@@ -2158,6 +2349,8 @@ static iproc subr2s[] = {
{s_assq, assq},
{s_assoc, assoc},
{s_quotient, lquotient},
+ /* {"rq", rq}, */
+ {s_rquotient, scm_round_quotient},
{s_remainder, lremainder},
{s_modulo, modulo},
{s_logtest, scm_logtest},