rand.c

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