+++ /dev/null
-/* -*-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 <limits.h>
-#include <stddef.h>
-#include <stdlib.h>
-#include <string.h>
-
-#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 <stdio.h>
-#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
projects / u / mdw / catacomb / commitdiff