3 * Simple linear congruential generator
5 * (c) 1999 Straylight/Edgeware
8 /*----- Licensing notice --------------------------------------------------*
10 * This file is part of Catacomb.
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.
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.
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,
28 /*----- Header files ------------------------------------------------------*/
35 #include <mLib/bits.h>
36 #include <mLib/macros.h>
42 /*----- Magic numbers -----------------------------------------------------*/
44 /* --- The generator parameters --- */
46 #define P LCRAND_P /* Modulus */
47 #define A LCRAND_A /* Multiplier (primitive mod @p@) */
48 #define C LCRAND_C /* Additive constant */
50 /* --- Precomputed values for modular reduction --- */
52 #define D 5 /* %$p = 2^{32} - d$% */
54 /* --- Other useful bits --- */
56 #define P256 4294967040u /* Highest multiple of 256 < %$p$% */
58 /*----- Main code ---------------------------------------------------------*/
62 * Arguments: @uint32 x@ = seed value
64 * Returns: New state of the generator.
66 * Use: Steps the generator. Returns %$ax + c \bmod p$%.
69 uint32 lcrand(uint32 x)
74 /* --- Unpack things into the arrays --- */
76 a[0] = U16(A); a[1] = U16(A >> 16);
77 xx[0] = U16(x); xx[1] = U16(x >> 16);
79 /* --- Multiply everything together --- *
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.
95 p = ((q < p) << 16) + (q >> 16);
102 p += (q >> 16) >> 16;
108 /* --- Now reduce mod p --- *
110 * I'm using shifts and adds to do the multiply step here.
114 STATIC_ASSERT(D == 5, "Shift sequence doesn't match prime");
128 y += (q >> 16) >> 16;
136 /* --- Now add on the constant --- */
147 /* --- @lcrand_range@ --- *
149 * Arguments: @uint32 *x@ = pointer to seed value (updated)
150 * @uint32 m@ = limit allowable
152 * Returns: A uniformly distributed pseudorandom integer in the interval
156 uint32 lcrand_range(uint32 *x, uint32 m)
159 uint32 r = P - P % m;
160 do xx = lcrand(xx); while (xx >= r);
165 /*----- Generic interface -------------------------------------------------*/
167 typedef struct gctx {
172 static void gdestroy(grand *r)
178 static int gmisc(grand *r, unsigned op, ...)
187 switch (va_arg(ap, unsigned)) {
190 case GRAND_SEEDUINT32:
200 g->x = va_arg(ap, unsigned);
202 case GRAND_SEEDUINT32:
203 g->x = va_arg(ap, uint32);
205 case GRAND_SEEDRAND: {
206 grand *rr = va_arg(ap, grand *);
208 do x = rr->ops->word(rr); while (x >= P || x == LCRAND_FIXEDPT);
220 static uint32 graw(grand *r)
227 static octet gbyte(grand *r)
231 do x = lcrand(x); while (x >= P256);
233 return (x / (P256 / 256));
236 static uint32 grange(grand *r, uint32 l)
239 return (lcrand_range(&g->x, l));
242 static const grand_ops gops = {
246 graw, gbyte, grand_defaultword, grange, grand_defaultfill
249 /* --- @lcrand_create@ --- *
251 * Arguments: @uint32 x@ = initial seed
253 * Returns: Pointer to a generic generator.
255 * Use: Constructs a generic generator interface over a linear
256 * congruential generator.
259 grand *lcrand_create(uint32 x)
261 gctx *g = CREATE(gctx);
267 /*----- Test rig ----------------------------------------------------------*/
271 #include <mLib/testrig.h>
273 static int verify(dstr *v)
275 uint32 x = *(uint32 *)v[0].buf;
276 uint32 y = *(uint32 *)v[1].buf;
277 uint32 z = lcrand(x);
281 "\n*** lcrand failed. lcrand(%lu) = %lu, expected %lu\n",
282 (unsigned long)x, (unsigned long)z, (unsigned long)y);
288 static test_chunk tests[] = {
289 { "lcrand", verify, { &type_uint32, &type_uint32, 0 } },
293 int main(int argc, char *argv[])
295 test_run(argc, argv, tests, SRCDIR"/t/lcrand");
301 /*----- That's all, folks -------------------------------------------------*/