1 | /**************************************************************************
|
---|
2 | *
|
---|
3 | * random.c -- pseudo random number generator
|
---|
4 | * Copyright (C) 1994 Chris Wallace ([email protected])
|
---|
5 | *
|
---|
6 | * This program is free software; you can redistribute it and/or modify
|
---|
7 | * it under the terms of the GNU General Public License as published by
|
---|
8 | * the Free Software Foundation; either version 2 of the License, or
|
---|
9 | * (at your option) any later version.
|
---|
10 | *
|
---|
11 | * This program is distributed in the hope that it will be useful,
|
---|
12 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
13 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
14 | * GNU General Public License for more details.
|
---|
15 | *
|
---|
16 | * You should have received a copy of the GNU General Public License
|
---|
17 | * along with this program; if not, write to the Free Software
|
---|
18 | * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
---|
19 | *
|
---|
20 | * $Id: random.c 439 1999-08-10 21:23:37Z sjboddie $
|
---|
21 | *
|
---|
22 | **************************************************************************/
|
---|
23 |
|
---|
24 | /*
|
---|
25 | $Log$
|
---|
26 | Revision 1.1 1999/08/10 21:16:57 sjboddie
|
---|
27 | renamed mg-1.3d directory mg
|
---|
28 |
|
---|
29 | Revision 1.1 1998/11/17 09:32:20 rjmcnab
|
---|
30 | *** empty log message ***
|
---|
31 |
|
---|
32 | */
|
---|
33 |
|
---|
34 | static char *RCSID = "$Id: random.c 439 1999-08-10 21:23:37Z sjboddie $";
|
---|
35 |
|
---|
36 | /*
|
---|
37 | * A random number generator called as a function by
|
---|
38 | * random (iseed) or irandm (iseed)
|
---|
39 | * The parameter should be a pointer to a 2-element long vector.
|
---|
40 | * The first function returns a double uniform in 0 .. 1.
|
---|
41 | * The second returns a long integer uniform in 0 .. 2**31-1
|
---|
42 | * Both update iseed[] in exactly the same way.
|
---|
43 | * iseed[] must be a 2-element integer vector.
|
---|
44 | * The initial value of the second element may be anything.
|
---|
45 | *
|
---|
46 | * The period of the random sequence is 2**32 * (2**32-1)
|
---|
47 | * The table mt[0:127] is defined by mt[i] = 69069 ** (128-i)
|
---|
48 | */
|
---|
49 |
|
---|
50 | #define MASK ((long) 593970775)
|
---|
51 | /* or in hex, 23674657 */
|
---|
52 |
|
---|
53 | #define SCALE ((double) 1.0 / (1024.0 * 1024.0 * 1024.0 * 2.0))
|
---|
54 | /* i.e. 2 to power -31 */
|
---|
55 |
|
---|
56 | static long mt [128] = {
|
---|
57 | 902906369,
|
---|
58 | 2030498053,
|
---|
59 | -473499623,
|
---|
60 | 1640834941,
|
---|
61 | 723406961,
|
---|
62 | 1993558325,
|
---|
63 | -257162999,
|
---|
64 | -1627724755,
|
---|
65 | 913952737,
|
---|
66 | 278845029,
|
---|
67 | 1327502073,
|
---|
68 | -1261253155,
|
---|
69 | 981676113,
|
---|
70 | -1785280363,
|
---|
71 | 1700077033,
|
---|
72 | 366908557,
|
---|
73 | -1514479167,
|
---|
74 | -682799163,
|
---|
75 | 141955545,
|
---|
76 | -830150595,
|
---|
77 | 317871153,
|
---|
78 | 1542036469,
|
---|
79 | -946413879,
|
---|
80 | -1950779155,
|
---|
81 | 985397153,
|
---|
82 | 626515237,
|
---|
83 | 530871481,
|
---|
84 | 783087261,
|
---|
85 | -1512358895,
|
---|
86 | 1031357269,
|
---|
87 | -2007710807,
|
---|
88 | -1652747955,
|
---|
89 | -1867214463,
|
---|
90 | 928251525,
|
---|
91 | 1243003801,
|
---|
92 | -2132510467,
|
---|
93 | 1874683889,
|
---|
94 | -717013323,
|
---|
95 | 218254473,
|
---|
96 | -1628774995,
|
---|
97 | -2064896159,
|
---|
98 | 69678053,
|
---|
99 | 281568889,
|
---|
100 | -2104168611,
|
---|
101 | -165128239,
|
---|
102 | 1536495125,
|
---|
103 | -39650967,
|
---|
104 | 546594317,
|
---|
105 | -725987007,
|
---|
106 | 1392966981,
|
---|
107 | 1044706649,
|
---|
108 | 687331773,
|
---|
109 | -2051306575,
|
---|
110 | 1544302965,
|
---|
111 | -758494647,
|
---|
112 | -1243934099,
|
---|
113 | -75073759,
|
---|
114 | 293132965,
|
---|
115 | -1935153095,
|
---|
116 | 118929437,
|
---|
117 | 807830417,
|
---|
118 | -1416222507,
|
---|
119 | -1550074071,
|
---|
120 | -84903219,
|
---|
121 | 1355292929,
|
---|
122 | -380482555,
|
---|
123 | -1818444007,
|
---|
124 | -204797315,
|
---|
125 | 170442609,
|
---|
126 | -1636797387,
|
---|
127 | 868931593,
|
---|
128 | -623503571,
|
---|
129 | 1711722209,
|
---|
130 | 381210981,
|
---|
131 | -161547783,
|
---|
132 | -272740131,
|
---|
133 | -1450066095,
|
---|
134 | 2116588437,
|
---|
135 | 1100682473,
|
---|
136 | 358442893,
|
---|
137 | -1529216831,
|
---|
138 | 2116152005,
|
---|
139 | -776333095,
|
---|
140 | 1265240893,
|
---|
141 | -482278607,
|
---|
142 | 1067190005,
|
---|
143 | 333444553,
|
---|
144 | 86502381,
|
---|
145 | 753481377,
|
---|
146 | 39000101,
|
---|
147 | 1779014585,
|
---|
148 | 219658653,
|
---|
149 | -920253679,
|
---|
150 | 2029538901,
|
---|
151 | 1207761577,
|
---|
152 | -1515772851,
|
---|
153 | -236195711,
|
---|
154 | 442620293,
|
---|
155 | 423166617,
|
---|
156 | -1763648515,
|
---|
157 | -398436623,
|
---|
158 | -1749358155,
|
---|
159 | -538598519,
|
---|
160 | -652439379,
|
---|
161 | 430550625,
|
---|
162 | -1481396507,
|
---|
163 | 2093206905,
|
---|
164 | -1934691747,
|
---|
165 | -962631983,
|
---|
166 | 1454463253,
|
---|
167 | -1877118871,
|
---|
168 | -291917555,
|
---|
169 | -1711673279,
|
---|
170 | 201201733,
|
---|
171 | -474645415,
|
---|
172 | -96764739,
|
---|
173 | -1587365199,
|
---|
174 | 1945705589,
|
---|
175 | 1303896393,
|
---|
176 | 1744831853,
|
---|
177 | 381957665,
|
---|
178 | 2135332261,
|
---|
179 | -55996615,
|
---|
180 | -1190135011,
|
---|
181 | 1790562961,
|
---|
182 | -1493191723,
|
---|
183 | 475559465,
|
---|
184 | 69069
|
---|
185 | };
|
---|
186 |
|
---|
187 | double
|
---|
188 | random (long is [2])
|
---|
189 | {
|
---|
190 | long it, leh, nit;
|
---|
191 |
|
---|
192 | it = is [0];
|
---|
193 | leh = is [1];
|
---|
194 | if (it <= 0)
|
---|
195 | it = (it + it) ^ MASK;
|
---|
196 | else
|
---|
197 | it = it + it;
|
---|
198 | nit = it - 1;
|
---|
199 | /* to ensure all-ones pattern omitted */
|
---|
200 | leh = leh * mt[nit & 127] + nit;
|
---|
201 | is [0] = it; is [1] = leh;
|
---|
202 | if (leh < 0) leh = ~leh;
|
---|
203 | return (SCALE * ((long) (leh | 1)));
|
---|
204 | }
|
---|
205 |
|
---|
206 |
|
---|
207 |
|
---|
208 | long
|
---|
209 | irandm (long is [2])
|
---|
210 | {
|
---|
211 | long it, leh, nit;
|
---|
212 |
|
---|
213 | it = is [0];
|
---|
214 | leh = is [1];
|
---|
215 | if (it <= 0)
|
---|
216 | it = (it + it) ^ MASK;
|
---|
217 | else
|
---|
218 | it = it + it;
|
---|
219 | nit = it - 1;
|
---|
220 | /* to ensure all-ones pattern omitted */
|
---|
221 | leh = leh * mt[nit & 127] + nit;
|
---|
222 | is [0] = it; is [1] = leh;
|
---|
223 | if (leh < 0) leh = ~leh;
|
---|
224 | return (leh);
|
---|
225 | }
|
---|