[u/mdw/catacomb] / lcrand.c

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