X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/blobdiff_plain/9c1437f372e62f0b3b3a7162aabee73bdc96ce4b..HEAD:/math/scaf.c

diff --git a/math/scaf.c b/math/scaf.c
index 778e5e34..59f9c6ee 100644
--- a/math/scaf.c
+++ b/math/scaf.c
@@ -31,6 +31,52 @@
 
 #include "scaf.h"
 
+/*----- Debugging utilties ------------------------------------------------*/
+
+#ifdef SCAF_DEBUG
+
+#include <stdio.h>
+
+#include "mp.h"
+#include "mpint.h"
+#include "mptext.h"
+
+static void scaf_dump(const char *what, const scaf_piece *x,
+		      size_t npiece, size_t piecewd)
+{
+  mp *y = MP_ZERO, *t = MP_NEW;
+  size_t i;
+  unsigned o = 0;
+
+  for (i = 0; i < npiece; i++) {
+    t = mp_fromuint64(t, x[i]);
+    t = mp_lsl(t, t, o);
+    y = mp_add(y, y, t);
+    o += piecewd;
+  }
+  printf(";; %s", what); MP_PRINT("", y); putchar('\n');
+  mp_drop(y); mp_drop(t);
+}
+
+static void scaf_dumpdbl(const char *what, const scaf_dblpiece *x,
+		      size_t npiece, size_t piecewd)
+{
+  mp *y = MP_ZERO, *t = MP_NEW;
+  size_t i;
+  unsigned o = 0;
+
+  for (i = 0; i < npiece; i++) {
+    t = mp_fromuint64(t, x[i]);
+    t = mp_lsl(t, t, o);
+    y = mp_add(y, y, t);
+    o += piecewd;
+  }
+  printf(";; %s", what); MP_PRINT("", y); putchar('\n');
+  mp_drop(y); mp_drop(t);
+}
+
+#endif
+
 /*----- Main code ---------------------------------------------------------*/
 
 /* --- @scaf_load@ --- *
@@ -159,7 +205,7 @@ void scaf_mul(scaf_dblpiece *z, const scaf_piece *x, const scaf_piece *y,
  *		@const scaf_piece *mu@ = scaled approximation to @1/l@
  *		@size_t npiece@ = number of pieces in @l@
  *		@unsigned piecewd@ = nominal width of a piece in bits
- *		@scaf_piece *scratch@ = @3*npiece + 1@ scratch pieces
+ *		@scaf_piece *scratch@ = @3*npiece@ scratch pieces
  *
  * Returns:	---
  *
@@ -184,7 +230,7 @@ void scaf_reduce(scaf_piece *z, const scaf_dblpiece *x,
 
   /* This here is the hard part.
    *
-   * Let w = PIECEWD, let n = NPIECE, and let B = 2^w.  Wwe must have
+   * Let w = PIECEWD, let n = NPIECE, and let B = 2^w.  We must have
    * B^(n-1) <= l < B^n.
    *
    * The argument MU contains pieces of the quantity µ = floor(B^2n/l), which
@@ -269,7 +315,7 @@ void scaf_reduce(scaf_piece *z, const scaf_dblpiece *x,
       u += z[j] + (l[j] ^ m);
       t[j] = u&m; u >>= piecewd;
     }
-    for (j = 0, u = -u; j < npiece; j++) z[i] = (t[i]&u) | (z[i]&~u);
+    for (j = 0, u = -u; j < npiece; j++) z[j] = (t[j]&u) | (z[j]&~u);
   }
 }