X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/4f29a7323a664994f633a7bcd5a1afc2914156fa..c65df27983057ec76ed0e72bb370f9a5ae7dad28:/mpx.c diff --git a/mpx.c b/mpx.c index d7ea70a..e122760 100644 --- a/mpx.c +++ b/mpx.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpx.c,v 1.14 2002/10/19 18:55:08 mdw Exp $ + * $Id: mpx.c,v 1.20 2004/04/08 01:36:15 mdw Exp $ * * Low-level multiprecision arithmetic * @@ -27,55 +27,6 @@ * MA 02111-1307, USA. */ -/*----- Revision history --------------------------------------------------* - * - * $Log: mpx.c,v $ - * Revision 1.14 2002/10/19 18:55:08 mdw - * Fix overflows in shift primitives. - * - * Revision 1.13 2002/10/19 17:56:50 mdw - * Fix bit operations. Test them (a bit) better. - * - * Revision 1.12 2002/10/06 22:52:50 mdw - * Pile of changes for supporting two's complement properly. - * - * Revision 1.11 2001/04/03 19:36:05 mdw - * Add some simple bitwise operations so that Perl can use them. - * - * Revision 1.10 2000/10/08 12:06:12 mdw - * Provide @mpx_ueq@ for rapidly testing equality of two integers. - * - * Revision 1.9 2000/06/26 07:52:50 mdw - * Portability fix for the bug fix. - * - * Revision 1.8 2000/06/25 12:59:02 mdw - * (mpx_udiv): Fix bug in quotient digit estimation. - * - * Revision 1.7 1999/12/22 15:49:07 mdw - * New function for division by a small integer. - * - * Revision 1.6 1999/11/20 22:43:44 mdw - * Integrate testing for MPX routines. - * - * Revision 1.5 1999/11/20 22:23:27 mdw - * Add function versions of some low-level macros with wider use. - * - * Revision 1.4 1999/11/17 18:04:09 mdw - * Add two's-complement functionality. Improve mpx_udiv a little by - * performing the multiplication of the divisor by q with the subtraction - * from r. - * - * Revision 1.3 1999/11/13 01:57:31 mdw - * Remove stray debugging code. - * - * Revision 1.2 1999/11/13 01:50:59 mdw - * Multiprecision routines finished and tested. - * - * Revision 1.1 1999/09/03 08:41:12 mdw - * Initial import. - * - */ - /*----- Header files ------------------------------------------------------*/ #include @@ -281,12 +232,12 @@ void mpx_storel2cn(const mpw *v, const mpw *vl, void *pp, size_t sz) bits -= 8; } b = U8(~b + c); - c = !b; + c = c && !b; *p++ = b; } while (p < q) { b = U8(~b + c); - c = !b; + c = c && !b; *p++ = b; b = 0; } @@ -318,7 +269,7 @@ void mpx_loadl2cn(mpw *v, mpw *vl, const void *pp, size_t sz) return; while (p < q) { n = U8(~(*p++) + c); - c = !n; + c = c && !n; w |= n << bits; bits += 8; if (bits >= MPW_BITS) { @@ -371,12 +322,12 @@ void mpx_storeb2cn(const mpw *v, const mpw *vl, void *pp, size_t sz) bits -= 8; } b = U8(~b + c); - c = !b; + c = c && !b; *--q = b; } while (q > p) { b = ~b + c; - c = !(b & 0xff); + c = c && !(b & 0xff); *--q = b; b = 0; } @@ -408,7 +359,7 @@ void mpx_loadb2cn(mpw *v, mpw *vl, const void *pp, size_t sz) return; while (q > p) { n = U8(~(*--q) + c); - c = !n; + c = c && !n; w |= n << bits; bits += 8; if (bits >= MPW_BITS) { @@ -522,6 +473,104 @@ void mpx_lsl(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n) done:; } +/* --- @mpx_lslc@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = source vector base and limit + * @size_t n@ = number of bit positions to shift by + * + * Returns: --- + * + * Use: Performs a logical shift left operation on an integer, only + * it fills in the bits with ones instead of zeroes. + */ + +void mpx_lslc(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n) +{ + size_t nw; + unsigned nb; + + /* --- Trivial special case --- */ + + if (n == 0) + MPX_COPY(dv, dvl, av, avl); + + /* --- Single bit shifting --- */ + + else if (n == 1) { + mpw w = 1; + while (av < avl) { + mpw t; + if (dv >= dvl) + goto done; + t = *av++; + *dv++ = MPW((t << 1) | w); + w = t >> (MPW_BITS - 1); + } + if (dv >= dvl) + goto done; + *dv++ = MPW(w); + MPX_ZERO(dv, dvl); + goto done; + } + + /* --- Break out word and bit shifts for more sophisticated work --- */ + + nw = n / MPW_BITS; + nb = n % MPW_BITS; + + /* --- Handle a shift by a multiple of the word size --- */ + + if (nb == 0) { + if (nw >= dvl - dv) + MPX_ONE(dv, dvl); + else { + MPX_COPY(dv + nw, dvl, av, avl); + MPX_ONE(dv, dv + nw); + } + } + + /* --- And finally the difficult case --- * + * + * This is a little convoluted, because I have to start from the end and + * work backwards to avoid overwriting the source, if they're both the same + * block of memory. + */ + + else { + mpw w; + size_t nr = MPW_BITS - nb; + size_t dvn = dvl - dv; + size_t avn = avl - av; + + if (dvn <= nw) { + MPX_ONE(dv, dvl); + goto done; + } + + if (dvn > avn + nw) { + size_t off = avn + nw + 1; + MPX_ZERO(dv + off, dvl); + dvl = dv + off; + w = 0; + } else { + avl = av + dvn - nw; + w = *--avl << nb; + } + + while (avl > av) { + mpw t = *--avl; + *--dvl = (t >> nr) | w; + w = t << nb; + } + + *--dvl = (MPW_MAX >> nr) | w; + MPX_ONE(dv, dvl); + } + +done:; +} + /* --- @mpx_lsr@ --- * * * Arguments: @mpw *dv, *dvl@ = destination vector base and limit @@ -546,7 +595,7 @@ void mpx_lsr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n) /* --- Single bit shifting --- */ else if (n == 1) { - mpw w = *av++ >> 1; + mpw w = av < avl ? *av++ >> 1 : 0; while (av < avl) { mpw t; if (dv >= dvl) @@ -773,6 +822,30 @@ void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, void mpx_uaddn(mpw *dv, mpw *dvl, mpw n) { MPX_UADDN(dv, dvl, n); } +/* --- @mpx_uaddnlsl@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination and first argument vector + * @mpw a@ = second argument + * @unsigned o@ = offset in bits + * + * Returns: --- + * + * Use: Computes %$d + 2^o a$%. If the result overflows then + * high-order bits are discarded, as usual. We must have + * @0 < o < MPW_BITS@. + */ + +void mpx_uaddnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o) +{ + mpd x = (mpd)a << o; + + while (x && dv < dvl) { + x += *dv; + *dv++ = MPW(x); + x >>= MPW_BITS; + } +} + /* --- @mpx_usub@ --- * * * Arguments: @mpw *dv, *dvl@ = destination vector base and limit @@ -827,6 +900,33 @@ void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, void mpx_usubn(mpw *dv, mpw *dvl, mpw n) { MPX_USUBN(dv, dvl, n); } +/* --- @mpx_uaddnlsl@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination and first argument vector + * @mpw a@ = second argument + * @unsigned o@ = offset in bits + * + * Returns: --- + * + * Use: Computes %$d + 2^o a$%. If the result overflows then + * high-order bits are discarded, as usual. We must have + * @0 < o < MPW_BITS@. + */ + +void mpx_usubnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o) +{ + mpw b = a >> (MPW_BITS - o); + a <<= o; + + if (dv < dvl) { + mpd x = (mpd)*dv - (mpd)a; + *dv++ = MPW(x); + if (x >> MPW_BITS) + b++; + MPX_USUBN(dv, dvl, b); + } +} + /* --- @mpx_umul@ --- * * * Arguments: @mpw *dv, *dvl@ = destination vector base and limit @@ -1021,7 +1121,7 @@ void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl, d = dvl[-1]; for (b = MPW_BITS / 2; b; b >>= 1) { - if (d < (MPW_MAX >> b)) { + if (d <= (MPW_MAX >> b)) { d <<= b; norm += b; } @@ -1429,6 +1529,31 @@ static int lsl(dstr *v) return (ok); } +static int lslc(dstr *v) +{ + mpw *a, *al; + int n = *(int *)v[1].buf; + mpw *c, *cl; + mpw *d, *dl; + int ok = 1; + + LOAD(a, al, &v[0]); + LOAD(c, cl, &v[2]); + ALLOC(d, dl, al - a + (n + MPW_BITS - 1) / MPW_BITS); + + mpx_lslc(d, dl, a, al, n); + if (!mpx_ueq(d, dl, c, cl)) { + fprintf(stderr, "\n*** lslc(%i) failed\n", n); + dumpmp(" a", a, al); + dumpmp("expected", c, cl); + dumpmp(" result", d, dl); + ok = 0; + } + + free(a); free(c); free(d); + return (ok); +} + static int lsr(dstr *v) { mpw *a, *al; @@ -1597,6 +1722,7 @@ static test_chunk defs[] = { { "2cl", twocl, { &type_hex, &type_hex, } }, { "2cb", twocb, { &type_hex, &type_hex, } }, { "lsl", lsl, { &type_hex, &type_int, &type_hex, 0 } }, + { "lslc", lslc, { &type_hex, &type_int, &type_hex, 0 } }, { "lsr", lsr, { &type_hex, &type_int, &type_hex, 0 } }, { "uadd", uadd, { &type_hex, &type_hex, &type_hex, 0 } }, { "usub", usub, { &type_hex, &type_hex, &type_hex, 0 } },