From: mdw Date: Wed, 17 Nov 1999 18:01:11 +0000 (+0000) Subject: Split into several parts. X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/commitdiff_plain/71ec78ce88339bda7c265d6dafa4982077acc901 Split into several parts. --- diff --git a/mp.c b/mp.c deleted file mode 100644 index fb5d627..0000000 --- a/mp.c +++ /dev/null @@ -1,1252 +0,0 @@ -/* -*-c-*- - * - * $Id: mp.c,v 1.1 1999/09/03 08:41:12 mdw Exp $ - * - * Multiprecision arithmetic - * - * (c) 1998 Mark Wooding - */ - -/*----- Licensing notice --------------------------------------------------* - * - * This file is part of Catacomb. - * - * Catacomb is free software; you can redistribute it and/or modify - * 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: mp.c,v $ - * Revision 1.1 1999/09/03 08:41:12 mdw - * Initial import. - * - */ - -/*----- Header files ------------------------------------------------------*/ - -#include -#include -#include -#include - -#include "alloc.c" -#include "bits.h" -#include "exc.h" - -/*----- Important types ---------------------------------------------------*/ - -/* For now... */ - -#ifndef TEST_RIG - -typedef unsigned short mpw; -#define MPW_BITS 16 -#define MPW_MAX 0xffffu -typedef unsigned int mpd; -#define MPD_BITS 32 -#define MPD_MAX 0xffffffffu - -#else - -typedef unsigned long mpw; -static unsigned mpw_bits; -static mpw mpw_max; -typedef unsigned long mpd; - -#define MPW_BITS mpw_bits -#define MPW_MAX mpw_max - -#endif - -/*----- Data structures ---------------------------------------------------*/ - -/*----- Some basic values which are important -----------------------------*/ - -static mpw mp__v[] = { 1, 2, 3, 4, 5, 10 }; -const mp mp_std[] = { - { 0, 0, 0, 0 }, /* 0 */ - { mp__v + 0, 0, 0, 1 }, /* 1 */ - { mp__v + 1, 0, 0, 1 }, /* 2 */ - { mp__v + 2, 0, 0, 1 }, /* 3 */ - { mp__v + 3, 0, 0, 1 }, /* 4 */ - { mp__v + 4, 0, 0, 1 }, /* 5 */ - { mp__v + 5, 0, 0, 1 }, /* 10 */ - { mp__v + 0, MPF_SIGN, 0, 1 } /* -1 */ -}; - -/*----- Memory management -------------------------------------------------*/ - -/* --- @mp_create@ --- * - * - * Arguments @mp *x@ = pointer to MP head - * - * Returns: --- - * - * Use: Initializes an MP ready for use. The initial value is zero. - */ - -void mp_create(mp *x) -{ - x->v = 0; - x->f = 0; - x->sz = 0; - x->len = 0; -} - -/* --- @mp_destroy@ --- * - * - * Arguments: @mp *x@ = pointer to MP head - * - * Returns: --- - * - * Use: Releases the memory used by an MP. - */ - -void mp_destroy(mp *x) { MP_DESTROY(x); } - -/* --- @mp_resize@ --- * - * - * Arguments: @mp *x@ = pointer to MP head - * @size_t sz@ = size required (in words) - * - * Returns: --- - * - * Use: Resizes the MP so that its word vector has space for - * exactly @sz@ words. - */ - -void mp_resize(mp *x, size_t sz) -{ - if (sz == 0) - mp_destroy(x); - else { - size_t min = sz > x->sz ? x->sz : sz; - mpw *v = xmalloc(sz * sizeof(mpw)); - if (min) - memcpy(v, x->v, min * sizeof(mpw)); - MP_BURN(x); - if (x->v) - free(x->v); - if (sz > min) - memset(v + min, 0, (sz - min) * sizeof(mpw)); - x->v = v; - x->sz = sz; - } -} - -/* --- @mp_norm@ --- * - * - * Arguments: @mp *x@ = pointer to MP head - * - * Returns: --- - * - * Use: `Normalizes' an MP. Fixes the @len@ field so that it's - * correct. Assumes that @len@ is either already correct or - * too large. - */ - -void mp_norm(mp *x) { MP_NORM(x); } - -/* --- @mp_dump@ --- * - * - * Arguments: @mp *x@ = pointer to MP head - * @FILE *fp@ = pointer to stream to write on - * - * Returns: --- - * - * Use: Dumps an MP to a stream. - */ - -void mp_dump(mp *x, FILE *fp) -{ - if (x->v) { - mpw *v = x->v, *vl = v + x->len; - while (v < vl) { - fprintf(fp, "%lx", (unsigned long)*v++); - if (v < vl) fputc('-', fp); - } - } else - fputs("<0>", fp); -} - -/* --- @mp_copy@ --- * - * - * Arguments: @mp *d@ = pointer to MP head for destination - * @const mp *s@ = pointer to MP head for source - * - * Returns: --- - * - * Use: Copies an MP. - */ - -void mp_copy(mp *d, const mp *s) -{ - if (d == s) { - MP_NORM(d); - return; - } - MP_BURN(d); - MP_ENSURE(d, s->len); - memcpy(d->v, s->v, s->len * sizeof(mpw)); - if (d->f & MPF_BURN && d->sz > s->len) - memset(d->v + s->len, 0, (d->sz - s->len) * sizeof(mpw)); - d->len = s->len; - d->f = s->f; - MP_NORM(d); -} - -/* --- @mp_bits@ --- * - * - * Arguments: @mp *x@ = pointer to MP head - * - * Returns: Length of the number in bits. - * - * Use: Calculates the number of bits required to represent a number. - * The number must be normalized. - */ - -unsigned long mp_bits(mp *x) -{ - if (!x->v) - return (0); - else { - unsigned long bits = MPW_BITS * (x->len - 1); - mpw w = x->v[x->len - 1]; - while (w) { - bits++; - w >>= 1; - } - return (bits); - } -} - -/* --- @mp_octets@ --- * - * - * Arguments: @mp *x@ = pointer to MP head - * - * Returns: Length of number in octets. - * - * Use: Calculates the number of octets required to represent a - * number. The number must be normalized. - */ - -size_t mp_octets(mp *x) -{ - return ((mp_bits(x) + 7) & 7); -} - -/*----- Loading and storing as binary data --------------------------------*/ - -/* --- @mp_storel@ --- * - * - * Arguments: @mp *x@ = pointer to MP head - * @octet *p@ = pointer to octet array - * @size_t sz@ = size of octet array - * - * Returns: --- - * - * Use: Stores an MP in an octet array, least significant octet - * first. High-end octets are silently discarded if there - * isn't enough space for them. - */ - -void mp_storel(mp *x, octet *p, size_t sz) -{ - mpw *v = x->v, *vl = x->v + x->len; - mpw n, w = 0; - octet *q = p + sz; - unsigned bits = 0; - - while (p < q) { - if (bits < 8) { - n = (v >= vl) ? 0 : *v++; - *p++ = (w | n << bits) & MASK8; - w = n >> (8 - bits); - bits += MPW_BITS - 8; - } else { - *p++ = w & MASK8; - w >>= 8; - bits -= 8; - } - } -} - -/* --- @mp_loadl@ --- * - * - * Arguments: @mp *x@ = pointer to MP head - * @const octet *p@ = pointer to octet array - * @size_t sz@ = size of octet array - * - * Returns: --- - * - * Use: Loads an MP in an octet array, least significant octet - * first. - */ - -void mp_loadl(mp *x, const octet *p, size_t sz) -{ - mpw *v; - mpw w = 0; - unsigned n; - const octet *q = p + sz; - unsigned bits = 0; - - MP_BURN(x); - MP_ENSURE(x, ((sz * 8 + MPW_BITS - 1) * MPW_BITS) / MPW_BITS); - v = x->v; - - while (p < q) { - n = *p++ & MASK8; - w |= n << bits; - bits += 8; - if (bits >= MPW_BITS) { - *v++ = w & MPW_MAX; - w = n >> (MPW_BITS - bits + 8); - bits -= MPW_BITS; - } - } - *v++ = w; - x->len = v - x->v; - x->f = 0; - MP_NORM(x); -} - -/* --- @mp_storeb@ --- * - * - * Arguments: @mp *x@ = pointer to MP head - * @octet *p@ = pointer to octet array - * @size_t sz@ = size of octet array - * - * Returns: --- - * - * Use: Stores an MP in an octet array, most significant octet - * first. High-end octets are silently discarded if there - * isn't enough space for them. - */ - -void mp_storeb(mp *x, octet *p, size_t sz) -{ - mpw *v = x->v, *vl = x->v + x->len; - mpw n, w = 0; - octet *q = p + sz; - unsigned bits = 0; - - while (q > p) { - if (bits < 8) { - n = (v >= vl) ? 0 : *v++; - *--q = (w | n << bits) & MASK8; - w = n >> (8 - bits); - bits += MPW_BITS - 8; - } else { - *--q= w & MASK8; - w >>= 8; - bits -= 8; - } - } -} - -/* --- @mp_loadb@ --- * - * - * Arguments: @mp *x@ = pointer to MP head - * @const octet *p@ = pointer to octet array - * @size_t sz@ = size of octet array - * - * Returns: --- - * - * Use: Loads an MP in an octet array, most significant octet - * first. - */ - -void mp_loadb(mp *x, const octet *p, size_t sz) -{ - mpw *v; - mpw w = 0; - unsigned n; - const octet *q = p + sz; - unsigned bits = 0; - - MP_BURN(x); - MP_ENSURE(x, ((sz * 8 + MPW_BITS - 1) * MPW_BITS) / MPW_BITS); - v = x->v; - - while (q > p) { - n = *--q & MASK8; - w |= n << bits; - bits += 8; - if (bits >= MPW_BITS) { - *v++ = w & MPW_MAX; - w = n >> (MPW_BITS - bits + 8); - bits -= MPW_BITS; - } - } - *v++ = w; - x->len = v - x->v; - x->f = 0; - MP_NORM(x); -} - -/*----- Iterating through bits --------------------------------------------*/ - -/* --- @mp_mkbitscan@ --- * - * - * Arguments: @mp_bitscan *sc@ = pointer to bitscan object - * @const mp *x@ = pointer to MP head - * - * Returns: --- - * - * Use: Initializes a bitscan object. - */ - -void mp_mkbitscan(mp_bitscan *sc, const mp *x) -{ - sc->x = x; - sc->bits = 0; - sc->i = 0; -} - -/* --- @mp_bstep@ --- * - * - * Arguments: @mp_bitscan *sc@ = pointer to bitscanner object - * - * Returns: Nonzero if there is another bit to read. - * - * Use: Steps on to the next bit, and tells the caller whether one - * exists. - */ - -int mp_bstep(mp_bitscan *sc) -{ - if (sc->bits) { - sc->w >>= 1; - sc->bits--; - return (1); - } - if (sc->i >= sc->x->len) - return (0); - sc->w = sc->x->v[sc->i++]; - sc->bits = MPW_BITS - 1; - return (1); -} - -/* --- @mp_bit@ --- * - * - * Arguments: @const mp_bitscan *sc@ = pointer to bitscanner - * - * Returns: Current bit value. - * - * Use: Returns the value of the current bit. - */ - -int mp_bit(const mp_bitscan *sc) { return (MP_BIT(sc)); } - -/*----- Shifting ----------------------------------------------------------*/ - -/* --- @mp_lsl@ --- * - * - * Arguments: @mp *d@ = pointer to MP head of destination - * @const mp *x@ = pointer to MP head of source - * @size_t n@ = number of bits to shift - * - * Returns: --- - * - * Use: Shifts a number left by a given number of bit positions. - */ - -void mp_lsl(mp *d, const mp *x, size_t n) -{ - size_t nw = n / MPW_BITS; - unsigned nb = n % MPW_BITS; - unsigned nr = MPW_BITS - nb; - size_t req; - - /* --- Trivial special case --- */ - - if (n == 0 || mp_ucmp(x, MP_ZERO) == 0) { - mp_copy(d, x); - goto copied; - } - - /* --- Decide on the about of memory needed in the destination --- */ - - req = x->len + nw; - if (nb && (x->v[x->len - 1] >> nr) != 0) - req++; - if (d != x) - MP_BURN(d); - MP_ENSURE(d, req); - - /* --- Handle single bit shifting --- */ - - if (n == 1) { - mpw *v = d->v; - const mpw *vx = x->v; - const mpw *vl; - mpw w = 0; - - vl = vx + x->len; - while (vx < vl) { - mpw t = *vx++; - *v++ = ((t << 1) | w) & MPW_MAX; - w = t >> (MPW_BITS - 1); - } - if (v < d->v + req) - *v++ = w & MPW_MAX; - goto done; - } - - /* --- Handle shifts by a multiple of the word size --- * - * - * This would be easy, except that C is irritating. Shifting an integer - * by an amount equal to the type's width yields undefined behaviour; - * in particular, under Intel it's a no-op. - */ - - if (!nb) { - mpw *v = d->v; - const mpw *vx = x->v; - memmove(v + nw, vx, (req - nw) * sizeof(mpw)); - memset(v, 0, nw * sizeof(mpw)); - goto done; - } - - /* --- Now do the difficult version --- */ - - { - mpw *v = d->v + req; - const mpw *vx = x->v + x->len; - const mpw *vl; - mpw w = *--vx; - - /* --- First shift the data over --- */ - - if (req > x->len + nw) - *--v = w >> nr; - vl = x->v; - while (vx > vl) { - mpw t = *--vx; - *--v = ((w << nb) | (t >> nr)) & MPW_MAX; - w = t; - } - - /* --- Deal with tail-end data --- */ - - vl = d->v; - if (v > vl) { - *--v = w << nb; - while (v > vl) - *--v = 0; - } - } - - /* --- Common end code --- */ - -done: - d->len = req; - d->f = x->f; -copied: - MP_NORM(d); -} - -/* --- @mp_lsr@ --- * - * - * Arguments: @mp *d@ = pointer to MP head of destination - * @const mp *x@ = pointer to MP head of source - * @size_t n@ = number of bits to shift - * - * Returns: --- - * - * Use: Shifts a number right by a given number of bit positions. - */ - -void mp_lsr(mp *d, const mp *x, size_t n) -{ - size_t nw = n / MPW_BITS; - unsigned nb = n % MPW_BITS; - unsigned nr = MPW_BITS - nb; - size_t req; - - /* --- Trivial special case --- */ - - if (n == 0 || mp_ucmp(x, MP_ZERO) == 0) { - mp_copy(d, x); - goto copied; - } - - /* --- Decide on the about of memory needed in the destination --- */ - - req = x->len - nw; - if ((x->v[x->len - 1] >> nb) == 0) - req--; - if (d != x) - MP_BURN(d); - MP_ENSURE(d, req); - - /* --- Handle single bit shifting --- */ - - if (n == 1) { - mpw *v = d->v + req; - const mpw *vx = x->v + x->len; - const mpw *vl; - mpw w = *--vx; - - if (req == x->len) - *--v = (w >> 1) & MPW_MAX; - vl = x->v; - while (vx > vl) { - mpw t = *--vx; - *--v = (w << (MPW_BITS - 1) | (t >> 1)) & MPW_MAX; - w = t; - } - goto done; - } - - /* --- Handle shifts by a multiple of the word size --- * - * - * This would be easy, except that C is irritating. Shifting an integer - * by an amount equal to the type's width yields undefined behaviour; - * in particular, under Intel it's a no-op. - */ - - if (!nb) { - mpw *v = d->v; - const mpw *vx = x->v; - memmove(v, vx + nw, (x->len - nw) * sizeof(mpw)); - goto done; - } - - /* --- Now do the difficult version --- */ - - { - mpw *v = d->v; - const mpw *vx = x->v + nw; - const mpw *vl; - mpw w = *vx++; - - /* --- First shift the data over --- */ - - vl = x->v + x->len; - while (vx < vl) { - mpw t = *vx++; - *v++ = ((w >> nb) | (t << nr)) & MPW_MAX; - w = t; - } - if (req == x->len - nw) { - *v++ = (w >> nb) & MPW_MAX; - } - } - - /* --- Common end code --- */ - -done: - d->len = req; - d->f = x->f; -copied: - MP_NORM(d); -} - -/*----- Adding and subtracting --------------------------------------------*/ - -/* --- @mp_uadd@ --- * - * - * Arguments: @const mp *d@ = pointers to MP head of destination - * @const mp *x, *y@ = pointers to MP heads of operands - * - * Returns: --- - * - * Use: Performs unsigned MP addition. - */ - -void mp_uadd(mp *d, const mp *x, const mp *y) -{ - mpd c; - mpw *v; - const mpw *vx, *vy; - const mpw *vxl, *vyl; - - /* --- Some trivial initialization --- */ - - if (d != x && d != y) - MP_BURN(d); - MP_ENSURE(d, (x->len > y->len ? x->len : y->len) + 1); - vx = x->v; vxl = vx + x->len; - vy = y->v; vyl = vy + y->len; - v = d->v; - c = 0; - - /* --- Start on the work --- */ - - while (vx < vxl || vy < vyl) { - if (vx < vxl) c += *vx++; - if (vy < vyl) c += *vy++; - *v++ = c & MPW_MAX; - c >>= MPW_BITS; - } - if (c) - *v++ = c & MPW_MAX; - - /* --- Tidy up --- */ - - d->len = v - d->v; - d->f = 0; - if (x->f & MPF_BURN || y->f & MPF_BURN) - MP_PARANOID(d); - MP_NORM(d); -} - -/* --- @mp_usub@ --- * - * - * Arguments: @const mp *d@ = pointers to MP head of destination - * @const mp *x, *y@ = pointers to MP heads of operands - * - * Returns: --- - * - * Use: Performs unsigned MP subtraction. - */ - -void mp_usub(mp *d, const mp *x, const mp *y) -{ - mpd c; - mpw *v; - const mpw *vx, *vy; - const mpw *vxl, *vyl; - - /* --- Some trivial initialization --- */ - - if (d != x && d != y) - MP_BURN(d); - MP_ENSURE(d, x->len); - vx = x->v; vxl = vx + x->len; - vy = y->v; vyl = vy + y->len; - v = d->v; - c = 0; - - /* --- Start on the work --- */ - - while (vx < vxl) { - c += *vx++; - if (vy < vyl) c -= *vy++; - *v++ = c & MPW_MAX; - if (c & ~MPW_MAX) - c = ~0ul; - else - c = 0; - } - - /* --- Tidy up --- */ - - d->len = v - d->v; - d->f = 0; - if (x->f & MPF_BURN || y->f & MPF_BURN) - MP_PARANOID(d); - MP_NORM(d); -} - -/* --- @mp_ucmp@ --- * - * - * Arguments: @const mp *x, *y@ = pointers to MP heads of operands - * - * Returns: Less than, equal to, or greater than zero. - * - * Use: Performs unsigned MP comparison. - */ - -int mp_ucmp(const mp *x, const mp *y) -{ - int i; - mpw a, b; - - /* --- Decide which to examine --- */ - - if (x->len > y->len) - i = x->len; - else - i = y->len; - - /* --- Loop through the data --- */ - - while (i > 0) { - i--; - a = (i < x->len ? x->v[i] : 0); - b = (i < y->len ? y->v[i] : 0); - if (a < b) - return (-1); - else if (a > b) - return (1); - } - - /* --- Finished --- */ - - return (0); -} - -/*----- Multiplying and dividing ------------------------------------------*/ - -/* --- @mp_umul@ --- * - * - * Arguments: @mp *d@ = pointer to MP head of destination - * @const mp *x, *y@ = pointes to MP heads of operands - * - * Returns: --- - * - * Use: Performs unsigned MP multiplication. - */ - -void mp_umul(mp *d, const mp *x, const mp *y) -{ - mpd c, f; - mpw *v = 0, *vbase; - const mpw *vx, *vy; - const mpw *vxl, *vyl; - - /* --- Check for special cases --- */ - - if (mp_ucmp(x, MP_ZERO) == 0 || mp_ucmp(y, MP_ZERO) == 0) { - mp_copy(d, MP_ZERO); - return; - } - - /* --- Some trivial initialization --- */ - - MP_BURN(d); - MP_ENSURE(d, x->len + y->len); - vxl = x->v + x->len; - vyl = y->v + y->len; - vbase = d->v; - memset(v, 0, (x->len + y->len) * sizeof(mpw)); - - /* --- Main loop --- */ - - for (vx = x->v; vx < vxl; vx++) { - c = 0; - f = *vx; - v = vbase++; - for (vy = y->v; vy < vyl; vy++) { - c += *v + f * *vy; - *v++ = c & MPW_MAX; - c >>= MPW_BITS; - } - *v++ = c & MPW_MAX; - } - - /* --- Tidying up --- */ - - d->len = v - d->v; - d->f = 0; - if (x->f & MPF_BURN || y->f & MPF_BURN) - MP_PARANOID(d); - MP_NORM(d); -} - -/* --- @mp_udiv@ --- * - * - * Arguments: @mp *q, *r@ = pointers to MP heads for quotient, remainder - * @const mp *x, *y@ = pointers to MP heads for operands - * - * Returns: --- - * - * Use: Performs unsigned MP division. - */ - -void mp_udiv(mp *q, mp *r, const mp *x, const mp *y) -{ - size_t n = x->len, t = y->len; - const mpw *vy, *ovy; - mpw *vx, *vq, *vxl, *va, *vb; - mpw yt, ytt; - mpd yd; - mp nx = MP_INIT, ny = MP_INIT, nq = MP_INIT; - mp tt = MP_INIT; - unsigned int shift; - - /* --- Fiddle with some pointers --- */ - - if (!r) - r = &nx; - if (!q) - q = &nq; - - /* --- Find the top word in @y@ --- */ - - vy = y->v; - for (;;) { - if (!t) - THROW(EXC_FAIL, "mp_udiv detected divide-by-zero"); - if ((yt = vy[t - 1]) != 0) - break; - t--; - } - - /* --- Check for a zero dividend --- */ - - if (mp_ucmp(x, MP_ZERO)) { - if (q) mp_copy(q, MP_ZERO); - if (r) mp_copy(r, MP_ZERO); - return; - } - - /* --- Test for some other trivial cases --- */ - - { - int i = mp_ucmp(x, y); - if (i < 0) { - if (r) mp_copy(r, x); - if (q) mp_copy(q, MP_ZERO); - return; - } else if (i == 0) { - if (r) mp_copy(r, MP_ZERO); - if (q) mp_copy(q, MP_ONE); - return; - } - } - - /* --- Normalize the divisor --- * - * - * Cheat. The original algorithm wants two-word values at least, so I just - * shift everything up by a word if necessary. - */ - - shift = 0; - while (yt < MPW_MAX / 2) { - yt <<= 1; - shift++; - } - if (t <= 1) - shift += MPW_BITS; - mp_lsl(r, x, shift); - mp_lsl(&ny, y, shift); - - n = r->len; - t = ny.len; - ytt = ny.v[t - 2]; - yd = ((mpd)yt << MPW_BITS) | (mpd)ytt; - - /* --- Initialize the quotient --- */ - - MP_ENSURE(q, n - t + 1); - vq = q->v + n - t; - memset(vq, 0, (n - t) * sizeof(mpw)); - - /* --- Shift the divisor up to match the dividend --- */ - - mp_lsl(&ny, (n - t) * MPW_BITS); - ovy = ny.v; - - /* --- Get the most significant quotient digit --- * - * - * Because of the normalization, this should only happen once. - */ - - while (mp_ucmp(x, y) >= 0) { - mp_usub(x, x, y); - (*vq)++; - } - - /* --- Now do the main loop --- */ - - vx = x->v + n; - vxl = x->v + t; - - while (vx > vxl) { - mpw xi, xii, xiii; - mpd qi; - mpd xd, yhi, ylo; - - /* --- Fetch the top words from @x@ --- */ - - vx--; - xi = vx[0]; - xii = vx[-1]; - xiii = vx[-2]; - xd = ((mpd)xi << MPW_BITS) | (mpd)xii; - - /* --- Get an approximation for @qi@ --- */ - - if (xi == yt) - qi = MPW_MAX; - else - qi = xd / yt; - - /* --- Work out how close the approximation is --- * - * - * This is more than a little ugly. - */ - - yhi = (mpd)yt * qi; - ylo = (mpd)ytt * qi; - yhi += ylo >> MPW_BITS; - ylo &= MPW_MAX; - - /* --- Now fix the approximation --- * - * - * Again, the normalization helps; this is never done more than twice. - */ - - while (!(yhi <= xd && ylo <= xiii)) { - qi--; - yhi -= yt; - if (ylo < ytt) - yhi--; - ylo = (ylo - ytt) & MPW_MAX; - } - - /* --- Subtract off a goodly big chunk --- */ - - -} - -/*----- Test rig ----------------------------------------------------------*/ - -#ifdef TEST_RIG - -#include -#include "dstr.h" -#include "testrig.h" - -/* --- Loading and storing numbers --- * - * - * The most reliable way I can think of for doing this is to convert both - * numbers (the MP and the original hexgorp) into binary text, and compare. - * This ought to be reliable, and it's sufficiently different from what the - * actual load and store routines are doing not to have any common bugs. - */ - -static void mpbindump(const mp *x, char *p) -{ - mp_bitscan sc; - for (mp_mkbitscan(&sc, x); mp_bstep(&sc); ) - *p++ = '0' + MP_BIT(&sc); - *p++ = 0; -} - -static void bindumpl(const octet *q, size_t sz, char *p) -{ - unsigned w; - int bits = 0; - - while (sz) { - w = *q++; - for (bits = 8; bits > 0; bits--) { - *p++ = '0' + (w & 1); - w >>= 1; - } - sz--; - } - *p++ = 0; -} - -static void bindumpb(const octet *q, size_t sz, char *p) -{ - unsigned w; - int bits = 0; - - q += sz; - while (sz) { - w = *--q; - for (bits = 8; bits > 0; bits--) { - *p++ = '0' + (w & 1); - w >>= 1; - } - sz--; - } - *p++ = 0; -} - -static int bincmp(const char *p, const char *q) -{ - for (;;) { - if (!*p && !*q) - return (0); - else if (!*p && *q == '0') - q++; - else if (*p == '0' && !*q) - p++; - else if (*p == *q) - p++, q++; - else - return (-1); - } -} - -static int lscheck(const char *reason, int w, dstr *d, mp *x, - const char *bufa, const char *bufb) -{ - if (bincmp(bufa, bufb)) { - printf("\nmismatch in %s (width = %i):\n" - "\tvalue = ", reason, w); - type_hex.dump(d, stdout); - printf("\n\tmp = "); - mp_dump(x, stdout); - printf("\n\texpected = %s\n" - "\tcalculated = %s\n", - bufb, bufa); - return (1); - } - return (0); -} - -static int loadstore(dstr *d) -{ - octet *p = (octet *)d[0].buf; - size_t sz = d[0].len; - char bufa[1024], bufb[1024]; - octet bufc[64]; - mp x; - int i; - int ok = 1; - - mp_wmax = 0x7ful; - for (i = 8; ok && i <= 32; i++) { - mpw_bits = i; - mpw_max = (mp_wmax << 1) + 1; - - mp_create(&x); mp_loadl(&x, p, sz); - mpbindump(&x, bufa); bindumpl(p, sz, bufb); - if (lscheck("mp_loadl", i, d, &x, bufa, bufb)) ok = 0; - mp_storeb(&x, bufc, sizeof(bufc)); - bindumpb(bufc, sizeof(bufc), bufa); - if (lscheck("mp_storeb", i, d, &x, bufa, bufb)) ok = 0; - mp_destroy(&x); - - mp_create(&x); mp_loadb(&x, p, sz); - mpbindump(&x, bufa); bindumpb(p, sz, bufb); - if (lscheck("mp_loadb", i, d, &x, bufa, bufb)) ok = 0; - mp_storel(&x, bufc, sizeof(bufc)); - bindumpl(bufc, sizeof(bufc), bufa); - if (lscheck("mp_storel", i, d, &x, bufa, bufb)) ok = 0; - mp_destroy(&x); - } - - mpw_max = 0xffff; - mpw_bits = 16; - return (ok); -} - -/* --- Comparison test --- */ - -static int compare(dstr *d) -{ - mp x, y; - int r, s; - int ok = 1; - - mp_create(&x); mp_create(&y); - mp_loadb(&x, (octet *)d[0].buf, d[0].len); - mp_loadb(&y, (octet *)d[1].buf, d[1].len); - r = *(int *)d[2].buf; - s = mp_ucmp(&x, &y); - - if (r != s) { - ok = 0; - printf("\nfailed compare:" - "\n\tx = "); - type_hex.dump(&d[0], stdout); - printf("\n\ty = "); - type_hex.dump(&d[1], stdout); - printf("\n\texpected %i; found %i\n", r, s); - } - - return (ok); -} - -/* --- Addition and subtraction test --- */ - -static int addsub(dstr *d) -{ - mp t, x, y, z; - int ok = 1; - - mp_create(&t); mp_create(&x); mp_create(&y); mp_create(&z); - mp_loadb(&x, (octet *)d[0].buf, d[0].len); - mp_loadb(&y, (octet *)d[1].buf, d[1].len); - mp_loadb(&z, (octet *)d[2].buf, d[2].len); - - mp_uadd(&t, &x, &y); - if (mp_ucmp(&t, &z)) { - ok = 0; - printf("\nfailed add:"); - printf("\n\tx = "); type_hex.dump(&d[0], stdout); - printf("\n\ty = "); type_hex.dump(&d[1], stdout); - printf("\n\tz = "); type_hex.dump(&d[2], stdout); - fputc('\n', stdout); - } - - mp_usub(&t, &z, &x); - if (mp_ucmp(&t, &y)) { - ok = 0; - printf("\nfailed subtract:"); - printf("\n\tz = "); type_hex.dump(&d[2], stdout); - printf("\n\tx = "); type_hex.dump(&d[0], stdout); - printf("\n\ty = "); type_hex.dump(&d[1], stdout); - fputc('\n', stdout); - } - - return (ok); -} - -/* --- Shifting --- */ - -static int shift(dstr *d) -{ - char bufa[1024], bufb[1024]; - mp x, y; - int n = *(int *)d[1].buf; - char *p; - int i; - int ok = 1; - - mp_create(&x); - mp_create(&y); - - mp_loadb(&x, (octet *)d[0].buf, d[0].len); - p = bufa; - for (i = 0; i < n; i++) - *p++ = '0'; - mpbindump(&x, p); - - mp_lsl(&y, &x, n); - mpbindump(&y, bufb); - if (lscheck("lsl", n, d, &x, bufb, bufa)) - ok = 0; - - mp_lsr(&y, &x, n); - mpbindump(&y, bufb); - if (lscheck("lsr", n, d, &x, bufb, bufa + 2 * n)) - ok = 0; - - return (ok); -} - -/* --- Test driver stub --- */ - -static test_chunk mp__defs[] = { - { "mp-loadstore", loadstore, { &type_hex, 0 } }, - { "mp-cmp", compare, { &type_hex, &type_hex, &type_int, 0 } }, - { "mp-addsub", addsub, { &type_hex, &type_hex, &type_hex, 0 } }, - { "mp-shift", shift, { &type_hex, &type_int, 0 } }, - { 0, 0, { 0 } } -}; - -int main(int argc, char *argv[]) -{ - test_run(argc, argv, mp__defs, SRCDIR"/tests/mp"); - return (0); -} - -#endif - -/*----- That's all, folks -------------------------------------------------*/ - -#endif