Bad bug makes all previous testing worthless.

[storin] / lcrand.c

/* -*-c-*-
 *
 * $Id: lcrand.c,v 1.1 2000/05/21 11:28:30 mdw Exp $
 *
 * Simple linear congruential generator
 *
 * (c) 1999 Straylight/Edgeware
 * (c) 2000 Mark Wooding
 */

/*----- Licensing notice --------------------------------------------------* 
 *
 * Copyright (c) 2000 Mark Wooding
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 
 * 2, Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 
 * 3. The name of the authors may not be used to endorse or promote
 *    products derived from this software without specific prior written
 *    permission.
 * 
 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED
 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN
 * NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 * 
 * Instead of accepting the above terms, you may redistribute and/or modify
 * this software under the terms of either the GNU General Public License,
 * or the GNU Library General Public License, published by the Free
 * Software Foundation; either version 2 of the License, or (at your
 * option) any later version.
 */

/*----- Revision history --------------------------------------------------* 
 *
 * $Log: lcrand.c,v $
 * Revision 1.1  2000/05/21 11:28:30  mdw
 * Initial check-in.
 *
 * --- Past lives (Catacomb) --- *
 *
 * Revision 1.2  1999/12/13 15:34:01  mdw
 * Add support for seeding from a generic pseudorandom source.
 *
 * Revision 1.1  1999/12/10 23:15:27  mdw
 * Noncryptographic random number generator.
 *
 */

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

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

#include "bits.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));
}

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