5 * Simple linear congruential generator
7 * (c) 1999 Straylight/Edgeware
8 * (c) 2000 Mark Wooding
11 /*----- Licensing notice --------------------------------------------------*
13 * Copyright (c) 2000 Mark Wooding
14 * All rights reserved.
16 * Redistribution and use in source and binary forms, with or without
17 * modification, are permitted provided that the following conditions are
20 * 1. Redistributions of source code must retain the above copyright
21 * notice, this list of conditions and the following disclaimer.
23 * 2, Redistributions in binary form must reproduce the above copyright
24 * notice, this list of conditions and the following disclaimer in the
25 * documentation and/or other materials provided with the distribution.
27 * 3. The name of the authors may not be used to endorse or promote
28 * products derived from this software without specific prior written
31 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESS OR IMPLIED
32 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
33 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
34 * NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
35 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
36 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
37 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
38 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
39 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
40 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
41 * POSSIBILITY OF SUCH DAMAGE.
43 * Instead of accepting the above terms, you may redistribute and/or modify
44 * this software under the terms of either the GNU General Public License,
45 * or the GNU Library General Public License, published by the Free
46 * Software Foundation; either version 2 of the License, or (at your
47 * option) any later version.
50 /*----- Revision history --------------------------------------------------*
53 * Revision 1.2 2000/07/02 15:21:20 mdw
56 * Revision 1.1 2000/05/21 11:28:30 mdw
59 * --- Past lives (Catacomb) --- *
61 * Revision 1.2 1999/12/13 15:34:01 mdw
62 * Add support for seeding from a generic pseudorandom source.
64 * Revision 1.1 1999/12/10 23:15:27 mdw
65 * Noncryptographic random number generator.
69 /*----- Header files ------------------------------------------------------*/
78 /*----- Magic numbers -----------------------------------------------------*/
80 /* --- The generator parameters --- */
82 #define P LCRAND_P /* Modulus */
83 #define A LCRAND_A /* Multiplier (primitive mod @p@) */
84 #define C LCRAND_C /* Additive constant */
86 /* --- Precomputed values for modular reduction --- */
88 #define D 5 /* %$p = 2^{32} - d$% */
90 /* --- Other useful bits --- */
92 #define P256 4294967040u /* Highest multiple of 256 < %$p$% */
94 /*----- Main code ---------------------------------------------------------*/
98 * Arguments: @uint32 x@ = seed value
100 * Returns: New state of the generator.
102 * Use: Steps the generator. Returns %$ax + c \bmod p$%.
105 uint32 lcrand(uint32 x)
110 /* --- Unpack things into the arrays --- */
112 a[0] = U16(A); a[1] = U16(A >> 16);
113 xx[0] = U16(x); xx[1] = U16(x >> 16);
115 /* --- Multiply everything together --- *
117 * This is plain old long multiplication, although it looks a bit strange.
118 * I set up the top and bottom partial products directly where they're
119 * supposed to be. The cross terms I add together, with the low 16 bits in
120 * @q@ and the high 32 bits in @p@. These I then add into the product.
126 yy[0] = a[0] * xx[0];
127 yy[1] = a[1] * xx[1];
130 q = p + a[1] * xx[0];
131 p = ((q < p) << 16) + (q >> 16);
138 p += (q >> 16) >> 16;
144 /* --- Now reduce mod p --- *
146 * I'm using shifts and adds to do the multiply step here. This needs to
147 * be changed if @D@ ever becomes something other than 5.
151 # error "Change shift sequence!"
168 y += (q >> 16) >> 16;
176 /* --- Now add on the constant --- */
187 /* --- @lcrand_range@ --- *
189 * Arguments: @uint32 *x@ = pointer to seed value (updated)
190 * @uint32 m@ = limit allowable
192 * Returns: A uniformly distributed pseudorandom integer in the interval
196 uint32 lcrand_range(uint32 *x, uint32 m)
199 uint32 r = P - P % m;
200 do xx = lcrand(xx); while (xx >= r);
202 return (xx / (r / m));
205 /*----- That's all, folks -------------------------------------------------*/