3 * $Id: lcrand.c,v 1.4 2000/12/06 20:31:06 mdw Exp $
5 * Simple linear congruential generator
7 * (c) 1999 Straylight/Edgeware
10 /*----- Licensing notice --------------------------------------------------*
12 * This file is part of Catacomb.
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.
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.
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,
30 /*----- Revision history --------------------------------------------------*
33 * Revision 1.4 2000/12/06 20:31:06 mdw
34 * Simplify uniform range transformation.
36 * Revision 1.3 2000/06/17 11:29:03 mdw
37 * Add the flags word to the generic generator.
39 * Revision 1.2 1999/12/13 15:34:01 mdw
40 * Add support for seeding from a generic pseudorandom source.
42 * Revision 1.1 1999/12/10 23:15:27 mdw
43 * Noncryptographic random number generator.
47 /*----- Header files ------------------------------------------------------*/
54 #include <mLib/bits.h>
60 /*----- Magic numbers -----------------------------------------------------*/
62 /* --- The generator parameters --- */
64 #define P LCRAND_P /* Modulus */
65 #define A LCRAND_A /* Multiplier (primitive mod @p@) */
66 #define C LCRAND_C /* Additive constant */
68 /* --- Precomputed values for modular reduction --- */
70 #define D 5 /* %$p = 2^{32} - d$% */
72 /* --- Other useful bits --- */
74 #define P256 4294967040u /* Highest multiple of 256 < %$p$% */
76 /*----- Main code ---------------------------------------------------------*/
80 * Arguments: @uint32 x@ = seed value
82 * Returns: New state of the generator.
84 * Use: Steps the generator. Returns %$ax + c \bmod p$%.
87 uint32
lcrand(uint32 x
)
92 /* --- Unpack things into the arrays --- */
94 a
[0] = U16(A
); a
[1] = U16(A
>> 16);
95 xx
[0] = U16(x
); xx
[1] = U16(x
>> 16);
97 /* --- Multiply everything together --- *
99 * This is plain old long multiplication, although it looks a bit strange.
100 * I set up the top and bottom partial products directly where they're
101 * supposed to be. The cross terms I add together, with the low 16 bits in
102 * @q@ and the high 32 bits in @p@. These I then add into the product.
108 yy
[0] = a
[0] * xx
[0];
109 yy
[1] = a
[1] * xx
[1];
112 q
= p
+ a
[1] * xx
[0];
113 p
= ((q
< p
) << 16) + (q
>> 16);
120 p
+= (q
>> 16) >> 16;
126 /* --- Now reduce mod p --- *
128 * I'm using shifts and adds to do the multiply step here. This needs to
129 * be changed if @D@ ever becomes something other than 5.
133 # error "Change shift sequence!"
150 y
+= (q
>> 16) >> 16;
158 /* --- Now add on the constant --- */
169 /* --- @lcrand_range@ --- *
171 * Arguments: @uint32 *x@ = pointer to seed value (updated)
172 * @uint32 m@ = limit allowable
174 * Returns: A uniformly distributed pseudorandom integer in the interval
178 uint32
lcrand_range(uint32
*x
, uint32 m
)
181 uint32 r
= P
- P
% m
;
182 do xx
= lcrand(xx
); while (xx
>= r
);
187 /*----- Generic interface -------------------------------------------------*/
189 typedef struct gctx
{
194 static void gdestroy(grand
*r
)
200 static int gmisc(grand
*r
, unsigned op
, ...)
209 switch (va_arg(ap
, unsigned)) {
212 case GRAND_SEEDUINT32
:
222 g
->x
= va_arg(ap
, unsigned);
224 case GRAND_SEEDUINT32
:
225 g
->x
= va_arg(ap
, uint32
);
227 case GRAND_SEEDRAND
: {
228 grand
*rr
= va_arg(ap
, grand
*);
230 do x
= rr
->ops
->word(rr
); while (x
>= P
|| x
== LCRAND_FIXEDPT
);
242 static uint32
graw(grand
*r
)
249 static octet
gbyte(grand
*r
)
253 do x
= lcrand(x
); while (x
>= P256
);
255 return (x
/ (P256
/ 256));
258 static uint32
grange(grand
*r
, uint32 l
)
261 return (lcrand_range(&g
->x
, l
));
264 static const grand_ops gops
= {
268 graw
, gbyte
, grand_word
, grange
, grand_fill
271 /* --- @lcrand_create@ --- *
273 * Arguments: @uint32 x@ = initial seed
275 * Returns: Pointer to a generic generator.
277 * Use: Constructs a generic generator interface over a linear
278 * congruential generator.
281 grand
*lcrand_create(uint32 x
)
283 gctx
*g
= CREATE(gctx
);
289 /*----- Test rig ----------------------------------------------------------*/
293 #include <mLib/testrig.h>
295 static int verify(dstr
*v
)
297 uint32 x
= *(uint32
*)v
[0].buf
;
298 uint32 y
= *(uint32
*)v
[1].buf
;
299 uint32 z
= lcrand(x
);
303 "\n*** lcrand failed. lcrand(%lu) = %lu, expected %lu\n",
304 (unsigned long)x
, (unsigned long)z
, (unsigned long)y
);
310 static test_chunk tests
[] = {
311 { "lcrand", verify
, { &type_uint32
, &type_uint32
, 0 } },
315 int main(int argc
, char *argv
[])
317 test_run(argc
, argv
, tests
, SRCDIR
"/tests/lcrand");
323 /*----- That's all, folks -------------------------------------------------*/