[u/mdw/catacomb] / mptext.c

/* -*-c-*-
 *
 * $Id: mptext.c,v 1.10 2001/06/16 13:22:39 mdw Exp $
 *
 * Textual representation of multiprecision numbers
 *
 * (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.
 */

/*----- Revision history --------------------------------------------------* 
 *
 * $Log: mptext.c,v $
 * Revision 1.10  2001/06/16 13:22:39  mdw
 * Added fast-track code for binary output bases, and tests.
 *
 * Revision 1.9  2001/02/03 16:05:17  mdw
 * Make flags be unsigned.  Improve the write algorithm: recurse until the
 * parts are one word long and use single-precision arithmetic from there.
 * Fix off-by-one bug when breaking the number apart.
 *
 * Revision 1.8  2000/12/06 20:32:42  mdw
 * Reduce binary bytes (to allow marker bits to be ignored).  Fix error
 * message string a bit.  Allow leading `+' signs.
 *
 * Revision 1.7  2000/07/15 10:01:08  mdw
 * Bug fix in binary input.
 *
 * Revision 1.6  2000/06/25 12:58:23  mdw
 * Fix the derivation of `depth' commentary.
 *
 * Revision 1.5  2000/06/17 11:46:19  mdw
 * New and much faster stack-based algorithm for reading integers.  Support
 * reading and writing binary integers in bases between 2 and 256.
 *
 * 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.
 *
 * Revision 1.2  1999/11/20 22:24:15  mdw
 * Use function versions of MPX_UMULN and MPX_UADDN.
 *
 * Revision 1.1  1999/11/17 18:02:16  mdw
 * New multiprecision integer arithmetic suite.
 *
 */

/*----- Header files ------------------------------------------------------*/

#include <ctype.h>
#include <limits.h>
#include <stdio.h>

#include "mp.h"
#include "mptext.h"
#include "paranoia.h"

/*----- Magical numbers ---------------------------------------------------*/

/* --- Maximum recursion depth --- *
 *
 * This is the number of bits in a @size_t@ object.  Why? 
 *
 * To see this, let %$b = \mathit{MPW\_MAX} + 1$% and let %$Z$% be the
 * largest @size_t@ value.  Then the largest possible @mp@ is %$M - 1$% where
 * %$M = b^Z$%.  Let %$r$% be a radix to read or write.  Since the recursion
 * squares the radix at each step, the highest number reached by the
 * recursion is %$d$%, where:
 *
 *   %$r^{2^d} = b^Z$%.
 *
 * Solving gives that %$d = \lg \log_r b^Z$%.  If %$r = 2$%, this is maximum,
 * so choosing %$d = \lg \lg b^Z = \lg (Z \lg b) = \lg Z + \lg \lg b$%.
 *
 * Expressing %$\lg Z$% as @CHAR_BIT * sizeof(size_t)@ yields an
 * overestimate, since a @size_t@ representation may contain `holes'.
 * Choosing to represent %$\lg \lg b$% by 10 is almost certainly sufficient
 * for `some time to come'.
 */

#define DEPTH (CHAR_BIT * sizeof(size_t) + 10)

/*----- Main code ---------------------------------------------------------*/

/* --- @mp_read@ --- *
 *
 * Arguments:   @mp *m@ = destination multiprecision number
 *              @int radix@ = base to assume for data (or zero to guess)
 *              @const mptext_ops *ops@ = pointer to operations block
 *              @void *p@ = data for the operations block
 *
 * Returns:     The integer read, or zero if it didn't work.
 *
 * Use:         Reads an integer from some source.  If the @radix@ is
 *              specified, the number is assumed to be given in that radix,
 *              with the letters `a' (either upper- or lower-case) upwards
 *              standing for digits greater than 9.  Otherwise, base 10 is
 *              assumed unless the number starts with `0' (octal), `0x' (hex)
 *              or `nnn_' (base `nnn').  An arbitrary amount of whitespace
 *              before the number is ignored.
 */

/* --- About the algorithm --- *
 *
 * The algorithm here is rather aggressive.  I maintain an array of
 * successive squarings of the radix, and a stack of partial results, each
 * with a counter attached indicating which radix square to multiply by.
 * Once the item at the top of the stack reaches the same counter level as
 * the next item down, they are combined together and the result is given a
 * counter level one higher than either of the results.
 *
 * Gluing the results together at the end is slightly tricky.  Pay attention
 * to the code.
 *
 * This is more complicated because of the need to handle the slightly
 * bizarre syntax.
 */

mp *mp_read(mp *m, int radix, const mptext_ops *ops, void *p)
{
  int ch;                               /* Current char being considered */
  unsigned f = 0;                       /* Flags about the current number */
  int r;                                /* Radix to switch over to */
  mpw rd;                               /* Radix as an @mp@ digit */
  mp rr;                                /* The @mp@ for the radix */
  unsigned nf = m ? m->f & MP_BURN : 0; /* New @mp@ flags */

  /* --- Stacks --- */

  mp *pow[DEPTH];                       /* List of powers */
  unsigned pows;                        /* Next index to fill */
  struct { unsigned i; mp *m; } s[DEPTH]; /* Main stack */
  unsigned sp;                          /* Current stack pointer */

  /* --- Flags --- */

#define f_neg 1u
#define f_ok 2u
#define f_start 4u

  /* --- Initialize the stacks --- */

  mp_build(&rr, &rd, &rd + 1);
  pow[0] = &rr;  
  pows = 1;

  sp = 0;

  /* --- Initialize the destination number --- */

  if (m)
    MP_DROP(m);

  /* --- Read an initial character --- */

  ch = ops->get(p);
  while (isspace(ch))
    ch = ops->get(p);

  /* --- Handle an initial sign --- */

  if (radix >= 0 && (ch == '-' || ch == '+')) {
    if (ch == '-')
      f |= f_neg;
    do ch = ops->get(p); while isspace(ch);
  }

  /* --- If the radix is zero, look for leading zeros --- */

  if (radix > 0) {
    assert(((void)"ascii radix must be <= 36", radix <= 36));
    rd = radix;
    r = -1;
  } else if (radix < 0) {
    rd = -radix;
    assert(((void)"binary radix must fit in a byte", rd < UCHAR_MAX));
    r = -1;
  } else if (ch != '0') {
    rd = 10;
    r = 0;
  } else {
    ch = ops->get(p);
    if (ch == 'x') {
      ch = ops->get(p);
      rd = 16;
    } else {
      rd = 8;
      f |= f_ok;
    }
    r = -1;
  }

  /* --- Use fast algorithm for binary radix --- *
   *
   * This is the restart point after having parsed a radix number from the
   * input.  We check whether the radix is binary, and if so use a fast
   * algorithm which just stacks the bits up in the right order.
   */

restart:
  switch (rd) {
    unsigned bit;

    case   2: bit = 1; goto bin;
    case   4: bit = 2; goto bin;
    case   8: bit = 3; goto bin;
    case  16: bit = 4; goto bin;
    case  32: bit = 5; goto bin;
    case  64: bit = 6; goto bin;
    case 128: bit = 7; goto bin;
    default:
      break;

  /* --- The fast binary algorithm --- *
   *
   * We stack bits up starting at the top end of a word.  When one word is
   * full, we write it to the integer, and start another with the left-over
   * bits.  When the array in the integer is full, we resize using low-level
   * calls and copy the current data to the top end.  Finally, we do a single
   * bit-shift when we know where the end of the number is.
   */

  bin: {
    mpw a = 0;
    unsigned b = MPW_BITS;
    size_t len, n;
    mpw *v;

    m = mp_dest(MP_NEW, 1, nf);
    len = n = m->sz;
    n = len;
    v = m->v + n;
    for (;; ch = ops->get(p)) {
      unsigned x;

      if (ch < 0)
        break;

      /* --- Check that the character is a digit and in range --- */

      if (radix < 0)
        x = ch % rd;
      else {
        if (!isalnum(ch))
          break;
        if (ch >= '0' && ch <= '9')
          x = ch - '0';
        else {
          ch = tolower(ch);
          if (ch >= 'a' && ch <= 'z')   /* ASCII dependent! */
            x = ch - 'a' + 10;
          else
            break;
        }
      }
      if (x >= rd)
        break;

      /* --- Feed the digit into the accumulator --- */

      f |= f_ok;
      if (!x && !(f & f_start))
        continue;
      f |= f_start;
      if (b > bit) {
        b -= bit;
        a |= MPW(x) << b;
      } else {
        a |= MPW(x) >> (bit - b);
        b += MPW_BITS - bit;
        *--v = MPW(a);
        n--;
        if (!n) {
          n = len;
          len <<= 1;
          v = mpalloc(m->a, len);
          memcpy(v + n, m->v, MPWS(n));
          mpfree(m->a, m->v);
          m->v = v;
          v = m->v + n;
        }
        a = (b < MPW_BITS) ? MPW(x) << b : 0;
      }
    }

    /* --- Finish up --- */

    if (!(f & f_ok)) {
      mp_drop(m);
      m = 0;
    } else {
      *--v = MPW(a);
      n--;
      m->sz = len;
      m->vl = m->v + len;
      m->f &= ~MP_UNDEF;
      m = mp_lsr(m, m, (unsigned long)n * MPW_BITS + b);
    }
    goto done;
  }}

  /* --- Time to start --- */

  for (;; ch = ops->get(p)) {
    unsigned x;

    if (ch < 0)
      break;

    /* --- An underscore indicates a numbered base --- */

    if (ch == '_' && r > 0 && r <= 36) {
      unsigned i;

      /* --- Clear out the stacks --- */

      for (i = 1; i < pows; i++)
        MP_DROP(pow[i]);
      pows = 1;
      for (i = 0; i < sp; i++)
        MP_DROP(s[i].m);
      sp = 0;

      /* --- Restart the search --- */

      rd = r;
      r = -1;
      f &= ~f_ok;
      ch = ops->get(p);
      goto restart;
    }

    /* --- Check that the character is a digit and in range --- */

    if (radix < 0)
      x = ch % rd;
    else {
      if (!isalnum(ch))
        break;
      if (ch >= '0' && ch <= '9')
        x = ch - '0';
      else {
        ch = tolower(ch);
        if (ch >= 'a' && ch <= 'z')     /* ASCII dependent! */
          x = ch - 'a' + 10;
        else
          break;
      }
    }

    /* --- Sort out what to do with the character --- */

    if (x >= 10 && r >= 0)
      r = -1;
    if (x >= rd)
      break;

    if (r >= 0)
      r = r * 10 + x;

    /* --- Stick the character on the end of my integer --- */

    assert(((void)"Number is too unimaginably huge", sp < DEPTH));
    s[sp].m = m = mp_new(1, nf);
    m->v[0] = x;
    s[sp].i = 0;

    /* --- Now grind through the stack --- */

    while (sp > 0 && s[sp - 1].i == s[sp].i) {

      /* --- Combine the top two items --- */

      sp--;
      m = s[sp].m;
      m = mp_mul(m, m, pow[s[sp].i]);
      m = mp_add(m, m, s[sp + 1].m);
      s[sp].m = m;
      MP_DROP(s[sp + 1].m);
      s[sp].i++;

      /* --- Make a new radix power if necessary --- */

      if (s[sp].i >= pows) {
        assert(((void)"Number is too unimaginably huge", pows < DEPTH));
        pow[pows] = mp_sqr(MP_NEW, pow[pows - 1]);
        pows++;
      }
    }
    f |= f_ok;
    sp++;
  }

  ops->unget(ch, p);

  /* --- If we're done, compute the rest of the number --- */

  if (f & f_ok) {
    if (!sp)
      return (MP_ZERO);
    else {
      mp *z = MP_ONE;
      sp--;

      while (sp > 0) {

        /* --- Combine the top two items --- */

        sp--;
        m = s[sp].m;
        z = mp_mul(z, z, pow[s[sp + 1].i]);
        m = mp_mul(m, m, z);
        m = mp_add(m, m, s[sp + 1].m);
        s[sp].m = m;
        MP_DROP(s[sp + 1].m);

        /* --- Make a new radix power if necessary --- */

        if (s[sp].i >= pows) {
          assert(((void)"Number is too unimaginably huge", pows < DEPTH));
          pow[pows] = mp_sqr(MP_NEW, pow[pows - 1]);
          pows++;
        }
      }
      MP_DROP(z);
      m = s[0].m;
    }
  } else {
    unsigned i;
    for (i = 0; i < sp; i++)
      MP_DROP(s[i].m);
  }

  /* --- Clear the radix power list --- */

  {
    unsigned i;
    for (i = 1; i < pows; i++)
      MP_DROP(pow[i]);
  }

  /* --- Bail out if the number was bad --- */

done:
  if (!(f & f_ok))    
    return (0);

  /* --- Set the sign and return --- */

  if (f & f_neg)
    m->f |= MP_NEG;
  return (m);

#undef f_start
#undef f_neg
#undef f_ok
}

/* --- @mp_write@ --- *
 *
 * Arguments:   @mp *m@ = pointer to a multi-precision integer
 *              @int radix@ = radix to use when writing the number out
 *              @const mptext_ops *ops@ = pointer to an operations block
 *              @void *p@ = data for the operations block
 *
 * Returns:     Zero if it worked, nonzero otherwise.
 *
 * Use:         Writes a large integer in textual form.
 */

/* --- Simple case --- *
 *
 * Use a fixed-sized buffer and single-precision arithmetic 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(mpw n, int radix, unsigned z,
                  const mptext_ops *ops, void *p)
{
  int rc = 0;
  char buf[64];
  unsigned i = sizeof(buf);
  int rd = radix > 0 ? radix : -radix;

  do {
    int ch;
    mpw x;

    x = n % rd;
    n /= rd;
    if (radix < 0)
      ch = x;
    else if (x < 10)
      ch = '0' + x;
    else
      ch = 'a' + x - 10;
    buf[--i] = ch;
    if (z)
      z--;
  } while (i && n);

  if (n)
    rc = simple(n, radix, z, ops, p);
  else {
    char zbuf[32];
    memset(zbuf, (radix < 0) ? 0 : '0', sizeof(zbuf));
    while (!rc && z >= sizeof(zbuf)) {
      rc = ops->put(zbuf, sizeof(zbuf), p);
      z -= sizeof(zbuf);
    }
    if (!rc && z)
      rc = ops->put(zbuf, z, p);
  }
  if (!rc)
    rc = ops->put(buf + i, sizeof(buf) - i, p);
  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) < 2)
    return (simple(MP_LEN(m) ? m->v[0] : 0, radix, z, ops, p));

  assert(i);
  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);
}

/* --- Binary case --- *
 *
 * Special case for binary output.  Goes much faster.
 */

static int binary(mp *m, int bit, int radix, const mptext_ops *ops, void *p)
{
  mpw *v;
  mpw a;
  int rc = 0;
  unsigned b;
  unsigned mask;
  unsigned long n;
  unsigned f = 0;
  char buf[8], *q;
  unsigned x;
  int ch;

#define f_out 1u

  /* --- Work out where to start --- */

  n = mp_bits(m);
  n += bit - (n % bit);
  b = n % MPW_BITS;
  n /= MPW_BITS;
  
  if (n > MP_LEN(m)) {
    n--;
    b += MPW_BITS;
  }

  v = m->v + n;
  a = *v;
  mask = (1 << bit) - 1;
  q = buf;

  /* --- Main code --- */

  for (;;) {
    if (b > bit) {
      b -= bit;
      x = a >> b;
    } else {
      x = a << (bit - b);
      b += MPW_BITS - bit;
      if (v == m->v)
        break;
      a = *--v;
      if (b < MPW_BITS)
        x |= a >> b;
    }
    x &= mask;
    if (!x && !(f & f_out))
      continue;

    if (radix < 0)
      ch = x;
    else if (x < 10)
      ch = '0' + x;
    else
      ch = 'a' + x - 10;
    *q++ = ch;
    if (q >= buf + sizeof(buf)) {
      if ((rc = ops->put(buf, sizeof(buf), p)) != 0)
        goto done;
      q = buf;
    }
    f |= f_out;
  }

  x &= mask;
  if (radix < 0)
    ch = x;
  else if (x < 10)
    ch = '0' + x;
  else
    ch = 'a' + x - 10;
  *q++ = ch;
  rc = ops->put(buf, q - buf, p);

done:
  mp_drop(m);
  return (rc);

#undef f_out
}

/* --- Main driver code --- */

int mp_write(mp *m, int radix, const mptext_ops *ops, void *p)
{
  int rc;

  /* --- Set various things up --- */

  m = MP_COPY(m);
  MP_SPLIT(m);

  /* --- Check the radix for sensibleness --- */

  if (radix > 0)
    assert(((void)"ascii radix must be <= 36", radix <= 36));
  else if (radix < 0)
    assert(((void)"binary radix must fit in a byte", -radix < UCHAR_MAX));
  else
    assert(((void)"radix can't be zero in mp_write", 0));

  /* --- If the number is negative, sort that out --- */

  if (m->f & MP_NEG) {
    if (ops->put("-", 1, p))
      return (EOF);
    m->f &= ~MP_NEG;
  }

  /* --- Handle binary radix --- */

  switch (radix) {
    case   2: case   -2: return (binary(m, 1, radix, ops, p));
    case   4: case   -4: return (binary(m, 2, radix, ops, p));
    case   8: case   -8: return (binary(m, 3, radix, ops, p));
    case  16: case  -16: return (binary(m, 4, radix, ops, p));
    case  32: case  -32: return (binary(m, 5, radix, ops, p));
              case  -64: return (binary(m, 6, radix, ops, p));
              case -128: return (binary(m, 7, radix, ops, p));
  }

  /* --- If the number is small, do it the easy way --- */

  if (MP_LEN(m) < 2)
    rc = simple(MP_LEN(m) ? m->v[0] : 0, 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.
   */

  else {
    mp *pr[DEPTH];
    size_t target = (MP_LEN(m) + 1) / 2;
    unsigned i = 0;
    mp *z = mp_new(1, 0);

    /* --- Set up the exponent table --- */

    z->v[0] = (radix > 0 ? radix : -radix);
    z->f = 0;
    for (;;) {
      assert(((void)"Number is too unimaginably huge", i < DEPTH));
      pr[i++] = z;
      if (MP_LEN(z) > target)
        break;
      z = mp_sqr(MP_NEW, z);
    }

    /* --- Write out the answer --- */

    rc = complicated(m, radix, pr, i - 1, 0, ops, p);

    /* --- Tidy away the array --- */

    while (i > 0)
      mp_drop(pr[--i]);
  }

  /* --- Tidying up code --- */

  MP_DROP(m);
  return (rc);
}

/*----- Test rig ----------------------------------------------------------*/

#ifdef TEST_RIG

#include <mLib/testrig.h>

static int verify(dstr *v)
{
  int ok = 1;
  int ib = *(int *)v[0].buf, ob = *(int *)v[2].buf;
  dstr d = DSTR_INIT;
  mp *m = mp_readdstr(MP_NEW, &v[1], 0, ib);
  if (m) {
    if (!ob) {
      fprintf(stderr, "*** unexpected successful parse\n"
                      "*** input [%2i] =     ", ib);
      if (ib < 0)
        type_hex.dump(&v[1], stderr);
      else
        fputs(v[1].buf, stderr);
      mp_writedstr(m, &d, 10);
      fprintf(stderr, "\n*** (value = %s)\n", d.buf);
      ok = 0;
    } else {
      mp_writedstr(m, &d, ob);
      if (d.len != v[3].len || memcmp(d.buf, v[3].buf, d.len) != 0) {
        fprintf(stderr, "*** failed read or write\n"
                        "*** input [%2i]      = ", ib);
        if (ib < 0)
          type_hex.dump(&v[1], stderr);
        else
          fputs(v[1].buf, stderr);
        fprintf(stderr, "\n*** output [%2i]     = ", ob);
        if (ob < 0)
          type_hex.dump(&d, stderr);
        else
          fputs(d.buf, stderr);
        fprintf(stderr, "\n*** expected [%2i]   = ", ob);
        if (ob < 0)
          type_hex.dump(&v[3], stderr);
        else
          fputs(v[3].buf, stderr);
        fputc('\n', stderr);
        ok = 0;
      }
    }
    mp_drop(m);
  } else {
    if (ob) {
      fprintf(stderr, "*** unexpected parse failure\n"
                      "*** input [%i]    = ", ib);
      if (ib < 0)
        type_hex.dump(&v[1], stderr);
      else
        fputs(v[1].buf, stderr);
      fprintf(stderr, "\n*** expected [%i]   = ", ob);
      if (ob < 0)
        type_hex.dump(&v[3], stderr);
      else
        fputs(v[3].buf, stderr);
      fputc('\n', stderr);
      ok = 0;
    }
  }

  dstr_destroy(&d);
  assert(mparena_count(MPARENA_GLOBAL) == 0);
  return (ok);
}

static test_chunk tests[] = {
  { "mptext-ascii", verify,
    { &type_int, &type_string, &type_int, &type_string, 0 } },
  { "mptext-bin-in", verify,
    { &type_int, &type_hex, &type_int, &type_string, 0 } },
  { "mptext-bin-out", verify,
    { &type_int, &type_string, &type_int, &type_hex, 0 } },
  { 0, 0, { 0 } }
};

int main(int argc, char *argv[])
{
  sub_init();
  test_run(argc, argv, tests, SRCDIR "/tests/mptext");
  return (0);
}

#endif

/*----- That's all, folks -------------------------------------------------*/