[3365] | 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 |
|
---|
| 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 |
|
---|
| 42 | static 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 |
|
---|
| 173 | double
|
---|
| 174 | random (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 |
|
---|
| 194 | long
|
---|
| 195 | irandm (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 | }
|
---|