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 | **************************************************************************/
|
---|
21 |
|
---|
22 | /*
|
---|
23 | * A random number generator called as a function by
|
---|
24 | * random (iseed) or irandm (iseed)
|
---|
25 | * The parameter should be a pointer to a 2-element long vector.
|
---|
26 | * The first function returns a double uniform in 0 .. 1.
|
---|
27 | * The second returns a long integer uniform in 0 .. 2**31-1
|
---|
28 | * Both update iseed[] in exactly the same way.
|
---|
29 | * iseed[] must be a 2-element integer vector.
|
---|
30 | * The initial value of the second element may be anything.
|
---|
31 | *
|
---|
32 | * The period of the random sequence is 2**32 * (2**32-1)
|
---|
33 | * The table mt[0:127] is defined by mt[i] = 69069 ** (128-i)
|
---|
34 | */
|
---|
35 | #include "random.h"
|
---|
36 |
|
---|
37 | #define MASK ((mg_s_long) 593970775)
|
---|
38 | /* or in hex, 23674657 */
|
---|
39 |
|
---|
40 | #define SCALE ((double) 1.0 / (1024.0 * 1024.0 * 1024.0 * 2.0))
|
---|
41 | /* i.e. 2 to power -31 */
|
---|
42 |
|
---|
43 | static mg_s_long mt [128] = {
|
---|
44 | 902906369,
|
---|
45 | 2030498053,
|
---|
46 | -473499623,
|
---|
47 | 1640834941,
|
---|
48 | 723406961,
|
---|
49 | 1993558325,
|
---|
50 | -257162999,
|
---|
51 | -1627724755,
|
---|
52 | 913952737,
|
---|
53 | 278845029,
|
---|
54 | 1327502073,
|
---|
55 | -1261253155,
|
---|
56 | 981676113,
|
---|
57 | -1785280363,
|
---|
58 | 1700077033,
|
---|
59 | 366908557,
|
---|
60 | -1514479167,
|
---|
61 | -682799163,
|
---|
62 | 141955545,
|
---|
63 | -830150595,
|
---|
64 | 317871153,
|
---|
65 | 1542036469,
|
---|
66 | -946413879,
|
---|
67 | -1950779155,
|
---|
68 | 985397153,
|
---|
69 | 626515237,
|
---|
70 | 530871481,
|
---|
71 | 783087261,
|
---|
72 | -1512358895,
|
---|
73 | 1031357269,
|
---|
74 | -2007710807,
|
---|
75 | -1652747955,
|
---|
76 | -1867214463,
|
---|
77 | 928251525,
|
---|
78 | 1243003801,
|
---|
79 | -2132510467,
|
---|
80 | 1874683889,
|
---|
81 | -717013323,
|
---|
82 | 218254473,
|
---|
83 | -1628774995,
|
---|
84 | -2064896159,
|
---|
85 | 69678053,
|
---|
86 | 281568889,
|
---|
87 | -2104168611,
|
---|
88 | -165128239,
|
---|
89 | 1536495125,
|
---|
90 | -39650967,
|
---|
91 | 546594317,
|
---|
92 | -725987007,
|
---|
93 | 1392966981,
|
---|
94 | 1044706649,
|
---|
95 | 687331773,
|
---|
96 | -2051306575,
|
---|
97 | 1544302965,
|
---|
98 | -758494647,
|
---|
99 | -1243934099,
|
---|
100 | -75073759,
|
---|
101 | 293132965,
|
---|
102 | -1935153095,
|
---|
103 | 118929437,
|
---|
104 | 807830417,
|
---|
105 | -1416222507,
|
---|
106 | -1550074071,
|
---|
107 | -84903219,
|
---|
108 | 1355292929,
|
---|
109 | -380482555,
|
---|
110 | -1818444007,
|
---|
111 | -204797315,
|
---|
112 | 170442609,
|
---|
113 | -1636797387,
|
---|
114 | 868931593,
|
---|
115 | -623503571,
|
---|
116 | 1711722209,
|
---|
117 | 381210981,
|
---|
118 | -161547783,
|
---|
119 | -272740131,
|
---|
120 | -1450066095,
|
---|
121 | 2116588437,
|
---|
122 | 1100682473,
|
---|
123 | 358442893,
|
---|
124 | -1529216831,
|
---|
125 | 2116152005,
|
---|
126 | -776333095,
|
---|
127 | 1265240893,
|
---|
128 | -482278607,
|
---|
129 | 1067190005,
|
---|
130 | 333444553,
|
---|
131 | 86502381,
|
---|
132 | 753481377,
|
---|
133 | 39000101,
|
---|
134 | 1779014585,
|
---|
135 | 219658653,
|
---|
136 | -920253679,
|
---|
137 | 2029538901,
|
---|
138 | 1207761577,
|
---|
139 | -1515772851,
|
---|
140 | -236195711,
|
---|
141 | 442620293,
|
---|
142 | 423166617,
|
---|
143 | -1763648515,
|
---|
144 | -398436623,
|
---|
145 | -1749358155,
|
---|
146 | -538598519,
|
---|
147 | -652439379,
|
---|
148 | 430550625,
|
---|
149 | -1481396507,
|
---|
150 | 2093206905,
|
---|
151 | -1934691747,
|
---|
152 | -962631983,
|
---|
153 | 1454463253,
|
---|
154 | -1877118871,
|
---|
155 | -291917555,
|
---|
156 | -1711673279,
|
---|
157 | 201201733,
|
---|
158 | -474645415,
|
---|
159 | -96764739,
|
---|
160 | -1587365199,
|
---|
161 | 1945705589,
|
---|
162 | 1303896393,
|
---|
163 | 1744831853,
|
---|
164 | 381957665,
|
---|
165 | 2135332261,
|
---|
166 | -55996615,
|
---|
167 | -1190135011,
|
---|
168 | 1790562961,
|
---|
169 | -1493191723,
|
---|
170 | 475559465,
|
---|
171 | 69069
|
---|
172 | };
|
---|
173 |
|
---|
174 | double
|
---|
175 | random (mg_s_long is [2])
|
---|
176 | {
|
---|
177 | mg_s_long it, leh, nit;
|
---|
178 |
|
---|
179 | it = is [0];
|
---|
180 | leh = is [1];
|
---|
181 | if (it <= 0)
|
---|
182 | it = (it + it) ^ MASK;
|
---|
183 | else
|
---|
184 | it = it + it;
|
---|
185 | nit = it - 1;
|
---|
186 | /* to ensure all-ones pattern omitted */
|
---|
187 | leh = leh * mt[nit & 127] + nit;
|
---|
188 | is [0] = it; is [1] = leh;
|
---|
189 | if (leh < 0) leh = ~leh;
|
---|
190 | return (SCALE * ((mg_s_long) (leh | 1)));
|
---|
191 | }
|
---|
192 |
|
---|
193 |
|
---|
194 |
|
---|
195 | mg_s_long
|
---|
196 | irandm (mg_s_long is [2])
|
---|
197 | {
|
---|
198 | mg_s_long it, leh, nit;
|
---|
199 |
|
---|
200 | it = is [0];
|
---|
201 | leh = is [1];
|
---|
202 | if (it <= 0)
|
---|
203 | it = (it + it) ^ MASK;
|
---|
204 | else
|
---|
205 | it = it + it;
|
---|
206 | nit = it - 1;
|
---|
207 | /* to ensure all-ones pattern omitted */
|
---|
208 | leh = leh * mt[nit & 127] + nit;
|
---|
209 | is [0] = it; is [1] = leh;
|
---|
210 | if (leh < 0) leh = ~leh;
|
---|
211 | return (leh);
|
---|
212 | }
|
---|