src/mtrand.c

   1 /*
   2  * netsniff-ng - the packet sniffing beast
   3  * By Daniel Borkmann <daniel@netsniff-ng.org>
   4  * Copyright 2009, 2010 Daniel Borkmann.
   5  * Subject to the GPL, version 2.
   6  */
   7
   8 /*
   9  * Copyright (C) 1997-2004, Makoto Matsumoto, Takuji Nishimura, and
  10  * Eric Landry; All rights reserved.
  11  * Daniel Borkmann: Refactored, added initialization functions.
  12  *
  13  * Redistribution and use in source and binary forms, with or without
  14  * modification, are permitted provided that the following conditions
  15  * are met:
  16  *
  17  *   1. Redistributions of source code must retain the above copyright
  18  *      notice, this list of conditions and the following disclaimer.
  19  *
  20  *   2. Redistributions in binary form must reproduce the above copyright
  21  *      notice, this list of conditions and the following disclaimer
  22  *      in the documentation and/or other materials provided with the
  23  *      distribution.
  24  *
  25  *   3. The names of its contributors may not be used to endorse or
  26  *      promote products derived from this software without specific
  27  *      prior written permission.
  28  *
  29  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  30  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  31  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  32  * A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT
  33  * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  34  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  35  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  36  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  37  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  38  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  39  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  40  *
  41  * Any feedback is very welcome.
  42  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
  43  * email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
  44  *
  45  * Reference: M. Matsumoto and T. Nishimura, "Mersenne Twister:
  46  * A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number
  47  * Generator", ACM Transactions on Modeling and Computer Simulation,
  48  * Vol. 8, No. 1, January 1998, pp 3--30.
  49  */
  50
  51 #include <stdio.h>
  52 #include <stdlib.h>
  53 #include <time.h>
  54 #include <fcntl.h>
  55 #include <unistd.h>
  56
  57 #include "mtrand.h"
  58 #include "xio.h"
  59
  60 #define N           624
  61 #define M           397
  62 #define LEN_INIT    256
  63
  64 #define MATRIX_A    0x9908b0dfUL
  65 #define UPPER_MASK  0x80000000UL
  66 #define LOWER_MASK  0x7fffffffUL
  67
  68 static unsigned long x[N];
  69 static unsigned long *p0, *p1, *pm;
  70
  71 /*
  72  *  Initialize with a seed.
  73  *
  74  *  See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
  75  *  In the previous versions, MSBs of the seed affect only MSBs of
  76  *  the state.
  77  *  2002-01-09 modified by Makoto Matsumoto
  78  */
  79 void mt_init_by_seed_rand(unsigned long s)
  80 {
  81         int i;
  82
  83         x[0] = s & 0xffffffffUL;
  84
  85         for (i = 1; i < N; ++i) {
  86                 x[i] = (1812433253UL * (x[i - 1] ^ (x[i - 1] >> 30)) + i) &
  87                        0xffffffffUL;
  88         }
  89
  90         p0 = x;
  91         p1 = x + 1;
  92         pm = x + M;
  93 }
  94
  95 /*
  96  *  Initialize with time as seed.
  97  */
  98 void mt_init_by_seed_time(void)
  99 {
 100         int i;
 101
 102         x[0] = ((unsigned long) time(NULL)) & 0xffffffffUL;
 103
 104         for (i = 1; i < N; ++i) {
 105                 x[i] = (1812433253UL * (x[i - 1] ^ (x[i - 1] >> 30)) + i) &
 106                        0xffffffffUL;
 107         }
 108
 109         p0 = x;
 110         p1 = x + 1;
 111         pm = x + M;
 112 }
 113
 114 /*
 115  * Initialize by an array with array-length.
 116  */
 117 void mt_init_by_seed_array(unsigned long key[], int len)
 118 {
 119         int i, j, k;
 120
 121         mt_init_by_seed_rand(19650218UL);
 122
 123         i = 1;
 124         j = 0;
 125
 126         for (k = (N > len ? N : len); k; --k) {
 127                 /* Non linear */
 128                 x[i] = ((x[i] ^ ((x[i - 1] ^ (x[i - 1] >> 30)) *
 129                        1664525UL)) + key[j] + j) & 0xffffffffUL;
 130
 131                 if (++i >= N) {
 132                         x[0] = x[N - 1];
 133                         i = 1;
 134                 }
 135
 136                 if (++j >= len)
 137                         j = 0;
 138         }
 139
 140         for (k = N - 1; k; --k) {
 141                 /* Non linear */
 142                 x[i] = ((x[i] ^ ((x[i - 1] ^ (x[i - 1] >> 30)) *
 143                        1566083941UL)) - i) & 0xffffffffUL;
 144
 145                 if (++i >= N) {
 146                         x[0] = x[N - 1];
 147                         i = 1;
 148                 }
 149         }
 150
 151         x[0] = 0x80000000UL;
 152 }
 153
 154 /*
 155  * Initialize by an random array.
 156  */
 157 void mt_init_by_seed_rand_array(void)
 158 {
 159         int i;
 160         unsigned long k[LEN_INIT];
 161
 162         srand((unsigned int) time(NULL));
 163         for (i = 0; i < LEN_INIT; i++)
 164                 k[i] = rand();
 165
 166         mt_init_by_seed_array(k, LEN_INIT);
 167 }
 168
 169 /*
 170  * Initialize by an random array read from /dev/random
 171  */
 172 void mt_init_by_random_device(void)
 173 {
 174         int fd;
 175         unsigned long k[LEN_INIT];
 176
 177         fd = open_or_die("/dev/random", O_RDONLY);
 178         read_or_die(fd, k, sizeof(unsigned long) * LEN_INIT);
 179         close(fd);
 180
 181         mt_init_by_seed_array(k, LEN_INIT);
 182 }
 183
 184 /*
 185  * Generates a random number on the interval [0,0xffffffff]
 186  */
 187 unsigned long mt_rand_int32(void)
 188 {
 189         unsigned long y;
 190
 191         /* Default seed */
 192         if (p0 == NULL)
 193                 mt_init_by_seed_rand(5489UL);
 194
 195         /* Twisted feedback */
 196         y = *p0 = *pm++ ^ (((*p0 & UPPER_MASK) | (*p1 & LOWER_MASK)) >> 1) ^
 197                   (-(*p1 & 1) & MATRIX_A);
 198
 199         p0 = p1++;
 200
 201         if (pm == x + N)
 202                 pm = x;
 203         if (p1 == x + N)
 204                 p1 = x;
 205
 206         /* Temper */
 207         y ^= y >> 11;
 208         y ^= y << 7 & 0x9d2c5680UL;
 209         y ^= y << 15 & 0xefc60000UL;
 210         y ^= y >> 18;
 211
 212         return y;
 213 }
 214
 215 /*
 216  * Generates a random number on the interval [0,0x7fffffff]
 217  */
 218 long mt_rand_int31(void)
 219 {
 220         return (long) mt_rand_int32() >> 1;
 221 }
 222
 223 /*
 224  * Generates a random number on the real interval [0,1]
 225  */
 226 double mt_rand_real1(void)
 227 {
 228         return mt_rand_int32() * (1.0 / 4294967295.0);
 229         /* Divided by 2^32-1 */
 230 }
 231
 232 /*
 233  * Generates a random number on the real interval [0,1)
 234  */
 235 double mt_rand_real2(void)
 236 {
 237         return mt_rand_int32() * (1.0 / 4294967296.0);
 238         /* Divided by 2^32 */
 239 }
 240
 241 /*
 242  * Generates a random number on the real interval (0,1)
 243  */
 244 double mt_rand_real3(void)
 245 {
 246         return (((double) mt_rand_int32()) + 0.5) * (1.0 / 4294967296.0);
 247         /* Divided by 2^32 */
 248 }
 249
 250 /*
 251  * Generates a 53-bit random number on the real interval [0,1)
 252  */
 253 double mt_rand_res53(void)
 254 {
 255         unsigned long a = mt_rand_int32() >> 5, b = mt_rand_int32() >> 6;
 256         return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
 257 }
 258