[u/mdw/catacomb] / lcrand.c

/* -*-c-*-
 *
 * $Id: lcrand.c,v 1.3 2000/06/17 11:29:03 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.3  2000/06/17 11:29:03  mdw
 * Add the flags word to the generic generator.
 *
 * 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 <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:
	case GRAND_SEEDRAND:
	  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;
    case GRAND_SEEDRAND: {
      grand *rr = va_arg(ap, grand *);
      uint32 x;
      do x = rr->ops->word(rr); while (x >= P || x == LCRAND_FIXEDPT);
      g->x = x;
    } 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, 0,
  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 -------------------------------------------------*/
Commit	Line	Data
07871354	1	/* --c--
07871354	2	*
ec4825ef	3	* $Id: lcrand.c,v 1.3 2000/06/17 11:29:03 mdw Exp $
07871354	4	*
	5	* Simple linear congruential generator
	6	*
	7	* (c) 1999 Straylight/Edgeware
	8	*/
	9
	10	/----- Licensing notice --------------------------------------------------
	11	*
	12	* This file is part of Catacomb.
	13	*
	14	* Catacomb is free software; you can redistribute it and/or modify
	15	* it under the terms of the GNU Library General Public License as
	16	* published by the Free Software Foundation; either version 2 of the
	17	* License, or (at your option) any later version.
	18	*
	19	* Catacomb is distributed in the hope that it will be useful,
	20	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	21	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	22	* GNU Library General Public License for more details.
	23	*
	24	* You should have received a copy of the GNU Library General Public
	25	* License along with Catacomb; if not, write to the Free
	26	* Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
	27	* MA 02111-1307, USA.
	28	*/
	29
	30	/----- Revision history --------------------------------------------------
	31	*
	32	* $Log: lcrand.c,v $
ec4825ef	33	* Revision 1.3 2000/06/17 11:29:03 mdw
	34	* Add the flags word to the generic generator.
	35	*
4ab1268f	36	* Revision 1.2 1999/12/13 15:34:01 mdw
	37	* Add support for seeding from a generic pseudorandom source.
	38	*
07871354	39	* Revision 1.1 1999/12/10 23:15:27 mdw
	40	* Noncryptographic random number generator.
	41	*
	42	*/
	43
	44	/----- Header files ------------------------------------------------------/
	45
	46	#include <stdarg.h>
	47	#include <stdio.h>
	48	#include <stdlib.h>
	49	#include <string.h>
	50
	51	#include <mLib/bits.h>
	52	#include <mLib/sub.h>
	53
	54	#include "grand.h"
	55	#include "lcrand.h"
	56
	57	/----- Magic numbers -----------------------------------------------------/
	58
	59	/* --- The generator parameters --- */
	60
	61	#define P LCRAND_P /* Modulus */
	62	#define A LCRAND_A /* Multiplier (primitive mod @p@) */
	63	#define C LCRAND_C /* Additive constant */
	64
	65	/* --- Precomputed values for modular reduction --- */
	66
	67	#define D 5 /* %$p = 2^{32} - d$% */
	68
	69	/* --- Other useful bits --- */
	70
	71	#define P256 4294967040u /* Highest multiple of 256 < %$p$% */
	72
	73	/----- Main code ---------------------------------------------------------/
	74
	75	/* --- @lcrand@ --- *
	76	*
	77	* Arguments: @uint32 x@ = seed value
	78	*
	79	* Returns: New state of the generator.
	80	*
	81	* Use: Steps the generator. Returns %$ax + c \bmod p$%.
	82	*/
	83
	84	uint32 lcrand(uint32 x)
	85	{
	86	uint32 a[2], xx[2];
	87	uint32 yy[2];
	88
	89	/* --- Unpack things into the arrays --- */
	90
	91	a[0] = U16(A); a[1] = U16(A >> 16);
	92	xx[0] = U16(x); xx[1] = U16(x >> 16);
	93
	94	/* --- Multiply everything together --- *
	95	*
	96	* This is plain old long multiplication, although it looks a bit strange.
	97	* I set up the top and bottom partial products directly where they're
	98	* supposed to be. The cross terms I add together, with the low 16 bits in
	99	* @q@ and the high 32 bits in @p@. These I then add into the product.
	100	*/
	101
	102	{
103	uint32 p, q;
104
105	yy[0] = a[0] * xx[0];
106	yy[1] = a[1] * xx[1];
107
108	p = a[0] * xx[1];
109	q = p + a[1] * xx[0];
110	p = ((q < p) << 16) + (q >> 16);
111	q = U16(q) << 16;
112
113	q += yy[0];
114	if (q < yy[0])
115	p++;
116	else
117	p += (q >> 16) >> 16;
118	yy[0] = q;
119
120	yy[1] += p;
121	}
122
123	/* --- Now reduce mod p --- *
124	*
125	* I'm using shifts and adds to do the multiply step here. This needs to
126	* be changed if @D@ ever becomes something other than 5.
127	*/
128
129	#if D != 5
130	# error "Change shift sequence!"
131	#endif
132
133	{
134	uint32 q;
135
136	q = yy[1];
137	x = yy[0];
138
139	while (q) {
140	uint32 y, z;
141	y = q >> 30;
142	z = q << 2;
143	z += q;
144	if (z < q)
145	y++;
146	else
147	y += (q >> 16) >> 16;
148	q = y;
149	x += z;
150	if (x < z \|\| x > P)
151	x -= P;
152	}
153	}
154
155	/* --- Now add on the constant --- */
156
157	x += C;
158	if (x < C \|\| x >= P)
159	x -= P;
160
161	/* --- Done --- */
162
163	return (x);
164	}
165
166	/* --- @lcrand_range@ --- *
167	*
168	* Arguments: @uint32 *x@ = pointer to seed value (updated)
169	* @uint32 m@ = limit allowable
170	*
171	* Returns: A uniformly distributed pseudorandom integer in the interval
172	* %$[0, m)$%.
173	*/
174
175	uint32 lcrand_range(uint32 *x, uint32 m)
176	{
177	uint32 xx = *x;
178	uint32 r = P - P % m;
179	do xx = lcrand(xx); while (xx >= r);
180	*x = xx;
181	return (xx / (r / m));
182	}
183
184	/----- Generic interface -------------------------------------------------/
185
186	typedef struct gctx {
187	grand r;
188	uint32 x;
189	} gctx;
190
191	static void gdestroy(grand *r)
192	{
193	gctx g = (gctx )r;
194	DESTROY(g);
195	}
196
197	static int gmisc(grand *r, unsigned op, ...)
198	{
199	gctx g = (gctx )r;
200	va_list ap;
201	int rc = 0;
202	va_start(ap, op);
203
204	switch (op) {
205	case GRAND_CHECK:
206	switch (va_arg(ap, unsigned)) {
207	case GRAND_CHECK:
208	case GRAND_SEEDINT:
209	case GRAND_SEEDUINT32:
4ab1268f	210	case GRAND_SEEDRAND:
07871354	211	rc = 1;
	212	break;
	213	default:
	214	rc = 0;
	215	break;
	216	}
	217	break;
	218	case GRAND_SEEDINT:
	219	g->x = va_arg(ap, unsigned);
	220	break;
	221	case GRAND_SEEDUINT32:
	222	g->x = va_arg(ap, uint32);
	223	break;
4ab1268f	224	case GRAND_SEEDRAND: {
	225	grand rr = va_arg(ap, grand );
	226	uint32 x;
	227	do x = rr->ops->word(rr); while (x >= P \|\| x == LCRAND_FIXEDPT);
	228	g->x = x;
	229	} break;
07871354	230	default:
	231	GRAND_BADOP;
	232	break;
	233	}
	234
	235	va_end(ap);
	236	return (rc);
	237	}
	238
	239	static uint32 graw(grand *r)
	240	{
	241	gctx g = (gctx )r;
	242	g->x = lcrand(g->x);
	243	return (g->x);
	244	}
	245
	246	static octet gbyte(grand *r)
	247	{
	248	gctx g = (gctx )r;
	249	uint32 x = g->x;
	250	do x = lcrand(x); while (x >= P256);
	251	g->x = x;
	252	return (x / (P256 / 256));
	253	}
	254
	255	static uint32 grange(grand *r, uint32 l)
	256	{
	257	gctx g = (gctx )r;
	258	return (lcrand_range(&g->x, l));
	259	}
	260
	261	static const grand_ops gops = {
	262	"lcrand",
ec4825ef	263	LCRAND_P, 0,
07871354	264	gmisc, gdestroy,
	265	graw, gbyte, grand_word, grange, grand_fill
	266	};
	267
	268	/* --- @lcrand_create@ --- *
	269	*
	270	* Arguments: @uint32 x@ = initial seed
	271	*
	272	* Returns: Pointer to a generic generator.
	273	*
	274	* Use: Constructs a generic generator interface over a linear
	275	* congruential generator.
	276	*/
	277
	278	grand *lcrand_create(uint32 x)
	279	{
	280	gctx *g = CREATE(gctx);
	281	g->r.ops = &gops;
	282	g->x = x;
	283	return (&g->r);
	284	}
	285
	286	/----- Test rig ----------------------------------------------------------/
	287
	288	#ifdef TEST_RIG
	289
	290	#include <mLib/testrig.h>
	291
	292	static int verify(dstr *v)
	293	{
	294	uint32 x = (uint32 )v[0].buf;
	295	uint32 y = (uint32 )v[1].buf;
	296	uint32 z = lcrand(x);
	297	int ok = 1;
	298	if (y != z) {
	299	fprintf(stderr,
	300	"\n*** lcrand failed. lcrand(%lu) = %lu, expected %lu\n",
	301	(unsigned long)x, (unsigned long)z, (unsigned long)y);
	302	ok = 0;
	303	}
	304	return (ok);
	305	}
	306
	307	static test_chunk tests[] = {
	308	{ "lcrand", verify, { &type_uint32, &type_uint32, 0 } },
	309	{ 0, 0, { 0 } }
	310	};
	311
	312	int main(int argc, char *argv[])
	313	{
	314	test_run(argc, argv, tests, SRCDIR"/tests/lcrand");
	315	return (0);
	316	}
	317
	318	#endif
	319
	320	/----- That's all, folks -------------------------------------------------/