progs/perftest.c: Use from Glibc syscall numbers.

[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 <mLib/macros.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.
 */

static int char_digit(int ch, int radix)
{
  int r = radix < 0 ? -radix : radix;
  int d;

  if (ch < 0) return (-1);
  if (radix < 0) d = ch;
  else if ('0' <= ch && ch <= '9') d = ch - '0';
  else if ('a' <= ch && ch <= 'z') d = ch - 'a' + 10;
  else if ('A' <= ch && ch <= 'Z') d = ch - 'A' + (radix > 36 ? 36 : 10);
  else return (-1);
  if (d >= r) return (-1);
  return (d);
}

static mp *read_binary(int radix, unsigned bit, unsigned nf,
                       const mptext_ops *ops, void *p)
{
  mpw a = 0;
  unsigned b = MPW_BITS;
  int any = 0, nz = 0;
  int ch, d;
  size_t len, n;
  mpw *v;
  mp *m;

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

  m = mp_dest(MP_NEW, 1, nf);
  len = n = m->sz;
  n = len;
  v = m->v + n;

  for (;;) {
    ch = ops->get(p);
    if ((d = char_digit(ch, radix)) < 0) break;

    /* --- Ignore leading zeroes, but notice that the number is valid --- */

    any = 1;
    if (!d && !nz) continue;
    nz = 1;

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

    if (b > bit) {
      b -= bit;
      a |= MPW(d) << b;
    } else {
      a |= MPW(d) >> (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(d) << b : 0;
    }
  }

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

  ops->unget(ch, p);
  if (!any) { mp_drop(m); return (0); }

  *--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);

  return (m);
}

struct readstate {

  /* --- State for the general-base reader --- *
   *
   * There are two arrays.  The @pow@ array is set so that @pow[i]@ contains
   * %$R^{2^i}$% for @i < pows@.  The stack @s@ contains partial results:
   * each entry contains a value @m@ corresponding to %$2^i$% digits.
   * Inductively, an empty stack represents zero; if a stack represents %$x$%
   * then pushing a new entry on the top causes the stack to represent
   * %$R^{2^i} x + m$%.
   *
   * It is an invariant that each entry has a strictly smaller @i@ than the
   * items beneath it.  This is achieved by coaslescing entries at the top if
   * they have equal %$i$% values: if the top items are %$(m, i)$%, and
   * %$(M', i)$%, and the rest of the stack represents the integer %$x$%,
   * then %$R^{2^i} (R^{2^i} x + M) + m = R^{2^{i+1}} x + (R^{2^i} M + m)$%,
   * so we replace the top two items by %$((R^{2^i} M + m), i + 1)$%, and
   * repeat if necessary.
   */

  unsigned pows, sp;
  struct { unsigned i; mp *m; } s[DEPTH];
  mp *pow[DEPTH];
};

static void ensure_power(struct readstate *rs)
{
  /* --- Make sure we have the necessary %$R^{2^i}$% computed --- */

  if (rs->s[rs->sp].i >= rs->pows) {
    assert(rs->pows < DEPTH);
    rs->pow[rs->pows] = mp_sqr(MP_NEW, rs->pow[rs->pows - 1]);
    rs->pows++;
  }
}

static void read_digit(struct readstate *rs, unsigned nf, int d)
{
  mp *m = mp_new(1, nf);
  m->v[0] = d;

  /* --- Put the new digit on top --- */

  assert(rs->sp < DEPTH);
  rs->s[rs->sp].m = m;
  rs->s[rs->sp].i = 0;

  /* --- Restore the stack invariant --- */

  while (rs->sp && rs->s[rs->sp - 1].i <= rs->s[rs->sp].i) {
    assert(rs->sp > 0);
    ensure_power(rs);
    rs->sp--;

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

  /* --- Leave the stack pointer at an empty item --- */

  rs->sp++;
}

static mp *read_general(int radix, unsigned t, unsigned nf,
                        const mptext_ops *ops, void *p)
{
  struct readstate rs;
  unsigned char v[4];
  unsigned i;
  mpw r;
  int any = 0;
  int ch, d;
  mp rr;
  mp *m, *z, *n;

  /* --- Prepare the stack --- */

  r = radix < 0 ? -radix : radix;
  mp_build(&rr, &r, &r + 1);
  rs.pow[0] = &rr;
  rs.pows = 1;
  rs.sp = 0;

  /* --- If we've partially parsed some input then feed it in --- *
   *
   * Unfortunately, what we've got is backwards.  Fortunately there's a
   * fairly tight upper bound on how many digits @t@ might be, since we
   * aborted that loop once it got too large.
   */

  if (t) {
    i = 0;
    while (t) { assert(i < sizeof(v)); v[i++] = t%r; t /= r; }
    while (i) read_digit(&rs, nf, v[--i]);
    any = 1;
  }

  /* --- Read more stuff --- */

  for (;;) {
    ch = ops->get(p);
    if ((d = char_digit(ch, radix)) < 0) break;
    read_digit(&rs, nf, d); any = 1;
  }
  ops->unget(ch, p);

  /* --- Stitch all of the numbers together --- *
   *
   * This is not the same code as @read_digit@.  In particular, here we must
   * cope with the partial result being some inconvenient power of %$R$%,
   * rather than %$R^{2^i}$%.
   */

  if (!any) return (0);
  m = MP_ZERO; z = MP_ONE;
  while (rs.sp) {
    rs.sp--;
    ensure_power(&rs);
    n = rs.s[rs.sp].m;
    n = mp_mul(n, n, z);
    m = mp_add(m, m, n);
    z = mp_mul(z, z, rs.pow[rs.s[rs.sp].i]);
    MP_DROP(n);
  }
  for (i = 0; i < rs.pows; i++) MP_DROP(rs.pow[i]);
  MP_DROP(z);
  return (m);
}

mp *mp_read(mp *m, int radix, const mptext_ops *ops, void *p)
{
  unsigned t = 0;
  unsigned nf = 0;
  int ch, d, rd;

  unsigned f = 0;
#define f_neg 1u
#define f_ok 2u

  /* --- We don't actually need a destination so throw it away --- *
   *
   * But note the flags before we lose it entirely.
   */

  if (m) {
    nf = m->f & MP_BURN;
    MP_DROP(m);
  }

  /* --- Maintain a lookahead character --- */

  ch = ops->get(p);

  /* --- If we're reading text, skip leading space, and maybe a sign --- */

  if (radix >= 0) {
    while (ISSPACE(ch)) ch = ops->get(p);
    switch (ch) {
      case '-': f |= f_neg; /* and on */
      case '+': do ch = ops->get(p); while (ISSPACE(ch));
    }
  }

  /* --- If we don't have a fixed radix, then parse one from the input --- *
   *
   * This is moderately easy if the input starts with `0x' or similar.  If it
   * starts with `0' and something else, then it might be octal, or just a
   * plain old zero.  Finally, it might start with a leading `NN_', in which
   * case we carefully collect the decimal number until we're sure it's
   * either a radix prefix (in which case we accept it and start over) or it
   * isn't (in which case it's actually the start of a large number we need
   * to read).
   */

  if (radix == 0) {
    if (ch == '0') {
      ch = ops->get(p);
      switch (ch) {
        case 'x': case 'X': radix = 16; goto fetch;
        case 'o': case 'O': radix = 8; goto fetch;
        case 'b': case 'B': radix = 2; goto fetch;
        fetch: ch = ops->get(p); break;
        default: radix = 8; f |= f_ok; break;
      }
    } else {
      if ((d = char_digit(ch, 10)) < 0) { ops->unget(ch, p); return (0); }
      for (;;) {
        t = 10*t + d;
        ch = ops->get(p);
        if (t > 52) break;
        if ((d = char_digit(ch, 10)) < 0) break;
      }
      if (ch != '_' || t > 52) radix = 10;
      else {
        radix = t; t = 0;
        ch = ops->get(p);
      }
    }
  }

  /* --- We're now ready to dispatch to the correct handler --- */

  rd = radix < 0 ? -radix : radix;
  ops->unget(ch, p);
  switch (rd) {
    case   2: m = read_binary(radix,  1, nf, ops, p); break;
    case   4: m = read_binary(radix,  2, nf, ops, p); break;
    case   8: m = read_binary(radix,  3, nf, ops, p); break;
    case  16: m = read_binary(radix,  4, nf, ops, p); break;
    case  32: m = read_binary(radix,  5, nf, ops, p); break;
    case  64: m = read_binary(radix,  6, nf, ops, p); break;
    case 128: m = read_binary(radix,  7, nf, ops, p); break;
    default:  m = read_general(radix, t, nf, ops, p); break;
  }

  /* --- That didn't work --- *
   *
   * If we've already read something then return that.  Otherwise it's an
   * error.
   */

  if (!m) {
    if (f & f_ok) return (MP_ZERO);
    else return (0);
  }

  /* --- Negate the result if we should do that --- */

  if (f & f_neg) m = mp_neg(m, m);

  /* --- And we're all done --- */

  return (m);

#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)) {
        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)) {
    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 -------------------------------------------------*/