Pollard's rho algorithm for computing discrete logs.

[u/mdw/catacomb] / mptext.c

/* -*-c-*-
 *
 * $Id: mptext.c,v 1.6 2000/06/25 12:58:23 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.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 --- */

  enum {
    f_neg = 1u,
    f_ok = 2u
  };

  /* --- 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 (ch == '-') {
    f |= f_neg;
    ch = ops->get(p);
    while (isspace(ch))
      ch = ops->get(p);
  }

  /* --- 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;
  }

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

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

    /* --- 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;
      continue;
    }

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

    if (radix < 0)
      x = ch;
    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 --- */

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

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

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

/* --- @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 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);
  int rd = radix > 0 ? radix : -radix;

  do {
    int ch;
    mpw x;

    x = mpx_udivn(m->v, m->vl, m->v, m->vl, rd);
    MP_SHRINK(m);
    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 && 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)
{
  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;
  }

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

  else {
    mp *pr[DEPTH];
    size_t target = MP_LEN(m) / 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 [%i] = ", 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 [%i]    = ", ib);
	if (ib < 0)
	  type_hex.dump(&v[1], stderr);
	else
	  fputs(v[1].buf, stderr);
	fprintf(stderr, "\n*** output [%i]   = ", ob);
	if (ob < 0)
	  type_hex.dump(&d, stderr);
	else
	  fputs(d.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;
      }
    }
    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 -------------------------------------------------*/