Fix bit operations. Test them (a bit) better.

[u/mdw/catacomb] / mpx.h

/* -*-c-*-
 *
 * $Id: mpx.h,v 1.15 2002/10/19 17:56:50 mdw Exp $
 *
 * Low level multiprecision arithmetic
 *
 * (c) 1999 Straylight/Edgeware
 */

/*----- 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: mpx.h,v $
 * Revision 1.15  2002/10/19 17:56:50  mdw
 * Fix bit operations.  Test them (a bit) better.
 *
 * Revision 1.14  2002/10/09 00:36:03  mdw
 * Fix bounds on workspace for Karatsuba operations.
 *
 * Revision 1.13  2002/10/06 22:52:50  mdw
 * Pile of changes for supporting two's complement properly.
 *
 * Revision 1.12  2001/04/03 19:36:05  mdw
 * Add some simple bitwise operations so that Perl can use them.
 *
 * Revision 1.11  2000/10/08 15:48:35  mdw
 * Rename Karatsuba constants now that we have @gfx_kmul@ too.
 *
 * Revision 1.10  2000/10/08 12:06:12  mdw
 * Provide @mpx_ueq@ for rapidly testing equality of two integers.
 *
 * Revision 1.9  1999/12/22 15:49:07  mdw
 * New function for division by a small integer.
 *
 * Revision 1.8  1999/12/11 10:57:43  mdw
 * Karatsuba squaring algorithm.
 *
 * Revision 1.7  1999/12/11 01:51:28  mdw
 * Change Karatsuba parameters slightly.
 *
 * Revision 1.6  1999/12/10 23:23:51  mdw
 * Karatsuba-Ofman multiplication algorithm.
 *
 * 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:43  mdw
 * Add two's complement support.  Fix a bug in MPX_UMLAN.
 *
 * Revision 1.3  1999/11/13 01:51:29  mdw
 * Minor interface changes.  Should be stable now.
 *
 * Revision 1.2  1999/11/11 17:47:55  mdw
 * Minor changes for different `mptypes.h' format.
 *
 * Revision 1.1  1999/09/03 08:41:12  mdw
 * Initial import.
 *
 */

#ifndef CATACOMB_MPX_H
#define CATACOMB_MPX_H

#ifdef __cplusplus
  extern "C" {
#endif

/*----- The idea ----------------------------------------------------------*
 *
 * This file provides functions and macros which work on vectors of words as
 * unsigned multiprecision integers.  The interface works in terms of base
 * and limit pointers (i.e., a pointer to the start of a vector, and a
 * pointer just past its end) rather than base pointer and length, because
 * that requires more arithmetic and state to work on.
 *
 * The interfaces are slightly bizarre in other ways.  Try to use the
 * higher-level functions where you can: they're rather better designed to
 * actually be friendly and useful.
 */

/*----- Header files ------------------------------------------------------*/

#include <string.h>

#ifndef CATACOMB_MPW_H
#  include "mpw.h"
#endif

/*----- General manipulation ----------------------------------------------*/

/* --- @MPX_SHRINK@ --- *
 *
 * Arguments:   @const mpw *v@ = pointer to vector of words
 *              @const mpw *vl@ = (updated) current limit of vector
 *
 * Use:         Shrinks down the limit of a multiprecision integer vector.
 */

#define MPX_SHRINK(v, vl) do {                                          \
  const mpw *_vv = (v), *_vvl = (vl);                                   \
  while (_vvl > _vv && !_vvl[-1])                                       \
    _vvl--;                                                             \
  (vl) = (mpw *)_vvl;                                                   \
} while (0)

/* --- @MPX_BITS@ --- *
 *
 * Arguments:   @unsigned long b@ = result variable
 *              @const mpw *v@ = pointer to array of words
 *              @const mpw *vl@ = limit of vector (from @MPX_SHRINK@)
 *
 * Use:         Calculates the number of bits in a multiprecision value.
 */

#define MPX_BITS(b, v, vl) do {                                         \
  const mpw *_v = (v), *_vl = (vl);                                     \
  MPX_SHRINK(_v, _vl);                                                  \
  if (_v == _vl)                                                        \
    (b) = 0;                                                            \
  else {                                                                \
    unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1;                   \
    mpw _w = _vl[-1];                                                   \
    unsigned _k = MPW_BITS / 2;                                         \
    while (_k) {                                                        \
      if (_w >> _k) {                                                   \
        _w >>= _k;                                                      \
        _b += _k;                                                       \
      }                                                                 \
      _k >>= 1;                                                         \
    }                                                                   \
    (b) = _b;                                                           \
  }                                                                     \
} while (0)

/* --- @MPX_OCTETS@ --- *
 *
 * Arguments:   @size_t o@ = result variable
 *              @const mpw *v, *vl@ = pointer to array of words
 *
 * Use:         Calculates the number of octets in a multiprecision value.
 */

#define MPX_OCTETS(o, v, vl) do {                                       \
  unsigned long _bb;                                                    \
  MPX_BITS(_bb, (v), (vl));                                             \
  (o) = (_bb + 7) >> 3;                                                 \
} while (0)

/* --- @MPX_OCTETS2C@ --- *
 *
 * Arguments:   @size_t o@ = result variable
 *              @const mpw *v, *vl@ = pointer to array of words
 *
 * Use:         Calculates the number of octets in a multiprecision value, if
 *              you represent it as two's complement.
 */

#define MPX_OCTETS2C(o, v, vl) do {                                     \
  unsigned long _bb;                                                    \
  MPX_BITS(_bb, (v), (vl));                                             \
  (o) = (_bb >> 3) + 1;                                                 \
} while (0)

/* --- @MPX_COPY@ --- *
 *
 * Arguments:   @dv, dvl@ = destination vector base and limit
 *              @av, avl@ = source vector base and limit
 *
 * Use:         Copies a multiprecision integer.
 */

#define MPX_COPY(dv, dvl, av, avl) do {                                 \
  mpw *_dv = (dv), *_dvl = (dvl);                                       \
  size_t _dn = _dvl - _dv;                                              \
  const mpw *_av = (av), *_avl = (avl);                                 \
  size_t _an = _avl - _av;                                              \
  if (_av == _dv) {                                                     \
    if (_dvl > _avl)                                                    \
      memset(_dv, 0, MPWS(_dn - _an));                                  \
  } else if (_an >= _dn)                                                \
    memmove(_dv, _av, MPWS(_dn));                                       \
  else {                                                                \
    memmove(_dv, _av, MPWS(_an));                                       \
    memset(_dv + _an, 0, MPWS(_dn - _an));                              \
  }                                                                     \
} while (0)

/* --- @MPX_ZERO@ --- *
 *
 * Arguments:   @v, vl@ = base and limit of vector to clear
 *
 * Use:         Zeroes the area between the two vector pointers.
 */

#define MPX_ZERO(v, vl) do {                                            \
  mpw *_v = (v), *_vl = (vl);                                           \
  if (_v < _vl)                                                         \
    memset(_v, 0, MPWS(_vl - _v));                                      \
} while (0)

/*----- Loading and storing -----------------------------------------------*/

/* --- @mpx_storel@ --- *
 *
 * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 *              @void *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.
 */

extern void mpx_storel(const mpw */*v*/, const mpw */*vl*/,
                       void */*p*/, size_t /*sz*/);

/* --- @mpx_loadl@ --- *
 *
 * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 *              @const void *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.  High-end octets are ignored if there isn't enough
 *              space for them.
 */

extern void mpx_loadl(mpw */*v*/, mpw */*vl*/,
                      const void */*p*/, size_t /*sz*/);

/* --- @mpx_storeb@ --- *
 *
 * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 *              @void *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.
 */

extern void mpx_storeb(const mpw */*v*/, const mpw */*vl*/,
                       void */*p*/, size_t /*sz*/);

/* --- @mpx_loadb@ --- *
 *
 * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 *              @const void *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.  High-end octets are ignored if there isn't enough
 *              space for them.
 */

extern void mpx_loadb(mpw */*v*/, mpw */*vl*/,
                      const void */*p*/, size_t /*sz*/);

/* --- @mpx_storel2cn@ --- *
 *
 * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 *              @void *pp@ = pointer to octet array
 *              @size_t sz@ = size of octet array
 *
 * Returns:     ---
 *
 * Use:         Stores a negative MP in an octet array, least significant
 *              octet first, as two's complement.  High-end octets are
 *              silently discarded if there isn't enough space for them.
 *              This obviously makes the output bad.
 */

extern void mpx_storel2cn(const mpw */*v*/, const mpw */*vl*/,
                          void */*p*/, size_t /*sz*/);

/* --- @mpx_loadl2cn@ --- *
 *
 * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 *              @const void *pp@ = pointer to octet array
 *              @size_t sz@ = size of octet array
 *
 * Returns:     ---
 *
 * Use:         Loads a negative MP in an octet array, least significant
 *              octet first, as two's complement.  High-end octets are
 *              ignored if there isn't enough space for them.  This probably
 *              means you made the wrong choice coming here.
 */

extern void mpx_loadl2cn(mpw */*v*/, mpw */*vl*/,
                         const void */*p*/, size_t /*sz*/);

/* --- @mpx_storeb2cn@ --- *
 *
 * Arguments:   @const mpw *v, *vl@ = base and limit of source vector
 *              @void *pp@ = pointer to octet array
 *              @size_t sz@ = size of octet array
 *
 * Returns:     ---
 *
 * Use:         Stores a negative MP in an octet array, most significant
 *              octet first, as two's complement.  High-end octets are
 *              silently discarded if there isn't enough space for them,
 *              which probably isn't what you meant.
 */

extern void mpx_storeb2cn(const mpw */*v*/, const mpw */*vl*/,
                          void */*p*/, size_t /*sz*/);

/* --- @mpx_loadb2cn@ --- *
 *
 * Arguments:   @mpw *v, *vl@ = base and limit of destination vector
 *              @const void *pp@ = pointer to octet array
 *              @size_t sz@ = size of octet array
 *
 * Returns:     ---
 *
 * Use:         Loads a negative MP in an octet array, most significant octet
 *              first as two's complement.  High-end octets are ignored if
 *              there isn't enough space for them.  This probably means you
 *              chose this function wrongly.
 */

extern void mpx_loadb2cn(mpw */*v*/, mpw */*vl*/,
                         const void */*p*/, size_t /*sz*/);


/*----- Logical shifting --------------------------------------------------*/

/* --- @mpx_lsl@ --- *
 *
 * 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.
 */

extern void mpx_lsl(mpw */*dv*/, mpw */*dvl*/,
                    const mpw */*av*/, const mpw */*avl*/,
                    size_t /*n*/);

/* --- @mpx_lsr@ --- *
 *
 * 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 right operation on an integer.
 */

extern void mpx_lsr(mpw */*dv*/, mpw */*dvl*/,
                    const mpw */*av*/, const mpw */*avl*/,
                    size_t /*n*/);

/*----- Bitwise operations ------------------------------------------------*/

/* --- @mpx_bitop@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = destination vector
 *              @const mpw *av, *avl@ = first source vector
 *              @const mpw *bv, *bvl@ = second source vector
 *
 * Returns:     ---
 *
 * Use:         Provide the dyadic boolean functions.  The functions are
 *              named after the truth table they generate:
 *
 *                      a:      0011
 *                      b:      0101
 *                      @mpx_bitXXXX@
 */

#define MPX_DOBIN(what)                                                 \
  what(0000) what(0001) what(0010) what(0011)                           \
  what(0100) what(0101) what(0110) what(0111)                           \
  what(1000) what(1001) what(1010) what(1011)                           \
  what(1100) what(1101) what(1110) what(1111)

#define MPX_BITDECL(string)                                             \
  extern void mpx_bit##string(mpw */*dv*/, mpw */*dvl*/,                \
                              const mpw */*av*/, const mpw */*avl*/,    \
                              const mpw */*bv*/, const mpw */*bvl*/);
MPX_DOBIN(MPX_BITDECL)

/* --- @mpx_[n]and@, @mpx_[n]or@, @mpx_xor@ --- *
 *
 * Synonyms for the commonly-used functions above.
 */

#define mpx_and  mpx_bit0001
#define mpx_or   mpx_bit0111
#define mpx_nand mpx_bit1110
#define mpx_nor  mpx_bit1000
#define mpx_xor  mpx_bit0110

/* --- @mpx_not@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = destination vector
 *              @const mpw *av, *avl@ = first source vector
 *
 * Returns:     ---
 *
 * Use;         Bitwise NOT.
 */

extern void mpx_not(mpw */*dv*/, mpw */*dvl*/,
                    const mpw */*av*/, const mpw */*avl*/);

/*----- Unsigned arithmetic -----------------------------------------------*/

/* --- @mpx_2c@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = destination vector
 *              @const mpw *v, *vl@ = source vector
 *
 * Returns:     ---
 *
 * Use:         Calculates the two's complement of @v@.
 */

extern void mpx_2c(mpw */*dv*/, mpw */*dvl*/,
                   const mpw */*v*/, const mpw */*vl*/);

/* --- @mpx_ueq@ --- *
 *
 * Arguments:   @const mpw *av, *avl@ = first argument vector base and limit
 *              @const mpw *bv, *bvl@ = second argument vector base and limit
 *
 * Returns:     Nonzero if the two vectors are equal.
 *
 * Use:         Performs an unsigned integer test for equality.
 */

extern int mpx_ueq(const mpw */*av*/, const mpw */*avl*/,
                   const mpw */*bv*/, const mpw */*bvl*/);

/* --- @mpx_ucmp@ --- *
 *
 * Arguments:   @const mpw *av, *avl@ = first argument vector base and limit
 *              @const mpw *bv, *bvl@ = second argument vector base and limit
 *
 * Returns:     Less than, equal to, or greater than zero depending on
 *              whether @a@ is less than, equal to or greater than @b@,
 *              respectively.
 *
 * Use:         Performs an unsigned integer comparison.
 */

#define MPX_UCMP(av, avl, op, dv, dvl)                                  \
  (mpx_ucmp((av), (avl), (dv), (dvl)) op 0)

extern int mpx_ucmp(const mpw */*av*/, const mpw */*avl*/,
                    const mpw */*bv*/, const mpw */*bvl*/);

/* --- @mpx_uadd@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 *              @const mpw *av, *avl@ = first addend vector base and limit
 *              @const mpw *bv, *bvl@ = second addend vector base and limit
 *
 * Returns:     ---
 *
 * Use:         Performs unsigned integer addition.  If the result overflows
 *              the destination vector, high-order bits are discarded.  This
 *              means that two's complement addition happens more or less for
 *              free, although that's more a side-effect than anything else.
 *              The result vector may be equal to either or both source
 *              vectors, but may not otherwise overlap them.
 */

extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/,
                     const mpw */*av*/, const mpw */*avl*/,
                     const mpw */*bv*/, const mpw */*bvl*/);

/* --- @mpx_uaddn@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = source and destination base and limit
 *              @mpw n@ = other addend
 *
 * Returns:     ---
 *
 * Use:         Adds a small integer to a multiprecision number.
 */

#define MPX_UADDN(dv, dvl, n) do {                                      \
  mpw *_ddv = (dv), *_ddvl = (dvl);                                     \
  mpw _c = (n);                                                         \
                                                                        \
  while (_c && _ddv < _ddvl) {                                          \
    mpd _x = (mpd)*_ddv + (mpd)_c;                                      \
    *_ddv++ = MPW(_x);                                                  \
    _c = _x >> MPW_BITS;                                                \
  }                                                                     \
} while (0)

extern void mpx_uaddn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);

/* --- @mpx_usub@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 *              @const mpw *av, *avl@ = first argument vector base and limit
 *              @const mpw *bv, *bvl@ = second argument vector base and limit
 *
 * Returns:     ---
 *
 * Use:         Performs unsigned integer subtraction.  If the result
 *              overflows the destination vector, high-order bits are
 *              discarded.  This means that two's complement subtraction
 *              happens more or less for free, although that's more a side-
 *              effect than anything else.  The result vector may be equal to
 *              either or both source vectors, but may not otherwise overlap
 *              them.
 */

extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
                     const mpw */*av*/, const mpw */*avl*/,
                     const mpw */*bv*/, const mpw */*bvl*/);

/* --- @mpx_usubn@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = source and destination base and limit
 *              @n@ = subtrahend
 *
 * Returns:     ---
 *
 * Use:         Subtracts a small integer from a multiprecision number.
 */

#define MPX_USUBN(dv, dvl, n) do {                                      \
  mpw *_ddv = (dv), *_ddvl = (dvl);                                     \
  mpw _c = (n);                                                         \
                                                                        \
  while (_ddv < _ddvl) {                                                \
    mpd _x = (mpd)*_ddv - (mpd)_c;                                      \
    *_ddv++ = MPW(_x);                                                  \
    if (_x >> MPW_BITS)                                                 \
      _c = 1;                                                           \
    else                                                                \
      break;                                                            \
  }                                                                     \
} while (0)

extern void mpx_usubn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);

/* --- @mpx_umul@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 *              @const mpw *av, *avl@ = multiplicand vector base and limit
 *              @const mpw *bv, *bvl@ = multiplier vector base and limit
 *
 * Returns:     ---
 *
 * Use:         Performs unsigned integer multiplication.  If the result
 *              overflows the desination vector, high-order bits are
 *              discarded.  The result vector may not overlap the argument
 *              vectors in any way.
 */

extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/,
                     const mpw */*av*/, const mpw */*avl*/,
                     const mpw */*bv*/, const mpw */*bvl*/);

/* --- @mpx_umuln@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 *              @const mpw *av, *avl@ = multiplicand vector base and limit
 *              @mpw m@ = multiplier
 *
 * Returns:     ---
 *
 * Use:         Multiplies a multiprecision integer by a single-word value.
 *              The destination and source may be equal.  The destination
 *              is completely cleared after use.
 */

#define MPX_UMULN(dv, dvl, av, avl, m) do {                             \
  mpw *_dv = (dv), *_dvl = (dvl);                                       \
  const mpw *_av = (av), *_avl = (avl);                                 \
  mpw _c = 0;                                                           \
  mpd _m = (m);                                                         \
                                                                        \
  while (_av < _avl) {                                                  \
    mpd _x;                                                             \
    if (_dv >= _dvl)                                                    \
      break;                                                            \
    _x = (mpd)_m * (mpd)*_av++ + _c;                                    \
    *_dv++ = MPW(_x);                                                   \
    _c = _x >> MPW_BITS;                                                \
  }                                                                     \
  if (_dv < _dvl) {                                                     \
    *_dv++ = MPW(_c);                                                   \
    MPX_ZERO(_dv, _dvl);                                                \
  }                                                                     \
} while (0)

extern void mpx_umuln(mpw */*dv*/, mpw */*dvl*/,
                      const mpw */*av*/, const mpw */*avl*/, mpw m);

/* --- @mpx_umlan@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = destination/accumulator base and limit
 *              @const mpw *av, *avl@ = multiplicand vector base and limit
 *              @mpw m@ = multiplier
 *
 * Returns:     ---
 *
 * Use:         Multiplies a multiprecision integer by a single-word value
 *              and adds the result to an accumulator.
 */

#define MPX_UMLAN(dv, dvl, av, avl, m) do {                             \
  mpw *_dv = (dv), *_dvl = (dvl);                                       \
  const mpw *_av = (av), *_avl = (avl);                                 \
  mpw _cc = 0;                                                          \
  mpd _m = (m);                                                         \
                                                                        \
  while (_dv < _dvl && _av < _avl) {                                    \
    mpd _x;                                                             \
    _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc;                       \
    *_dv++ = MPW(_x);                                                   \
    _cc = _x >> MPW_BITS;                                               \
  }                                                                     \
  MPX_UADDN(_dv, _dvl, _cc);                                            \
} while (0)

extern void mpx_umlan(mpw */*dv*/, mpw */*dvl*/,
                      const mpw */*av*/, const mpw */*avl*/, mpw m);

/* --- @mpx_usqr@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 *              @const mpw *av, *av@ = source vector base and limit
 *
 * Returns:     ---
 *
 * Use:         Performs unsigned integer squaring.  The result vector must
 *              not overlap the source vector in any way.
 */

extern void mpx_usqr(mpw */*dv*/, mpw */*dvl*/,
                     const mpw */*av*/, const mpw */*avl*/);

/* --- @mpx_udiv@ --- *
 *
 * Arguments:   @mpw *qv, *qvl@ = quotient vector base and limit
 *              @mpw *rv, *rvl@ = dividend/remainder vector base and limit
 *              @const mpw *dv, *dvl@ = divisor vector base and limit
 *              @mpw *sv, *svl@ = scratch workspace 
 *
 * Returns:     ---
 *
 * Use:         Performs unsigned integer division.  If the result overflows
 *              the quotient vector, high-order bits are discarded.  (Clearly
 *              the remainder vector can't overflow.)  The various vectors
 *              may not overlap in any way.  Yes, I know it's a bit odd
 *              requiring the dividend to be in the result position but it
 *              does make some sense really.  The remainder must have
 *              headroom for at least two extra words.  The scratch space
 *              must be at least one word larger than the divisor.
 */

extern void mpx_udiv(mpw */*qv*/, mpw */*qvl*/, mpw */*rv*/, mpw */*rvl*/,
                     const mpw */*dv*/, const mpw */*dvl*/,
                     mpw */*sv*/, mpw */*svl*/);

/* --- @mpx_udivn@ --- *
 *
 * Arguments:   @mpw *qv, *qvl@ = storage for the quotient (may overlap
 *                      dividend)
 *              @const mpw *rv, *rvl@ = dividend
 *              @mpw d@ = single-precision divisor
 *
 * Returns:     Remainder after divison.
 *
 * Use:         Performs a single-precision division operation.
 */

extern mpw mpx_udivn(mpw */*qv*/, mpw */*qvl*/,
                     const mpw */*rv*/, const mpw */*rvl*/, mpw /*d*/);

/*----- Karatsuba multiplication algorithms -------------------------------*/

/* --- @MPK_THRESH@ --- *
 *
 * This is the limiting length for using Karatsuba algorithms.  It's best to
 * use the simpler classical multiplication method on numbers smaller than
 * this.  It is unsafe to make this constant less than four (i.e., the
 * algorithms will fail).
 */

#define MPK_THRESH 16

/* --- @mpx_kmul@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = pointer to destination buffer
 *              @const mpw *av, *avl@ = pointer to first argument
 *              @const mpw *bv, *bvl@ = pointer to second argument
 *              @mpw *sv, *svl@ = pointer to scratch workspace
 *
 * Returns:     ---
 *
 * Use:         Multiplies two multiprecision integers using Karatsuba's
 *              algorithm.  This is rather faster than traditional long
 *              multiplication (e.g., @mpx_umul@) on large numbers, although
 *              more expensive on small ones.
 *
 *              The destination must be three times as large as the larger
 *              argument.  The scratch space must be five times as large as
 *              the larger argument.
 */

extern void mpx_kmul(mpw */*dv*/, mpw */*dvl*/,
                     const mpw */*av*/, const mpw */*avl*/,
                     const mpw */*bv*/, const mpw */*bvl*/,
                     mpw */*sv*/, mpw */*svl*/);

/* --- @mpx_ksqr@ --- *
 *
 * Arguments:   @mpw *dv, *dvl@ = pointer to destination buffer
 *              @const mpw *av, *avl@ = pointer to first argument
 *              @mpw *sv, *svl@ = pointer to scratch workspace
 *
 * Returns:     ---
 *
 * Use:         Squares a multiprecision integers using something similar to
 *              Karatsuba's multiplication algorithm.  This is rather faster
 *              than traditional long multiplication (e.g., @mpx_umul@) on
 *              large numbers, although more expensive on small ones, and
 *              rather simpler than full-blown Karatsuba multiplication.
 *
 *              The destination must be three times as large as the larger
 *              argument.  The scratch space must be five times as large as
 *              the larger argument.
 */

extern void mpx_ksqr(mpw */*dv*/, mpw */*dvl*/,
                     const mpw */*av*/, const mpw */*avl*/,
                     mpw */*sv*/, mpw */*svl*/);

/*----- That's all, folks -------------------------------------------------*/

#ifdef __cplusplus
  }
#endif

#endif