Merge branch 'master' of git.distorted.org.uk:~mdw/publish/public-git/catacomb
[u/mdw/catacomb] / mpx.h
diff --git a/mpx.h b/mpx.h
index 0f1a91d..19f5cc7 100644 (file)
--- a/mpx.h
+++ b/mpx.h
@@ -1,13 +1,13 @@
 /* -*-c-*-
  *
- * $Id: mpx.h,v 1.2 1999/11/11 17:47:55 mdw Exp $
+ * $Id: mpx.h,v 1.18 2004/04/08 01:36:15 mdw Exp $
  *
  * Low level multiprecision arithmetic
  *
  * (c) 1999 Straylight/Edgeware
  */
 
-/*----- Licensing notice --------------------------------------------------* 
+/*----- Licensing notice --------------------------------------------------*
  *
  * This file is part of Catacomb.
  *
  * it under the terms of the GNU Library General Public License as
  * published by the Free Software Foundation; either version 2 of the
  * License, or (at your option) any later version.
- * 
+ *
  * Catacomb is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU Library General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU Library General Public
  * License along with Catacomb; if not, write to the Free
  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  * MA 02111-1307, USA.
  */
 
-/*----- Revision history --------------------------------------------------* 
- *
- * $Log: mpx.h,v $
- * 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 MPX_H
-#define MPX_H
+#ifndef CATACOMB_MPX_H
+#define CATACOMB_MPX_H
 
 #ifdef __cplusplus
   extern "C" {
 
 #include <string.h>
 
-#ifndef MPTYPES_H
-#  include "mptypes.h"
+#ifndef CATACOMB_MPW_H
+#  include "mpw.h"
 #endif
 
-/*----- Basic macros ------------------------------------------------------*/
-
-/* --- @MPW@ --- *
- *
- * Arguments:   @x@ = an unsigned value
- *
- * Use:         Expands to the value of @x@ masked and typecast to a
- *              multiprecision integer word.
- */
-
-#define MPW(x) ((mpw)((x) & MPW_MAX))
-
-/* --- @MPWS@ --- *
- *
- * Arguments:   @n@ = number of words
- *
- * Use:         Expands to the number of bytes occupied by a given number of
- *              words.
- */
-
-#define MPWS(n) ((n) * sizeof(mpw))
-
 /*----- General manipulation ----------------------------------------------*/
 
 /* --- @MPX_SHRINK@ --- *
  */
 
 #define MPX_SHRINK(v, vl) do {                                         \
-  const mpw *_v = (v), *_vl = (vl);                                    \
-  while (_vl > _v && *--_vl)                                           \
-    ;                                                                  \
-  (vl) = _vl;                                                          \
+  const mpw *_vv = (v), *_vvl = (vl);                                  \
+  while (_vvl > _vv && !_vvl[-1])                                      \
+    _vvl--;                                                            \
+  (vl) = (mpw *)_vvl;                                                  \
 } while (0)
 
 /* --- @MPX_BITS@ --- *
 
 #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;                  \
+    unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1;                  \
     mpw _w = _vl[-1];                                                  \
-    unsigned _k = MPW_BITS / 2;                                                \
+    unsigned _k = MPW_P2;                                              \
     while (_k) {                                                       \
       if (_w >> _k) {                                                  \
        _w >>= _k;                                                      \
 /* --- @MPX_OCTETS@ --- *
  *
  * Arguments:  @size_t o@ = result variable
- *             @const mpw *v@ = pointer to array of words
- *             @size_t len@ = length of the array (from @MPX_LEN@)
+ *             @const mpw *v, *vl@ = pointer to array of words
  *
  * Use:                Calculates the number of octets in a multiprecision value.
  */
 
-#define MPX_OCTETS(o, v, len) do {                                     \
-  const mpw *_v = (v), *_vl = (vl);                                    \
-  if (_v == _vl)                                                       \
-    (o) = 0;                                                           \
-  else {                                                               \
-    _size_t _o = (MPW_BITS / 8) * (_vl - _v - 1);                      \
-    mpw _w = _vl[-1];                                                  \
-    unsigned _k = MPW_BITS / 2;                                                \
-    while (_k > 3) {                                                   \
-      if (_w >> _k) {                                                  \
-       _w >>= _k;                                                      \
-       _o += _k - 3;                                                   \
-      }                                                                        \
-      _k >>= 1;                                                                \
-    }                                                                  \
-    if (_w)                                                            \
-      _o++;                                                            \
-    (o) = _o;                                                          \
-  }                                                                    \
+#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@ --- *
  * Use:                Copies a multiprecision integer.
  */
 
-#define MPX_COPY(dv, dvl, av, dvl) do {                                        \
-  mpw *_dv = (dv);                                                     \
-  size_t _dn = (dvl) - _dv;                                            \
-  const mpw *_av = (av);                                               \
-  size_t _an = (avl) - _av;                                            \
+#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(_avl, 0, MPWS(_dn - _an));                                        \
+      memset(_dv, 0, MPWS(_dn - _an));                                 \
   } else if (_an >= _dn)                                               \
     memmove(_dv, _av, MPWS(_dn));                                      \
   else {                                                               \
  * Use:                Zeroes the area between the two vector pointers.
  */
 
-#define MPX_ZERO(v, vl) {                                              \
+#define MPX_ZERO(v, vl) do {                                           \
   mpw *_v = (v), *_vl = (vl);                                          \
-  memset(_v, 0, MPWS(_vl - _v));                                       \
+  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
- *             @octet *p@ = pointer to octet array
+ *             @void *p@ = pointer to octet array
  *             @size_t sz@ = size of octet array
  *
  * Returns:    ---
  */
 
 extern void mpx_storel(const mpw */*v*/, const mpw */*vl*/,
-                      octet */*p*/, size_t /*sz*/);
+                      void */*p*/, size_t /*sz*/);
 
 /* --- @mpx_loadl@ --- *
  *
  * Arguments:  @mpw *v, *vl@ = base and limit of destination vector
- *             @const octet *p@ = pointer to octet array
+ *             @const void *p@ = pointer to octet array
  *             @size_t sz@ = size of octet array
  *
  * Returns:    ---
@@ -231,12 +213,12 @@ extern void mpx_storel(const mpw */*v*/, const mpw */*vl*/,
  */
 
 extern void mpx_loadl(mpw */*v*/, mpw */*vl*/,
-                     const octet */*p*/, size_t /*sz*/);
+                     const void */*p*/, size_t /*sz*/);
 
 /* --- @mpx_storeb@ --- *
  *
  * Arguments:  @const mpw *v, *vl@ = base and limit of source vector
- *             @octet *p@ = pointer to octet array
+ *             @void *p@ = pointer to octet array
  *             @size_t sz@ = size of octet array
  *
  * Returns:    ---
@@ -247,12 +229,12 @@ extern void mpx_loadl(mpw */*v*/, mpw */*vl*/,
  */
 
 extern void mpx_storeb(const mpw */*v*/, const mpw */*vl*/,
-                      octet */*p*/, size_t /*sz*/);
+                      void */*p*/, size_t /*sz*/);
 
 /* --- @mpx_loadb@ --- *
  *
  * Arguments:  @mpw *v, *vl@ = base and limit of destination vector
- *             @const octet *p@ = pointer to octet array
+ *             @const void *p@ = pointer to octet array
  *             @size_t sz@ = size of octet array
  *
  * Returns:    ---
@@ -263,7 +245,76 @@ extern void mpx_storeb(const mpw */*v*/, const mpw */*vl*/,
  */
 
 extern void mpx_loadb(mpw */*v*/, mpw */*vl*/,
-                     const octet */*p*/, size_t /*sz*/);
+                     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 --------------------------------------------------*/
 
@@ -282,6 +333,22 @@ 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
@@ -297,8 +364,88 @@ 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
@@ -337,6 +484,45 @@ 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
@@ -348,7 +534,7 @@ extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/,
  * 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, althuogh that's more a side-
+ *             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.
@@ -358,11 +544,74 @@ extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
                     const mpw */*av*/, const mpw */*avl*/,
                     const mpw */*bv*/, const mpw */*bvl*/);
 
-/* --- @MPX_UMULN@ --- *
+/* --- @mpx_usubn@ --- *
  *
- * Arguments:  @dv, dvl@ = destination vector base and limit
- *             @av, avl@ = multiplicand vector base and limit
- *             @m@ = multiplier
+ * 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
@@ -379,7 +628,7 @@ extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
     mpd _x;                                                            \
     if (_dv >= _dvl)                                                   \
       break;                                                           \
-    _x = _m * *_av++ + c;                                              \
+    _x = (mpd)_m * (mpd)*_av++ + _c;                                   \
     *_dv++ = MPW(_x);                                                  \
     _c = _x >> MPW_BITS;                                               \
   }                                                                    \
@@ -389,11 +638,16 @@ extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
   }                                                                    \
 } while (0)
 
-/* --- @MPX_UMLAN@ --- *
+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
  *
- * Arguments:  @dv, dvl@ = destination/accumulator vector base and limit
- *             @av, avl@ = multiplicand vector base and limit
- *             @m@ = multiplier
+ * Returns:    ---
  *
  * Use:                Multiplies a multiprecision integer by a single-word value
  *             and adds the result to an accumulator.
@@ -402,46 +656,41 @@ extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
 #define MPX_UMLAN(dv, dvl, av, avl, m) do {                            \
   mpw *_dv = (dv), *_dvl = (dvl);                                      \
   const mpw *_av = (av), *_avl = (avl);                                        \
-  mpw _c = 0;                                                          \
+  mpw _cc = 0;                                                         \
   mpd _m = (m);                                                                \
                                                                        \
-  while (_av < _avl) {                                                 \
+  while (_dv < _dvl && _av < _avl) {                                   \
     mpd _x;                                                            \
-    if (_dv >= _dvl)                                                   \
-      break;                                                           \
-    _x = *_dv + _m * *_av++ + _c;                                      \
+    _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc;                      \
     *_dv++ = MPW(_x);                                                  \
-    _c = _x >> MPW_BITS;                                               \
-  }                                                                    \
-  if (_dv < _dvl) {                                                    \
-    *_dv++ = MPW(_c);                                                  \
-    MPX_ZERO(_dv, _dvl);                                               \
+    _cc = _x >> MPW_BITS;                                              \
   }                                                                    \
+  MPX_UADDN(_dv, _dvl, _cc);                                           \
 } while (0)
 
-/* --- @mpx_umul@ --- *
+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, *avl@ = multiplicand vector base and limit
- *             @const mpw *bv, *bvl@ = multiplier vector base and limit
+ *             @const mpw *av, *av@ = source 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.
+ * Use:                Performs unsigned integer squaring.  The result vector must
+ *             not overlap the source vector in any way.
  */
 
-extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/,
-                    const mpw */*av*/, const mpw */*avl*/,
-                    const mpw */*bv*/, const mpw */*bvl*/);
+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:    ---
  *
@@ -450,11 +699,88 @@ extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/,
  *             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.
+ *             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*/);
+                    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 -------------------------------------------------*/