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