--- /dev/null
+/* -*-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/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"/t/lcrand");
+ return (0);
+}
+
+#endif
+
+/*----- That's all, folks -------------------------------------------------*/
projects / u / mdw / catacomb / blobdiff