math/mpreduce.h: Missing include files.
[u/mdw/catacomb] / rand / fibrand.c
1 /* -*-c-*-
2 *
3 * Fibonacci generator
4 *
5 * (c) 1999 Straylight/Edgeware
6 */
7
8 /*----- Licensing notice --------------------------------------------------*
9 *
10 * This file is part of Catacomb.
11 *
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.
16 *
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.
21 *
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,
25 * MA 02111-1307, USA.
26 */
27
28 /*----- Header files ------------------------------------------------------*/
29
30 #include <stdarg.h>
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <string.h>
34
35 #include <mLib/bits.h>
36 #include <mLib/sub.h>
37
38 #include "fibrand.h"
39 #include "grand.h"
40 #include "lcrand.h"
41
42 /*----- Main code ---------------------------------------------------------*/
43
44 /* --- @fibrand_step@ --- *
45 *
46 * Arguments: @fibrand *f@ = pointer to Fibonacci generator context
47 *
48 * Returns: Next output from generator.
49 *
50 * Use: Steps the generator. Returns
51 * %$x_{i - 24} + x_{i - 55} \bmod 2^{32}$%.
52 */
53
54 uint32 fibrand_step(fibrand *f)
55 {
56 unsigned i = f->i;
57 unsigned j = i + (FIB_SZ - FIB_TAP);
58 uint32 x;
59 if (j >= FIB_SZ)
60 j -= FIB_SZ;
61 x = f->x[i] = U32(f->x[i] + f->x[j]);
62 i++;
63 if (i >= FIB_SZ)
64 i = 0;
65 f->i = i;
66 return (x);
67 }
68
69 /* --- @fibrand_seed@ --- *
70 *
71 * Arguments: @fibrand *f@ = pointer to Fibonacci generator context
72 * @grand *r@ = random number generator to extract words from
73 *
74 * Returns: ---
75 *
76 * Use: Initializes a Fibonacci generator using word outputs from the
77 * given random number source @r@.
78 */
79
80 void fibrand_seed(fibrand *f, grand *r)
81 {
82 int i;
83 unsigned p = 0;
84
85 for (i = 0; i < FIB_SZ; i++)
86 p |= f->x[i] = r->ops->word(r);
87 if (!(p & 1)) {
88 i = r->ops->range(r, FIB_SZ);
89 f->x[i] |= 1;
90 }
91 f->i = 0;
92 }
93
94 /* --- @fibrand_lcseed@ --- *
95 *
96 * Arguments: @fibrand *f@ = pointer to Fibonacci generator context
97 * @uint32 seed@ = seed value
98 *
99 * Returns: ---
100 *
101 * Use: Initializes a Fibonacci generator using outputs from the
102 * @lcrand@ generator seeded from @seed@. This is faster than
103 * using a generic @lcrand@-based generator and @fibrand_rseed@
104 * because it uses raw outputs rather than uniformly distributed
105 * 32-bit words.
106 */
107
108 void fibrand_lcseed(fibrand *f, uint32 seed)
109 {
110 int i;
111 unsigned p = 0;
112
113 for (i = 0; i < FIB_SZ; i++)
114 p |= f->x[i] = seed = lcrand(seed);
115 if (!(p & 1)) {
116 i = lcrand_range(&seed, FIB_SZ);
117 f->x[i] |= 1;
118 }
119 f->i = 0;
120 }
121
122 /* --- @fibrand_range@ --- *
123 *
124 * Arguments: @fibrand *f@ = pointer to Fibonacci generator context
125 * @uint32 m@ = limit
126 *
127 * Returns: A uniformly distributed pseudorandom integer in the interval
128 * %$[0, m)$%.
129 */
130
131 uint32 fibrand_range(fibrand *f, uint32 m)
132 {
133 uint32 r = 0xffffffff - (0xffffffff % m);
134 uint32 x;
135
136 /* --- Now generate numbers until a good one comes along --- */
137
138 do x = fibrand_step(f); while (x >= r);
139 return (x % m);
140 }
141
142 /*----- Generic interface -------------------------------------------------*/
143
144 typedef struct gctx {
145 grand r;
146 fibrand f;
147 } gctx;
148
149 static void gdestroy(grand *r)
150 {
151 gctx *g = (gctx *)r;
152 DESTROY(g);
153 }
154
155 static int gmisc(grand *r, unsigned op, ...)
156 {
157 gctx *g = (gctx *)r;
158 va_list ap;
159 int rc = 0;
160 va_start(ap, op);
161
162 switch (op) {
163 case GRAND_CHECK:
164 switch (va_arg(ap, unsigned)) {
165 case GRAND_CHECK:
166 case GRAND_SEEDINT:
167 case GRAND_SEEDUINT32:
168 case GRAND_SEEDRAND:
169 rc = 1;
170 break;
171 default:
172 rc = 0;
173 break;
174 }
175 break;
176 case GRAND_SEEDINT:
177 fibrand_lcseed(&g->f, va_arg(ap, unsigned));
178 break;
179 case GRAND_SEEDUINT32:
180 fibrand_lcseed(&g->f, va_arg(ap, uint32));
181 break;
182 case GRAND_SEEDRAND:
183 fibrand_seed(&g->f, va_arg(ap, grand *));
184 break;
185 default:
186 GRAND_BADOP;
187 break;
188 }
189
190 va_end(ap);
191 return (rc);
192 }
193
194 static octet gbyte(grand *r)
195 {
196 gctx *g = (gctx *)r;
197 return (U8(fibrand_step(&g->f)));
198 }
199
200 static uint32 gword(grand *r)
201 {
202 gctx *g = (gctx *)r;
203 return (fibrand_step(&g->f));
204 }
205
206 static uint32 grange(grand *r, uint32 l)
207 {
208 gctx *g = (gctx *)r;
209 return (fibrand_range(&g->f, l));
210 }
211
212 static void gfill(grand *r, void *p, size_t sz)
213 {
214 gctx *g = (gctx *)r;
215 octet *q = p;
216 while (sz) {
217 *q++ = U8(fibrand_step(&g->f));
218 sz--;
219 }
220 }
221
222 static const grand_ops gops = {
223 "fibrand",
224 0, 0,
225 gmisc, gdestroy,
226 gword, gbyte, gword, grange, gfill
227 };
228
229 /* --- @fibrand_create@ --- *
230 *
231 * Arguments: @uint32 seed@ = initial seed
232 *
233 * Returns: Pointer to a generic generator.
234 *
235 * Use: Constructs a generic generator interface over a Fibonacci
236 * generator. The generator is seeded using @fibrand_lcseed@.
237 */
238
239 grand *fibrand_create(uint32 seed)
240 {
241 gctx *g = CREATE(gctx);
242 g->r.ops = &gops;
243 fibrand_lcseed(&g->f, seed);
244 return (&g->r);
245 }
246
247 /*----- That's all, folks -------------------------------------------------*/