X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/ba6e6b64033b1f9de49feccb5c9cd438354481f7..0f00dc4c8eb47e67bc0f148c2dd109f73a451e0a:/math/mpx.h

diff --git a/math/mpx.h b/math/mpx.h
new file mode 100644
index 0000000..4d6288e
--- /dev/null
+++ b/math/mpx.h
@@ -0,0 +1,789 @@
+/* -*-c-*-
+ *
+ * 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.
+ */
+
+#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_P2;						\
+    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)
+
+/* --- @MPX_ONE@ --- *
+ *
+ * Arguments:	@v, vl@ = base and limit of vector to clear
+ *
+ * Use:		Fills the area between the two vector pointers with ones.
+ */
+
+#define MPX_ONE(v, vl) do {						\
+  mpw * _v = (v);							\
+  const mpw *_vl = (vl);						\
+  while (_v < _vl)							\
+    *_v++ = MPW_MAX;							\
+} 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_lslc@ --- *
+ *
+ * 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, only
+ *		it fills in the bits with ones instead of zeroes.
+ */
+
+extern void mpx_lslc(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_uaddnlsl@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination and first argument vector
+ *		@mpw a@ = second argument
+ *		@unsigned o@ = offset in bits
+ *
+ * Returns:	---
+ *
+ * Use:		Computes %$d + 2^o a$%.  If the result overflows then
+ *		high-order bits are discarded, as usual.  We must have
+ *		@0 < o < MPW_BITS@.
+ */
+
+extern void mpx_uaddnlsl(mpw */*dv*/, mpw */*dvl*/,
+			 mpw /*a*/, unsigned /*o*/);
+
+/* --- @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_usubnlsl@ --- *
+ *
+ * Arguments:	@mpw *dv, *dvl@ = destination and first argument vector
+ *		@mpw a@ = second argument
+ *		@unsigned o@ = offset in bits
+ *
+ * Returns:	---
+ *
+ * Use:		Computes %$d - 2^o a$%.  If the result overflows then
+ *		high-order bits are discarded, as usual, so you get two's
+ *		complement.  Which might be what you wanted...  We must have
+ *		@0 < o < MPW_BITS@.
+ */
+
+extern void mpx_usubnlsl(mpw */*dv*/, mpw */*dvl*/,
+			 mpw /*a*/, unsigned /*o*/);
+
+/* --- @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