root/gsdl/trunk/trunk/mg/lib/random.c @ 16583

Revision 16583, 4.4 KB (checked in by davidb, 11 years ago)

Undoing change commited in r16582

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