[u/mdw/catacomb] / lcrand.c

/* -*-c-*-
 *
 * $Id: lcrand.c,v 1.5 2004/04/08 01:36:15 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.
 */

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