diff options
Diffstat (limited to 'subr.c')
| -rwxr-xr-x[-rw-r--r--] | subr.c | 471 |
1 files changed, 332 insertions, 139 deletions
@@ -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}, |