/* -*-c-*-
*
- * $Id: mpx.c,v 1.17 2004/03/27 00:04:46 mdw Exp $
+ * $Id$
*
* Low-level multiprecision arithmetic
*
* (c) 1999 Straylight/Edgeware
*/
-/*----- Licensing notice --------------------------------------------------*
+/*----- Licensing notice --------------------------------------------------*
*
* This file is part of Catacomb.
*
* it under the terms of the GNU Library General Public License as
* published by the Free Software Foundation; either version 2 of the
* License, or (at your option) any later version.
- *
+ *
* Catacomb is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Library General Public License for more details.
- *
+ *
* You should have received a copy of the GNU Library General Public
* License along with Catacomb; if not, write to the Free
* Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
* 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 <assert.h>
}
/* --- Break out word and bit shifts for more sophisticated work --- */
-
+
nw = n / MPW_BITS;
nb = n % MPW_BITS;
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);
}
}
/* --- Break out word and bit shifts for more sophisticated work --- */
-
+
nw = n / MPW_BITS;
nb = n % MPW_BITS;
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);
}
/* --- 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)
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++;
MPX_SHRINK(bv, bvl);
/* --- Deal with a multiply by zero --- */
-
+
if (bv == bvl) {
MPX_ZERO(dv, dvl);
return;
*/
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@ --- *
*
*/
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@ --- *
*
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;
}
} while (0)
#define MAX(x, y) ((x) > (y) ? (x) : (y))
-
+
static void dumpbits(const char *msg, const void *pp, size_t sz)
{
const octet *p = pp;
if (!ok)
dumpbits("input data", v->buf, v->len);
- free(m);
+ xfree(m);
dstr_destroy(&d);
return (ok);
}
dumpbits("neg", v[1].buf, v[1].len);
}
- free(m);
+ xfree(m);
dstr_destroy(&d);
return (ok);
dumpbits("neg", v[1].buf, v[1].len);
}
- free(m);
+ xfree(m);
dstr_destroy(&d);
return (ok);
mpx_lsl(d, dl, a, al, n);
if (!mpx_ueq(d, dl, c, cl)) {
fprintf(stderr, "\n*** lsl(%i) failed\n", n);
- dumpmp(" a", a, al);
+ dumpmp(" a", a, al);
dumpmp("expected", c, cl);
dumpmp(" result", d, dl);
ok = 0;
}
- free(a); free(c); free(d);
+ xfree(a); xfree(c); xfree(d);
return (ok);
}
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(" a", a, al);
dumpmp("expected", c, cl);
dumpmp(" result", d, dl);
ok = 0;
}
- free(a); free(c); free(d);
+ xfree(a); xfree(c); xfree(d);
return (ok);
}
mpx_lsr(d, dl, a, al, n);
if (!mpx_ueq(d, dl, c, cl)) {
fprintf(stderr, "\n*** lsr(%i) failed\n", n);
- dumpmp(" a", a, al);
+ dumpmp(" a", a, al);
dumpmp("expected", c, cl);
dumpmp(" result", d, dl);
ok = 0;
}
- free(a); free(c); free(d);
+ xfree(a); xfree(c); xfree(d);
return (ok);
}
mpx_uadd(d, dl, a, al, b, bl);
if (!mpx_ueq(d, dl, c, cl)) {
fprintf(stderr, "\n*** uadd failed\n");
- dumpmp(" a", a, al);
- dumpmp(" b", b, bl);
+ dumpmp(" a", a, al);
+ dumpmp(" b", b, bl);
dumpmp("expected", c, cl);
dumpmp(" result", d, dl);
ok = 0;
}
- free(a); free(b); free(c); free(d);
+ xfree(a); xfree(b); xfree(c); xfree(d);
return (ok);
}
mpx_usub(d, dl, a, al, b, bl);
if (!mpx_ueq(d, dl, c, cl)) {
fprintf(stderr, "\n*** usub failed\n");
- dumpmp(" a", a, al);
- dumpmp(" b", b, bl);
+ dumpmp(" a", a, al);
+ dumpmp(" b", b, bl);
dumpmp("expected", c, cl);
dumpmp(" result", d, dl);
ok = 0;
}
- free(a); free(b); free(c); free(d);
+ xfree(a); xfree(b); xfree(c); xfree(d);
return (ok);
}
mpx_umul(d, dl, a, al, b, bl);
if (!mpx_ueq(d, dl, c, cl)) {
fprintf(stderr, "\n*** umul failed\n");
- dumpmp(" a", a, al);
- dumpmp(" b", b, bl);
+ dumpmp(" a", a, al);
+ dumpmp(" b", b, bl);
dumpmp("expected", c, cl);
dumpmp(" result", d, dl);
ok = 0;
}
- free(a); free(b); free(c); free(d);
+ xfree(a); xfree(b); xfree(c); xfree(d);
return (ok);
}
mpx_usqr(d, dl, a, al);
if (!mpx_ueq(d, dl, c, cl)) {
fprintf(stderr, "\n*** usqr failed\n");
- dumpmp(" a", a, al);
+ dumpmp(" a", a, al);
dumpmp("expected", c, cl);
dumpmp(" result", d, dl);
ok = 0;
}
- free(a); free(c); free(d);
+ xfree(a); xfree(c); xfree(d);
return (ok);
}
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);
}