+++ /dev/null
-/* -*-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 -------------------------------------------------*/
projects / u / mdw / catacomb / blobdiff