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         mt_init_by_seed_rand((unsigned long) time(NULL));
  51 }
  52
  53 void mt_init_by_seed_array(unsigned long key[], int len)
  54 {
  55         int i, j, k;
  56         mt_init_by_seed_rand(19650218UL);
  57         i = 1;
  58         j = 0;
  59         for (k = (N > len ? N : len); k; --k) {
  60                 /* Non linear */
  61                 x[i] = ((x[i] ^ ((x[i - 1] ^ (x[i - 1] >> 30)) *
  62                        1664525UL)) + key[j] + j) & 0xffffffffUL;
  63                 if (++i >= N) {
  64                         x[0] = x[N - 1];
  65                         i = 1;
  66                 }
  67                 if (++j >= len)
  68                         j = 0;
  69         }
  70         for (k = N - 1; k; --k) {
  71                 /* Non linear */
  72                 x[i] = ((x[i] ^ ((x[i - 1] ^ (x[i - 1] >> 30)) *
  73                        1566083941UL)) - i) & 0xffffffffUL;
  74                 if (++i >= N) {
  75                         x[0] = x[N - 1];
  76                         i = 1;
  77                 }
  78         }
  79         x[0] = 0x80000000UL;
  80 }
  81
  82 void mt_init_by_seed_rand_array(void)
  83 {
  84         int i;
  85         unsigned long k[LEN_INIT];
  86         srand((unsigned int) time(NULL));
  87         for (i = 0; i < LEN_INIT; i++)
  88                 k[i] = rand();
  89         mt_init_by_seed_array(k, LEN_INIT);
  90 }
  91
  92 void mt_init_by_random_device(void)
  93 {
  94         int fd;
  95         unsigned long k[LEN_INIT];
  96         fd = open_or_die("/dev/random", O_RDONLY);
  97         read_or_die(fd, k, sizeof(unsigned long) * LEN_INIT);
  98         close(fd);
  99         mt_init_by_seed_array(k, LEN_INIT);
 100 }
 101
 102 unsigned long mt_rand_int32(void)
 103 {
 104         /* Interval [0,0xffffffff] */
 105         unsigned long y;
 106         /* Default seed */
 107         if (p0 == NULL)
 108                 mt_init_by_seed_rand(5489UL);
 109         /* Twisted feedback */
 110         y = *p0 = *pm++ ^ (((*p0 & UPPER_MASK) | (*p1 & LOWER_MASK)) >> 1) ^
 111                   (-(*p1 & 1) & MATRIX_A);
 112         p0 = p1++;
 113         if (pm == x + N)
 114                 pm = x;
 115         if (p1 == x + N)
 116                 p1 = x;
 117         /* Temper */
 118         y ^= y >> 11;
 119         y ^= y << 7 & 0x9d2c5680UL;
 120         y ^= y << 15 & 0xefc60000UL;
 121         y ^= y >> 18;
 122         return y;
 123 }
 124
 125 long mt_rand_int31(void)
 126 {
 127         /* Interval [0,0x7fffffff] */
 128         return (long) mt_rand_int32() >> 1;
 129 }
 130
 131 double mt_rand_real1(void)
 132 {
 133         /* Interval [0,1]; Divided by 2^32-1 */
 134         return mt_rand_int32() * (1.0 / 4294967295.0);
 135 }
 136
 137 double mt_rand_real2(void)
 138 {
 139         /* Interval [0,1); Divided by 2^32 */
 140         return mt_rand_int32() * (1.0 / 4294967296.0);
 141 }
 142
 143 double mt_rand_real3(void)
 144 {
 145         /* Interval (0,1); Divided by 2^32 */
 146         return (((double) mt_rand_int32()) + 0.5) * (1.0 / 4294967296.0);
 147 }
 148
 149 double mt_rand_res53(void)
 150 {
 151         /* 53-bit random number on the real interval [0,1) */
 152         unsigned long a = mt_rand_int32() >> 5, b = mt_rand_int32() >> 6;
 153         return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
 154 }
 155