More changes. Still embryonic.
[u/mdw/catacomb] / mpx.h
diff --git a/mpx.h b/mpx.h
index 9429e91..509da02 100644 (file)
--- a/mpx.h
+++ b/mpx.h
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mpx.h,v 1.3 1999/11/13 01:51:29 mdw Exp $
+ * $Id: mpx.h,v 1.8 1999/12/11 10:57:43 mdw Exp $
  *
  * Low level multiprecision arithmetic
  *
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mpx.h,v $
+ * 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.
  *
@@ -41,8 +56,8 @@
  *
  */
 
-#ifndef MPX_H
-#define MPX_H
+#ifndef CATACOMB_MPX_H
+#define CATACOMB_MPX_H
 
 #ifdef __cplusplus
   extern "C" {
@@ -65,7 +80,7 @@
 
 #include <string.h>
 
-#ifndef MPW_H
+#ifndef CATACOMB_MPW_H
 #  include "mpw.h"
 #endif
 
@@ -282,6 +297,19 @@ extern void mpx_lsr(mpw */*dv*/, mpw */*dvl*/,
 
 /*----- 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_ucmp@ --- *
  *
  * Arguments:  @const mpw *av, *avl@ = first argument vector base and limit
@@ -320,10 +348,12 @@ extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/,
                     const mpw */*av*/, const mpw */*avl*/,
                     const mpw */*bv*/, const mpw */*bvl*/);
 
-/* --- @MPX_UADDN@ --- *
+/* --- @mpx_uaddn@ --- *
  *
- * Arguments:  @dv, dvl@ = source and destination vector base and limit
- *             @n@ = other addend
+ * Arguments:  @mpw *dv, *dvl@ = source and destination base and limit
+ *             @mpw n@ = other addend
+ *
+ * Returns:    ---
  *
  * Use:                Adds a small integer to a multiprecision number.
  */
@@ -339,6 +369,8 @@ extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/,
   }                                                                    \
 } while (0)
 
+extern void mpx_uaddn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);
+
 /* --- @mpx_usub@ --- *
  *
  * Arguments:  @mpw *dv, *dvl@ = destination vector base and limit
@@ -360,10 +392,12 @@ extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
                     const mpw */*av*/, const mpw */*avl*/,
                     const mpw */*bv*/, const mpw */*bvl*/);
 
-/* --- @MPX_USUBN@ --- *
+/* --- @mpx_usubn@ --- *
  *
- * Arguments:  @@dv, dvl@ = destination vector base and limit
- *             @n@ = other addend
+ * Arguments:  @mpw *dv, *dvl@ = source and destination base and limit
+ *             @n@ = subtrahend
+ *
+ * Returns:    ---
  *
  * Use:                Subtracts a small integer from a multiprecision number.
  */
@@ -382,6 +416,8 @@ extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
   }                                                                    \
 } while (0)
 
+extern void mpx_usubn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);
+
 /* --- @mpx_umul@ --- *
  *
  * Arguments:  @mpw *dv, *dvl@ = destination vector base and limit
@@ -400,11 +436,13 @@ extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/,
                     const mpw */*av*/, const mpw */*avl*/,
                     const mpw */*bv*/, const mpw */*bvl*/);
 
-/* --- @MPX_UMULN@ --- *
+/* --- @mpx_umuln@ --- *
  *
- * Arguments:  @dv, dvl@ = destination vector base and limit
- *             @av, avl@ = multiplicand vector base and limit
- *             @m@ = multiplier
+ * 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
@@ -431,11 +469,16 @@ extern void mpx_umul(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:  @dv, dvl@ = destination/accumulator vector base and limit
- *             @av, avl@ = multiplicand vector base and limit
- *             @m@ = multiplier
+ * 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.
@@ -444,20 +487,21 @@ extern void mpx_umul(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 = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _c;                       \
+    _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc;                      \
     *_dv++ = MPW(_x);                                                  \
-    _c = _x >> MPW_BITS;                                               \
+    _cc = _x >> MPW_BITS;                                              \
   }                                                                    \
-  MPX_UADDN(_dv, _dvl, _c);                                            \
+  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
@@ -488,14 +532,81 @@ extern void mpx_usqr(mpw */*dv*/, mpw */*dvl*/,
  *             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 two words larger than twice the size of the
- *             divisor.
+ *             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*/);
 
+/*----- Karatsuba multiplication algorithms -------------------------------*/
+
+/* --- @KARATSUBA_CUTOFF@ --- *
+ *
+ * This is the limiting length for using Karatsuba algorithms.  It's best to
+ * use the simpler classical multiplication method on numbers smaller than
+ * this.
+ */
+
+#define KARATSUBA_CUTOFF 16
+
+/* --- @KARATSUBA_SLOP@ --- *
+ *
+ * The extra number of words required as scratch space by the Karatsuba
+ * routines.  This is a (generous) guess, since the actual amount of space
+ * required is proportional to the recursion depth.
+ */
+
+#define KARATSUBA_SLOP 32
+
+/* --- @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 and scratch buffers must be twice as large as
+ *             the larger argument.  The scratch space must be twice as
+ *             large as the larger argument, plus the magic number
+ *             @KARATSUBA_SLOP@.
+ */
+
+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 twice as large as the argument.  The
+ *             scratch space must be twice as large as the argument, plus
+ *             the magic number @KARATSUBA_SLOP@.
+ */
+
+extern void mpx_ksqr(mpw */*dv*/, mpw */*dvl*/,
+                    const mpw */*av*/, const mpw */*avl*/,
+                    mpw */*sv*/, mpw */*svl*/);
+
 /*----- That's all, folks -------------------------------------------------*/
 
 #ifdef __cplusplus