X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/f46efa79cd2bb9adc81541f1218965f85a6b2eef..fe6657c961b01ec72e9f35f4c3d96b11b31cf09c:/mpx.c diff --git a/mpx.c b/mpx.c index f1cbbd9..01264b1 100644 --- a/mpx.c +++ b/mpx.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpx.c,v 1.17 2004/03/27 00:04:46 mdw Exp $ + * $Id$ * * Low-level multiprecision arithmetic * @@ -27,64 +27,6 @@ * MA 02111-1307, USA. */ -/*----- Revision history --------------------------------------------------* - * - * $Log: mpx.c,v $ - * Revision 1.17 2004/03/27 00:04:46 mdw - * Implement efficient reduction for pleasant-looking primes. - * - * Revision 1.16 2003/05/16 09:09:24 mdw - * Fix @mp_lsl2c@. Turns out to be surprisingly tricky. - * - * Revision 1.15 2002/10/20 01:12:31 mdw - * Two's complement I/O fixes. - * - * 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 @@ -520,11 +462,11 @@ void mpx_lsl(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n) while (avl > av) { mpw t = *--avl; - *--dvl = (t >> nr) | w; + *--dvl = MPW((t >> nr) | w); w = t << nb; } - *--dvl = w; + *--dvl = MPW(w); MPX_ZERO(dv, dvl); } @@ -618,11 +560,11 @@ void mpx_lslc(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n) while (avl > av) { mpw t = *--avl; - *--dvl = (t >> nr) | w; + *--dvl = MPW((t >> nr) | w); w = t << nb; } - *--dvl = (MPW_MAX >> nr) | w; + *--dvl = MPW((MPW_MAX >> nr) | w); MPX_ONE(dv, dvl); } @@ -653,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) @@ -977,7 +919,7 @@ void mpx_usubnlsl(mpw *dv, mpw *dvl, mpw a, unsigned o) a <<= o; if (dv < dvl) { - mpd x = (mpd)*dv - (mpd)a; + mpd x = (mpd)*dv - MPW(a); *dv++ = MPW(x); if (x >> MPW_BITS) b++; @@ -1053,9 +995,7 @@ void mpx_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, */ void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m) -{ - MPX_UMULN(dv, dvl, av, avl, m); -} + { MPX_UMULN(dv, dvl, av, avl, m); } /* --- @mpx_umlan@ --- * * @@ -1070,9 +1010,7 @@ void mpx_umuln(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m) */ void mpx_umlan(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, mpw m) -{ - MPX_UMLAN(dv, dvl, av, avl, m); -} + { MPX_UMLAN(dv, dvl, av, avl, m); } /* --- @mpx_usqr@ --- * * @@ -1178,8 +1116,8 @@ void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl, unsigned b; d = dvl[-1]; - for (b = MPW_BITS / 2; b; b >>= 1) { - if (d < (MPW_MAX >> b)) { + for (b = MPW_P2; b; b >>= 1) { + if (d <= (MPW_MAX >> b)) { d <<= b; norm += b; } @@ -1479,7 +1417,7 @@ static int loadstore(dstr *v) if (!ok) dumpbits("input data", v->buf, v->len); - free(m); + xfree(m); dstr_destroy(&d); return (ok); } @@ -1517,7 +1455,7 @@ static int twocl(dstr *v) dumpbits("neg", v[1].buf, v[1].len); } - free(m); + xfree(m); dstr_destroy(&d); return (ok); @@ -1556,7 +1494,7 @@ static int twocb(dstr *v) dumpbits("neg", v[1].buf, v[1].len); } - free(m); + xfree(m); dstr_destroy(&d); return (ok); @@ -1583,7 +1521,7 @@ static int lsl(dstr *v) ok = 0; } - free(a); free(c); free(d); + xfree(a); xfree(c); xfree(d); return (ok); } @@ -1608,7 +1546,7 @@ static int lslc(dstr *v) ok = 0; } - free(a); free(c); free(d); + xfree(a); xfree(c); xfree(d); return (ok); } @@ -1633,7 +1571,7 @@ static int lsr(dstr *v) ok = 0; } - free(a); free(c); free(d); + xfree(a); xfree(c); xfree(d); return (ok); } @@ -1660,7 +1598,7 @@ static int uadd(dstr *v) ok = 0; } - free(a); free(b); free(c); free(d); + xfree(a); xfree(b); xfree(c); xfree(d); return (ok); } @@ -1687,7 +1625,7 @@ static int usub(dstr *v) ok = 0; } - free(a); free(b); free(c); free(d); + xfree(a); xfree(b); xfree(c); xfree(d); return (ok); } @@ -1714,7 +1652,7 @@ static int umul(dstr *v) ok = 0; } - free(a); free(b); free(c); free(d); + xfree(a); xfree(b); xfree(c); xfree(d); return (ok); } @@ -1738,7 +1676,7 @@ static int usqr(dstr *v) ok = 0; } - free(a); free(c); free(d); + xfree(a); xfree(c); xfree(d); return (ok); } @@ -1771,7 +1709,7 @@ static int udiv(dstr *v) ok = 0; } - free(a); free(b); free(r); free(q); free(s); free(qq); + xfree(a); xfree(b); xfree(r); xfree(q); xfree(s); xfree(qq); return (ok); }