chiark / gitweb /
Version bump.
[storin] / lcrand.c
CommitLineData
e6e0e332
MW
1/* -*-c-*-
2 *
6b2d9d76 3 * $Id: lcrand.c,v 1.2 2000/07/02 15:21:20 mdw Exp $
e6e0e332
MW
4 *
5 * Simple linear congruential generator
6 *
7 * (c) 1999 Straylight/Edgeware
8 * (c) 2000 Mark Wooding
9 */
10
11/*----- Licensing notice --------------------------------------------------*
12 *
13 * Copyright (c) 2000 Mark Wooding
14 * All rights reserved.
15 *
16 * Redistribution and use in source and binary forms, with or without
17 * modification, are permitted provided that the following conditions are
18 * met:
19 *
20 * 1. Redistributions of source code must retain the above copyright
21 * notice, this list of conditions and the following disclaimer.
22 *
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.
26 *
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
29 * permission.
30 *
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
6b2d9d76 34 * NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
e6e0e332
MW
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.
42 *
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.
48 */
49
50/*----- Revision history --------------------------------------------------*
51 *
52 * $Log: lcrand.c,v $
6b2d9d76
MW
53 * Revision 1.2 2000/07/02 15:21:20 mdw
54 * Fix licence text.
55 *
e6e0e332
MW
56 * Revision 1.1 2000/05/21 11:28:30 mdw
57 * Initial check-in.
58 *
59 * --- Past lives (Catacomb) --- *
60 *
61 * Revision 1.2 1999/12/13 15:34:01 mdw
62 * Add support for seeding from a generic pseudorandom source.
63 *
64 * Revision 1.1 1999/12/10 23:15:27 mdw
65 * Noncryptographic random number generator.
66 *
67 */
68
69/*----- Header files ------------------------------------------------------*/
70
71#include <stdio.h>
72#include <stdlib.h>
73#include <string.h>
74
75#include "bits.h"
76#include "lcrand.h"
77
78/*----- Magic numbers -----------------------------------------------------*/
79
80/* --- The generator parameters --- */
81
82#define P LCRAND_P /* Modulus */
83#define A LCRAND_A /* Multiplier (primitive mod @p@) */
84#define C LCRAND_C /* Additive constant */
85
86/* --- Precomputed values for modular reduction --- */
87
88#define D 5 /* %$p = 2^{32} - d$% */
89
90/* --- Other useful bits --- */
91
92#define P256 4294967040u /* Highest multiple of 256 < %$p$% */
93
94/*----- Main code ---------------------------------------------------------*/
95
96/* --- @lcrand@ --- *
97 *
98 * Arguments: @uint32 x@ = seed value
99 *
100 * Returns: New state of the generator.
101 *
102 * Use: Steps the generator. Returns %$ax + c \bmod p$%.
103 */
104
105uint32 lcrand(uint32 x)
106{
107 uint32 a[2], xx[2];
108 uint32 yy[2];
109
110 /* --- Unpack things into the arrays --- */
111
112 a[0] = U16(A); a[1] = U16(A >> 16);
113 xx[0] = U16(x); xx[1] = U16(x >> 16);
114
115 /* --- Multiply everything together --- *
116 *
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.
121 */
122
123 {
124 uint32 p, q;
125
126 yy[0] = a[0] * xx[0];
127 yy[1] = a[1] * xx[1];
128
129 p = a[0] * xx[1];
130 q = p + a[1] * xx[0];
131 p = ((q < p) << 16) + (q >> 16);
132 q = U16(q) << 16;
133
134 q += yy[0];
135 if (q < yy[0])
136 p++;
137 else
138 p += (q >> 16) >> 16;
139 yy[0] = q;
140
141 yy[1] += p;
142 }
143
144 /* --- Now reduce mod p --- *
145 *
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.
148 */
149
150#if D != 5
151# error "Change shift sequence!"
152#endif
153
154 {
155 uint32 q;
156
157 q = yy[1];
158 x = yy[0];
159
160 while (q) {
161 uint32 y, z;
162 y = q >> 30;
163 z = q << 2;
164 z += q;
165 if (z < q)
166 y++;
167 else
168 y += (q >> 16) >> 16;
169 q = y;
170 x += z;
171 if (x < z || x > P)
172 x -= P;
173 }
174 }
175
176 /* --- Now add on the constant --- */
177
178 x += C;
179 if (x < C || x >= P)
180 x -= P;
181
182 /* --- Done --- */
183
184 return (x);
185}
186
187/* --- @lcrand_range@ --- *
188 *
189 * Arguments: @uint32 *x@ = pointer to seed value (updated)
190 * @uint32 m@ = limit allowable
191 *
192 * Returns: A uniformly distributed pseudorandom integer in the interval
193 * %$[0, m)$%.
194 */
195
196uint32 lcrand_range(uint32 *x, uint32 m)
197{
198 uint32 xx = *x;
199 uint32 r = P - P % m;
200 do xx = lcrand(xx); while (xx >= r);
201 *x = xx;
202 return (xx / (r / m));
203}
204
205/*----- That's all, folks -------------------------------------------------*/