1 /* $OpenBSD: drand.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */
2
3 /*
4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 *
6 * Permission to use, copy, modify, and distribute this software for any
7 * purpose with or without fee is hereby granted, provided that the above
8 * copyright notice and this permission notice appear in all copies.
9 *
10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17 */
18
19 /* drand.c
20 *
21 * Pseudorandom number generator
22 *
23 *
24 *
25 * SYNOPSIS:
26 *
27 * double y, drand();
28 *
29 * drand( &y );
30 *
31 *
32 *
33 * DESCRIPTION:
34 *
35 * Yields a random number 1.0 <= y < 2.0.
36 *
37 * The three-generator congruential algorithm by Brian
38 * Wichmann and David Hill (BYTE magazine, March, 1987,
39 * pp 127-8) is used. The period, given by them, is
40 * 6953607871644.
41 *
42 * Versions invoked by the different arithmetic compile
43 * time options DEC, IBMPC, and MIEEE, produce
44 * approximately the same sequences, differing only in the
45 * least significant bits of the numbers. The UNK option
46 * implements the algorithm as recommended in the BYTE
47 * article. It may be used on all computers. However,
48 * the low order bits of a double precision number may
49 * not be adequately random, and may vary due to arithmetic
50 * implementation details on different computers.
51 *
52 * The other compile options generate an additional random
53 * integer that overwrites the low order bits of the double
54 * precision number. This reduces the period by a factor of
55 * two but tends to overcome the problems mentioned.
56 *
57 */
58
59 #include "mconf.h"
60
61
62 /* Three-generator random number algorithm
63 * of Brian Wichmann and David Hill
64 * BYTE magazine, March, 1987 pp 127-8
65 *
66 * The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
67 */
68
69 static int sx = 1;
70 static int sy = 10000;
71 static int sz = 3000;
72
73 static union {
74 double d;
75 unsigned short s[4];
76 } unkans;
77
78 /* This function implements the three
79 * congruential generators.
80 */
81
ranwh()82 static int ranwh()
83 {
84 int r, s;
85
86 /* sx = sx * 171 mod 30269 */
87 r = sx/177;
88 s = sx - 177 * r;
89 sx = 171 * s - 2 * r;
90 if( sx < 0 )
91 sx += 30269;
92
93
94 /* sy = sy * 172 mod 30307 */
95 r = sy/176;
96 s = sy - 176 * r;
97 sy = 172 * s - 35 * r;
98 if( sy < 0 )
99 sy += 30307;
100
101 /* sz = 170 * sz mod 30323 */
102 r = sz/178;
103 s = sz - 178 * r;
104 sz = 170 * s - 63 * r;
105 if( sz < 0 )
106 sz += 30323;
107 /* The results are in static sx, sy, sz. */
108 return 0;
109 }
110
111 /* drand.c
112 *
113 * Random double precision floating point number between 1 and 2.
114 *
115 * C callable:
116 * drand( &x );
117 */
118
drand(a)119 int drand( a )
120 double *a;
121 {
122 unsigned short r;
123 #ifdef DEC
124 unsigned short s, t;
125 #endif
126
127 /* This algorithm of Wichmann and Hill computes a floating point
128 * result:
129 */
130 ranwh();
131 unkans.d = sx/30269.0 + sy/30307.0 + sz/30323.0;
132 r = unkans.d;
133 unkans.d -= r;
134 unkans.d += 1.0;
135
136 /* if UNK option, do nothing further.
137 * Otherwise, make a random 16 bit integer
138 * to overwrite the least significant word
139 * of unkans.
140 */
141 #ifdef UNK
142 /* do nothing */
143 #else
144 ranwh();
145 r = sx * sy + sz;
146 #endif
147
148 #ifdef DEC
149 /* To make the numbers as similar as possible
150 * in all arithmetics, the random integer has
151 * to be inserted 3 bits higher up in a DEC number.
152 * An alternative would be put it 3 bits lower down
153 * in all the other number types.
154 */
155 s = unkans.s[2];
156 t = s & 07; /* save these bits to put in at the bottom */
157 s &= 0177770;
158 s |= (r >> 13) & 07;
159 unkans.s[2] = s;
160 t |= r << 3;
161 unkans.s[3] = t;
162 #endif
163
164 #ifdef IBMPC
165 unkans.s[0] = r;
166 #endif
167
168 #ifdef MIEEE
169 unkans.s[3] = r;
170 #endif
171
172 *a = unkans.d;
173 return 0;
174 }
175