Use clever recursive algorithm for writing numbers out.

diff --git a/mptext.c b/mptext.c
index b216fea..c5d20e4 100644 (file)
--- a/mptext.c
+++ b/mptext.c
@@ -1,6 +1,6 @@
 /* -*-c-*-
  *
- * $Id: mptext.c,v 1.3 1999/12/10 23:23:26 mdw Exp $
+ * $Id: mptext.c,v 1.4 1999/12/22 15:56:56 mdw Exp $
  *
  * Textual representation of multiprecision numbers
  *
@@ -30,6 +30,9 @@
 /*----- Revision history --------------------------------------------------* 
  *
  * $Log: mptext.c,v $
+ * Revision 1.4  1999/12/22 15:56:56  mdw
+ * Use clever recursive algorithm for writing numbers out.
+ *
  * Revision 1.3  1999/12/10 23:23:26  mdw
  * Allocate slightly less memory.
  *
@@ -46,17 +49,9 @@
 #include <ctype.h>
 #include <stdio.h>
 
-#include <mLib/darray.h>
-
 #include "mp.h"
 #include "mptext.h"
-
-/*----- Data structures ---------------------------------------------------*/
-
-#ifndef CHARV
-   DA_DECL(charv, char);
-#  define CHARV
-#endif
+#include "paranoia.h"
 
 /*----- Main code ---------------------------------------------------------*/
 
@@ -206,54 +201,157 @@ mp *mp_read(mp *m, int radix, const mptext_ops *ops, void *p)
  * Use:                Writes a large integer in textual form.
  */
 
+/* --- Simple case --- *
+ *
+ * Use a fixed-sized buffer and the simple single-precision division
+ * algorithm to pick off low-order digits.  Put each digit in a buffer,
+ * working backwards from the end.  If the buffer becomes full, recurse to
+ * get another one.  Ensure that there are at least @z@ digits by writing
+ * leading zeroes if there aren't enough real digits.
+ */
+
+static int simple(mp *m, int radix, unsigned z,
+                 const mptext_ops *ops, void *p)
+{
+  int rc = 0;
+  char buf[64];
+  unsigned i = sizeof(buf);
+
+  do {
+    int ch;
+    mpw x;
+
+    x = mpx_udivn(m->v, m->vl, m->v, m->vl, radix);
+    MP_SHRINK(m);
+    if (x < 10)
+      ch = '0' + x;
+    else
+      ch = 'a' + x - 10;
+    buf[--i] = ch;
+    if (z)
+      z--;
+  } while (i && MP_LEN(m));
+
+  if (MP_LEN(m))
+    rc = simple(m, radix, z, ops, p);
+  else {
+    static const char zero[32] = "00000000000000000000000000000000";
+    while (!rc && z >= sizeof(zero)) {
+      rc = ops->put(zero, sizeof(zero), p);
+      z -= sizeof(zero);
+    }
+    if (!rc && z)
+      rc = ops->put(zero, z, p);
+  }
+  if (!rc)
+    ops->put(buf + i, sizeof(buf) - i, p);
+  if (m->f & MP_BURN)
+    BURN(buf);
+  return (rc);
+}
+
+/* --- Complicated case --- *
+ *
+ * If the number is small, fall back to the simple case above.  Otherwise
+ * divide and take remainder by current large power of the radix, and emit
+ * each separately.  Don't emit a zero quotient.  Be very careful about
+ * leading zeroes on the remainder part, because they're deeply significant.
+ */
+
+static int complicated(mp *m, int radix, mp **pr, unsigned i, unsigned z,
+                      const mptext_ops *ops, void *p)
+{
+  int rc = 0;
+  mp *q = MP_NEW;
+  unsigned d = 1 << i;
+
+  if (MP_LEN(m) < 8)
+    return (simple(m, radix, z, ops, p));
+
+  mp_div(&q, &m, m, pr[i]);
+  if (!MP_LEN(q))
+    d = z;
+  else {
+    if (z > d)
+      z -= d;
+    else
+      z = 0;
+    rc = complicated(q, radix, pr, i - 1, z, ops, p);
+  }
+  if (!rc)
+    rc = complicated(m, radix, pr, i - 1, d, ops, p);
+  mp_drop(q);
+  return (rc);
+}
+
+/* --- Main driver code --- */
+
 int mp_write(mp *m, int radix, const mptext_ops *ops, void *p)
 {
-  charv v = DA_INIT;
-  mpw b = radix;
-  mp bb;
+  int rc;
 
   /* --- Set various things up --- */
 
   m = MP_COPY(m);
-  mp_build(&bb, &b, &b + 1);
+  MP_SPLIT(m);
 
   /* --- If the number is negative, sort that out --- */
 
   if (m->f & MP_NEG) {
     if (ops->put("-", 1, p))
       return (EOF);
-    MP_SPLIT(m);
-    m->f &= ~MP_NEG;
   }
 
-  /* --- Write digits to a temporary array --- */
-
-  do {
-    mp *q = MP_NEW, *r = MP_NEW;
-    int ch;
-    mpw x;
+  /* --- If the number is small, do it the easy way --- */
+
+  if (MP_LEN(m) < 8)
+    rc = simple(m, radix, 0, ops, p);
+
+  /* --- Use a clever algorithm --- *
+   *
+   * Square the radix repeatedly, remembering old results, until I get
+   * something more than half the size of the number @m@.  Use this to divide
+   * the number: the quotient and remainder will be approximately the same
+   * size, and I'll have split them on a digit boundary, so I can just emit
+   * the quotient and remainder recursively, in order.
+   *
+   * The array size copes with the largest number possibly representable on
+   * the host machine.  Such a large number shouldn't ever arise in real use.
+   */
+
+  else {
+    mp *pr[CHAR_BIT * sizeof(size_t)];
+    size_t target = MP_LEN(m) / 2;
+    unsigned i = 0;
+    mp *z = mp_create(1);
+
+    /* --- Set up the exponent table --- */
+
+    z->v[0] = radix;
+    z->f = 0;
+    for (;;) {
+      assert(((void)"Number is too unimaginably huge",
+             i < sizeof(pr) / sizeof(pr[0])));
+      pr[i++] = z;
+      if (MP_LEN(z) > target)
+       break;
+      z = mp_sqr(MP_NEW, z);
+    }
 
-    mp_div(&q, &r, m, &bb);
-    x = r->v[0];
-    if (x < 10)
-      ch = '0' + x;
-    else
-      ch = 'a' + x - 10;
-    DA_UNSHIFT(&v, ch);
-    MP_DROP(m);
-    MP_DROP(r);
-    m = q;
-  } while (MP_CMP(m, !=, MP_ZERO));
+    /* --- Write out the answer --- */
 
-  /* --- Finished that --- */
+    rc = complicated(m, radix, pr, i - 1, 0, ops, p);
 
-  MP_DROP(m);
+    /* --- Tidy away the array --- */
 
-  {
-    int rc = ops->put(DA(&v), DA_LEN(&v), p);
-    DA_DESTROY(&v);
-    return (rc ? EOF : 0);
+    while (i > 0)
+      mp_drop(pr[--i]);
   }
+
+  /* --- Tidying up code --- */
+
+  MP_DROP(m);
+  return (rc);
 }
 
 /*----- Test rig ----------------------------------------------------------*/