[catacomb] / rand / lcrand.c

/* -*-c-*-
 *
 * 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/macros.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.
   */

  {
    STATIC_ASSERT(D == 5, "Shift sequence doesn't match prime");
    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",
  0, LCRAND_P,
  gmisc, gdestroy,
  graw, gbyte, grand_defaultword, grange, grand_defaultfill
};

/* --- @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"/t/lcrand");
  return (0);
}

#endif

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