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  * Copyright (C) 1997-2004, Makoto Matsumoto, Takuji Nishimura, and
   6  * Eric Landry; All rights reserved. (3-clause BSD license)
   7  * Daniel Borkmann: Refactored, added initialization functions.
   8  * Subject to the GPL, version 2.
   9  * Reference: M. Matsumoto and T. Nishimura, "Mersenne Twister:
  10  * A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number
  11  * Generator", ACM Transactions on Modeling and Computer Simulation,
  12  * Vol. 8, No. 1, January 1998, pp 3--30.
  13  */
  14
  15 #include <stdio.h>
  16 #include <stdlib.h>
  17 #include <time.h>
  18 #include <fcntl.h>
  19 #include <unistd.h>
  20
  21 #include "mtrand.h"
  22 #include "xio.h"
  23
  24 #define N           624
  25 #define M           397
  26 #define LEN_INIT    256
  27
  28 #define MATRIX_A    0x9908b0dfUL
  29 #define UPPER_MASK  0x80000000UL
  30 #define LOWER_MASK  0x7fffffffUL
  31
  32 static unsigned long x[N];
  33 static unsigned long *p0, *p1, *pm;
  34
  35 void mt_init_by_seed_rand(unsigned long s)
  36 {
  37         int i;
  38         x[0] = s & 0xffffffffUL;
  39         for (i = 1; i < N; ++i) {
  40                 x[i] = (1812433253UL * (x[i - 1] ^ (x[i - 1] >> 30)) + i) &
  41                        0xffffffffUL;
  42         }
  43         p0 = x;
  44         p1 = x + 1;
  45         pm = x + M;
  46 }
  47
  48 void mt_init_by_seed_time(void)
  49 {
  50         int i;
  51         x[0] = ((unsigned long) time(NULL)) & 0xffffffffUL;
  52         for (i = 1; i < N; ++i) {
  53                 x[i] = (1812433253UL * (x[i - 1] ^ (x[i - 1] >> 30)) + i) &
  54                        0xffffffffUL;
  55         }
  56         p0 = x;
  57         p1 = x + 1;
  58         pm = x + M;
  59 }
  60
  61 void mt_init_by_seed_array(unsigned long key[], int len)
  62 {
  63         int i, j, k;
  64         mt_init_by_seed_rand(19650218UL);
  65         i = 1;
  66         j = 0;
  67         for (k = (N > len ? N : len); k; --k) {
  68                 /* Non linear */
  69                 x[i] = ((x[i] ^ ((x[i - 1] ^ (x[i - 1] >> 30)) *
  70                        1664525UL)) + key[j] + j) & 0xffffffffUL;
  71                 if (++i >= N) {
  72                         x[0] = x[N - 1];
  73                         i = 1;
  74                 }
  75                 if (++j >= len)
  76                         j = 0;
  77         }
  78         for (k = N - 1; k; --k) {
  79                 /* Non linear */
  80                 x[i] = ((x[i] ^ ((x[i - 1] ^ (x[i - 1] >> 30)) *
  81                        1566083941UL)) - i) & 0xffffffffUL;
  82                 if (++i >= N) {
  83                         x[0] = x[N - 1];
  84                         i = 1;
  85                 }
  86         }
  87         x[0] = 0x80000000UL;
  88 }
  89
  90 void mt_init_by_seed_rand_array(void)
  91 {
  92         int i;
  93         unsigned long k[LEN_INIT];
  94         srand((unsigned int) time(NULL));
  95         for (i = 0; i < LEN_INIT; i++)
  96                 k[i] = rand();
  97         mt_init_by_seed_array(k, LEN_INIT);
  98 }
  99
 100 void mt_init_by_random_device(void)
 101 {
 102         int fd;
 103         unsigned long k[LEN_INIT];
 104         fd = open_or_die("/dev/random", O_RDONLY);
 105         read_or_die(fd, k, sizeof(unsigned long) * LEN_INIT);
 106         close(fd);
 107         mt_init_by_seed_array(k, LEN_INIT);
 108 }
 109
 110 unsigned long mt_rand_int32(void)
 111 {
 112         /* Interval [0,0xffffffff] */
 113         unsigned long y;
 114         /* Default seed */
 115         if (p0 == NULL)
 116                 mt_init_by_seed_rand(5489UL);
 117         /* Twisted feedback */
 118         y = *p0 = *pm++ ^ (((*p0 & UPPER_MASK) | (*p1 & LOWER_MASK)) >> 1) ^
 119                   (-(*p1 & 1) & MATRIX_A);
 120         p0 = p1++;
 121         if (pm == x + N)
 122                 pm = x;
 123         if (p1 == x + N)
 124                 p1 = x;
 125         /* Temper */
 126         y ^= y >> 11;
 127         y ^= y << 7 & 0x9d2c5680UL;
 128         y ^= y << 15 & 0xefc60000UL;
 129         y ^= y >> 18;
 130         return y;
 131 }
 132
 133 long mt_rand_int31(void)
 134 {
 135         /* Interval [0,0x7fffffff] */
 136         return (long) mt_rand_int32() >> 1;
 137 }
 138
 139 double mt_rand_real1(void)
 140 {
 141         /* Interval [0,1]; Divided by 2^32-1 */
 142         return mt_rand_int32() * (1.0 / 4294967295.0);
 143 }
 144
 145 double mt_rand_real2(void)
 146 {
 147         /* Interval [0,1); Divided by 2^32 */
 148         return mt_rand_int32() * (1.0 / 4294967296.0);
 149 }
 150
 151 double mt_rand_real3(void)
 152 {
 153         /* Interval (0,1); Divided by 2^32 */
 154         return (((double) mt_rand_int32()) + 0.5) * (1.0 / 4294967296.0);
 155 }
 156
 157 double mt_rand_res53(void)
 158 {
 159         /* 53-bit random number on the real interval [0,1) */
 160         unsigned long a = mt_rand_int32() >> 5, b = mt_rand_int32() >> 6;
 161         return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
 162 }
 163