Noncryptographic random number generator.

[u/mdw/catacomb] / lcrand.c

/* -*-c-*-
 *
 * $Id: lcrand.c,v 1.1 1999/12/10 23:15:27 mdw Exp $
 *
 * Simple linear congruential generator
 *
 * (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: lcrand.c,v $
 * Revision 1.1  1999/12/10 23:15:27  mdw
 * Noncryptographic random number generator.
 *
 */

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

#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <mLib/bits.h>
#include <mLib/sub.h>

#include "grand.h"
#include "lcrand.h"

/*----- Magic numbers -----------------------------------------------------*/

/* --- The generator parameters --- */

#define P LCRAND_P			/* Modulus */
#define A LCRAND_A			/* Multiplier (primitive mod @p@) */
#define C LCRAND_C			/* Additive constant */

/* --- Precomputed values for modular reduction --- */

#define D 5				/* %$p = 2^{32} - d$% */

/* --- Other useful bits --- */

#define P256 4294967040u		/* Highest multiple of 256 < %$p$% */

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

/* --- @lcrand@ --- *
 *
 * Arguments:	@uint32 x@ = seed value
 *
 * Returns:	New state of the generator.
 *
 * Use:		Steps the generator.  Returns %$ax + c \bmod p$%.
 */

uint32 lcrand(uint32 x)
{
  uint32 a[2], xx[2];
  uint32 yy[2];

  /* --- Unpack things into the arrays --- */

  a[0] = U16(A); a[1] = U16(A >> 16);
  xx[0] = U16(x); xx[1] = U16(x >> 16);

  /* --- Multiply everything together --- *
   *
   * This is plain old long multiplication, although it looks a bit strange.
   * I set up the top and bottom partial products directly where they're
   * supposed to be.  The cross terms I add together, with the low 16 bits in
   * @q@ and the high 32 bits in @p@.  These I then add into the product.
   */

  {
    uint32 p, q;

    yy[0] = a[0] * xx[0];
    yy[1] = a[1] * xx[1];

    p = a[0] * xx[1];
    q = p + a[1] * xx[0];
    p = ((q < p) << 16) + (q >> 16);
    q = U16(q) << 16;

    q += yy[0];
    if (q < yy[0])
      p++;
    else
      p += (q >> 16) >> 16;
    yy[0] = q;

    yy[1] += p;
  }

  /* --- Now reduce mod p --- *
   *
   * I'm using shifts and adds to do the multiply step here.  This needs to
   * be changed if @D@ ever becomes something other than 5.
   */

#if D != 5
#  error "Change shift sequence!"
#endif

  {
    uint32 q;

    q = yy[1];
    x = yy[0];

    while (q) {
      uint32 y, z;
      y = q >> 30;
      z = q << 2;
      z += q;
      if (z < q)
	y++;
      else
	y += (q >> 16) >> 16;
      q = y;
      x += z;
      if (x < z || x > P)
	x -= P;
    }
  }

  /* --- Now add on the constant --- */

  x += C;
  if (x < C || x >= P)
    x -= P;

  /* --- Done --- */

  return (x);
}

/* --- @lcrand_range@ --- *
 *
 * Arguments:	@uint32 *x@ = pointer to seed value (updated)
 *		@uint32 m@ = limit allowable
 *
 * Returns:	A uniformly distributed pseudorandom integer in the interval
 *		%$[0, m)$%.
 */

uint32 lcrand_range(uint32 *x, uint32 m)
{
  uint32 xx = *x;
  uint32 r = P - P % m;
  do xx = lcrand(xx); while (xx >= r);
  *x = xx;
  return (xx / (r / m));
}

/*----- Generic interface -------------------------------------------------*/

typedef struct gctx {
  grand r;
  uint32 x;
} gctx;

static void gdestroy(grand *r)
{
  gctx *g = (gctx *)r;
  DESTROY(g);
}

static int gmisc(grand *r, unsigned op, ...)
{
  gctx *g = (gctx *)r;
  va_list ap;
  int rc = 0;
  va_start(ap, op);

  switch (op) {
    case GRAND_CHECK:
      switch (va_arg(ap, unsigned)) {
	case GRAND_CHECK:
	case GRAND_SEEDINT:
	case GRAND_SEEDUINT32:
	  rc = 1;
	  break;
	default:
	  rc = 0;
	  break;
      }
      break;
    case GRAND_SEEDINT:
      g->x = va_arg(ap, unsigned);
      break;
    case GRAND_SEEDUINT32:
      g->x = va_arg(ap, uint32);
      break;
    default:
      GRAND_BADOP;
      break;
  }

  va_end(ap);
  return (rc);
}

static uint32 graw(grand *r)
{
  gctx *g = (gctx *)r;
  g->x = lcrand(g->x);
  return (g->x);
}

static octet gbyte(grand *r)
{
  gctx *g = (gctx *)r;
  uint32 x = g->x;
  do x = lcrand(x); while (x >= P256);
  g->x = x;
  return (x / (P256 / 256));
}  

static uint32 grange(grand *r, uint32 l)
{
  gctx *g = (gctx *)r;
  return (lcrand_range(&g->x, l));
}

static const grand_ops gops = {
  "lcrand",
  LCRAND_P,
  gmisc, gdestroy,
  graw, gbyte, grand_word, grange, grand_fill
};

/* --- @lcrand_create@ --- *
 *
 * Arguments:	@uint32 x@ = initial seed
 *
 * Returns:	Pointer to a generic generator.
 *
 * Use:		Constructs a generic generator interface over a linear
 *		congruential generator.
 */

grand *lcrand_create(uint32 x)
{
  gctx *g = CREATE(gctx);
  g->r.ops = &gops;
  g->x = x;
  return (&g->r);
}

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

#ifdef TEST_RIG

#include <mLib/testrig.h>

static int verify(dstr *v)
{
  uint32 x = *(uint32 *)v[0].buf;
  uint32 y = *(uint32 *)v[1].buf;
  uint32 z = lcrand(x);
  int ok = 1;
  if (y != z) {
    fprintf(stderr,
	    "\n*** lcrand failed.  lcrand(%lu) = %lu, expected %lu\n",
	    (unsigned long)x, (unsigned long)z, (unsigned long)y);
    ok = 0;
  }
  return (ok);
}

static test_chunk tests[] = {
  { "lcrand", verify, { &type_uint32, &type_uint32, 0 } },
  { 0, 0, { 0 } }
};

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

#endif

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