math/mptext.c: Reformat and refactor output functions.

[catacomb] / math / mptext.c

/* -*-c-*-
 *
 * 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.
 */

/*----- 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 = \textit{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)

/*----- Input -------------------------------------------------------------*/

/* --- @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);
  if (radix >= 0) {
    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 <= 62", radix <= 62));
    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);
    switch (ch) {
      case 'x':
	rd = 16;
	goto prefix;
      case 'o':
	rd = 8;
	goto prefix;
      case 'b':
	rd = 2;
	goto prefix;
      prefix:
	ch = ops->get(p);
	break;
      default:
	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 {
	  if (rd <= 36)
	    ch = tolower(ch);
	  if (ch >= 'a' && ch <= 'z')	/* ASCII dependent! */
	    x = ch - 'a' + 10;
	  else if (ch >= 'A' && ch <= 'Z')
	    x = ch - 'A' + 36;
	  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);
    }
    ops->unget(ch, p);
    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 <= 62) {
      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 {
	if (rd <= 36)
	  ch = tolower(ch);
	if (ch >= 'a' && ch <= 'z')	/* ASCII dependent! */
	  x = ch - 'a' + 10;
	else if (ch >= 'A' && ch <= 'Z')
	  x = ch - 'A' + 36;
	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;
  MP_SHRINK(m);
  return (m);

#undef f_start
#undef f_neg
#undef f_ok
}

/*----- Output ------------------------------------------------------------*/

/* --- @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.
 */

static int digit_char(int d, int radix)
{
  if (radix < 0) return (d);
  else if (d < 10) return (d + '0');
  else if (d < 26) return (d - 10 + 'a');
  else return (d - 36 + 'A');
}

/* --- 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 write_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;
  mpw x;

  do {
    x = n % rd; n /= rd;
    buf[--i] = digit_char(x, radix);
    if (z) z--;
  } while (i && n);

  if (n)
    rc = write_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 write_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 (write_simple(MP_LEN(m) ? m->v[0] : 0, radix, z, ops, p));

  assert(i);
  mp_div(&q, &m, m, pr[i]);
  if (MP_ZEROP(q)) d = z;
  else {
    if (z > d) z -= d;
    else z = 0;
    rc = write_complicated(q, radix, pr, i - 1, z, ops, p);
  }
  if (!rc) rc = write_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 write_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;

#define f_out 1u

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

  n = mp_bits(m);
  if (n % bit) 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;

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

  x &= mask;
  *q++ = digit_char(x, radix);
  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;
  mp *pr[DEPTH];
  size_t target;
  unsigned i = 0;
  mp *z;

  if (MP_EQ(m, MP_ZERO))
    return (ops->put(radix > 0 ? "0" : "\0", 1, p));

  /* --- 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 <= 62", radix <= 62));
  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 (MP_NEGP(m)) {
    assert(radix > 0);
    if (ops->put("-", 1, p)) return (EOF);
    m->f &= ~MP_NEG;
  }

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

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

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

  if (MP_LEN(m) < 2)
    rc = write_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 {
    target = (MP_LEN(m) + 1) / 2;
    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 = write_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;
  size_t off = 0;
  mp *m = mp_readdstr(MP_NEW, &v[1], &off, 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 [%2i]	  = ", ib);
      if (ib < 0)
	type_hex.dump(&v[1], stderr);
      else
	fputs(v[1].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;
    }
  }

  if (v[1].len - off != v[4].len ||
      memcmp(v[1].buf + off, v[4].buf, v[4].len) != 0) {
    fprintf(stderr, "*** leftovers incorrect\n"
		    "*** input [%2i]	= ", ib);
    if (ib < 0)
      type_hex.dump(&v[1], stderr);
    else
      fputs(v[1].buf, stderr);
    fprintf(stderr, "\n*** expected `%s'\n"
		    "*** found `%s'\n",
	    v[4].buf, v[1].buf + off);
    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, &type_string, 0 } },
  { "mptext-bin-in", verify,
    { &type_int, &type_hex, &type_int, &type_string, &type_string, 0 } },
  { "mptext-bin-out", verify,
    { &type_int, &type_string, &type_int, &type_hex, &type_string, 0 } },
  { 0, 0, { 0 } }
};

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

#endif

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