src/random.c

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