src/randomkit.c

   1 /* Random kit 1.3 */
   2
   3 /*
   4  * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org)
   5  *
   6  * The rk_random and rk_seed functions algorithms and the original design of
   7  * the Mersenne Twister RNG:
   8  *
   9  *   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
  10  *   All rights reserved.
  11  *
  12  *   Redistribution and use in source and binary forms, with or without
  13  *   modification, are permitted provided that the following conditions
  14  *   are met:
  15  *
  16  *   1. Redistributions of source code must retain the above copyright
  17  *   notice, this list of conditions and the following disclaimer.
  18  *
  19  *   2. Redistributions in binary form must reproduce the above copyright
  20  *   notice, this list of conditions and the following disclaimer in the
  21  *   documentation and/or other materials provided with the distribution.
  22  *
  23  *   3. The names of its contributors may not be used to endorse or promote
  24  *   products derived from this software without specific prior written
  25  *   permission.
  26  *
  27  *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  28  *   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  29  *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  30  *   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
  31  *   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
  32  *   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  33  *   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  34  *   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
  35  *   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
  36  *   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
  37  *   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  38  *
  39  * Original algorithm for the implementation of rk_interval function from
  40  * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by
  41  * Magnus Jonsson.
  42  *
  43  * Constants used in the rk_double implementation by Isaku Wada.
  44  *
  45  * Permission is hereby granted, free of charge, to any person obtaining a
  46  * copy of this software and associated documentation files (the
  47  * "Software"), to deal in the Software without restriction, including
  48  * without limitation the rights to use, copy, modify, merge, publish,
  49  * distribute, sublicense, and/or sell copies of the Software, and to
  50  * permit persons to whom the Software is furnished to do so, subject to
  51  * the following conditions:
  52  *
  53  * The above copyright notice and this permission notice shall be included
  54  * in all copies or substantial portions of the Software.
  55  *
  56  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  57  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  58  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  59  * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
  60  * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  61  * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  62  * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  63  */
  64
  65 /* static char const rcsid[] =
  66   "@(#) $Jeannot: randomkit.c,v 1.28 2005/07/21 22:14:09 js Exp $"; */
  67 #include <stddef.h>
  68 #include <stdio.h>
  69 #include <stdlib.h>
  70 #include <errno.h>
  71 #include <limits.h>
  72 #include <math.h>
  73
  74 #ifdef _WIN32
  75 /*
  76  * Windows
  77  * XXX: we have to use this ugly defined(__GNUC__) because it is not easy to
  78  * detect the compiler used in distutils itself
  79  */
  80 #if (defined(__GNUC__) && defined(NPY_NEEDS_MINGW_TIME_WORKAROUND))
  81
  82 /*
  83  * FIXME: ideally, we should set this to the real version of MSVCRT. We need
  84  * something higher than 0x601 to enable _ftime64 and co
  85  */
  86 #define __MSVCRT_VERSION__ 0x0700
  87 #include <time.h>
  88 #include <sys/timeb.h>
  89
  90 /*
  91  * mingw msvcr lib import wrongly export _ftime, which does not exist in the
  92  * actual msvc runtime for version >= 8; we make it an alias to _ftime64, which
  93  * is available in those versions of the runtime
  94  */
  95 #define _FTIME(x) _ftime64((x))
  96 #else
  97 #include <time.h>
  98 #include <sys/timeb.h>
  99 #define _FTIME(x) _ftime((x))
 100 #endif
 101
 102 #ifndef RK_NO_WINCRYPT
 103 /* Windows crypto */
 104 #ifndef _WIN32_WINNT
 105 #define _WIN32_WINNT 0x0400
 106 #endif
 107 #include <windows.h>
 108 #include <wincrypt.h>
 109 #endif
 110
 111 #else
 112 /* Unix */
 113 #include <time.h>
 114 #include <sys/time.h>
 115 #include <unistd.h>
 116 #endif
 117
 118 #include "randomkit.h"
 119
 120 #ifndef RK_DEV_URANDOM
 121 #define RK_DEV_URANDOM "/dev/urandom"
 122 #endif
 123
 124 #ifndef RK_DEV_RANDOM
 125 #define RK_DEV_RANDOM "/dev/random"
 126 #endif
 127
 128 char *rk_strerror[RK_ERR_MAX] =
 129 {
 130     "no error",
 131     "random device unvavailable"
 132 };
 133
 134 /* static functions */
 135 static unsigned long rk_hash(unsigned long key);
 136
 137 void
 138 rk_seed(unsigned long seed, rk_state *state)
 139 {
 140     int pos;
 141     seed &= 0xffffffffUL;
 142
 143     /* Knuth's PRNG as used in the Mersenne Twister reference implementation */
 144     for (pos = 0; pos < RK_STATE_LEN; pos++) {
 145         state->key[pos] = seed;
 146         seed = (1812433253UL * (seed ^ (seed >> 30)) + pos + 1) & 0xffffffffUL;
 147     }
 148     state->pos = RK_STATE_LEN;
 149     state->gauss = 0;
 150     state->has_gauss = 0;
 151     state->has_binomial = 0;
 152 }
 153
 154 /* Thomas Wang 32 bits integer hash function */
 155 unsigned long
 156 rk_hash(unsigned long key)
 157 {
 158     key += ~(key << 15);
 159     key ^=  (key >> 10);
 160     key +=  (key << 3);
 161     key ^=  (key >> 6);
 162     key += ~(key << 11);
 163     key ^=  (key >> 16);
 164     return key;
 165 }
 166
 167 rk_error
 168 rk_randomseed(rk_state *state)
 169 {
 170 #ifndef _WIN32
 171     struct timeval tv;
 172 #else
 173     struct _timeb  tv;
 174 #endif
 175     int i;
 176
 177     if (rk_devfill(state->key, sizeof(state->key), 0) == RK_NOERR) {
 178         /* ensures non-zero key */
 179         state->key[0] |= 0x80000000UL;
 180         state->pos = RK_STATE_LEN;
 181         state->gauss = 0;
 182         state->has_gauss = 0;
 183         state->has_binomial = 0;
 184
 185         for (i = 0; i < 624; i++) {
 186             state->key[i] &= 0xffffffffUL;
 187         }
 188         return RK_NOERR;
 189     }
 190
 191 #ifndef _WIN32
 192     gettimeofday(&tv, NULL);
 193     rk_seed(rk_hash(getpid()) ^ rk_hash(tv.tv_sec) ^ rk_hash(tv.tv_usec)
 194             ^ rk_hash(clock()), state);
 195 #else
 196     _FTIME(&tv);
 197     rk_seed(rk_hash(tv.time) ^ rk_hash(tv.millitm) ^ rk_hash(clock()), state);
 198 #endif
 199
 200     return RK_ENODEV;
 201 }
 202
 203 /* Magic Mersenne Twister constants */
 204 #define N 624
 205 #define M 397
 206 #define MATRIX_A 0x9908b0dfUL
 207 #define UPPER_MASK 0x80000000UL
 208 #define LOWER_MASK 0x7fffffffUL
 209
 210 /* Slightly optimised reference implementation of the Mersenne Twister */
 211 unsigned long
 212 rk_random(rk_state *state)
 213 {
 214     unsigned long y;
 215
 216     if (state->pos == RK_STATE_LEN) {
 217         int i;
 218
 219         for (i = 0; i < N - M; i++) {
 220             y = (state->key[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK);
 221             state->key[i] = state->key[i+M] ^ (y>>1) ^ (-(y & 1) & MATRIX_A);
 222         }
 223         for (; i < N - 1; i++) {
 224             y = (state->key[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK);
 225             state->key[i] = state->key[i+(M-N)] ^ (y>>1) ^ (-(y & 1) & MATRIX_A);
 226         }
 227         y = (state->key[N - 1] & UPPER_MASK) | (state->key[0] & LOWER_MASK);
 228         state->key[N - 1] = state->key[M - 1] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A);
 229
 230         state->pos = 0;
 231     }
 232     y = state->key[state->pos++];
 233
 234     /* Tempering */
 235     y ^= (y >> 11);
 236     y ^= (y << 7) & 0x9d2c5680UL;
 237     y ^= (y << 15) & 0xefc60000UL;
 238     y ^= (y >> 18);
 239
 240     return y;
 241 }
 242
 243 long
 244 rk_long(rk_state *state)
 245 {
 246     return rk_ulong(state) >> 1;
 247 }
 248
 249 unsigned long
 250 rk_ulong(rk_state *state)
 251 {
 252 #if ULONG_MAX <= 0xffffffffUL
 253     return rk_random(state);
 254 #else
 255     return (rk_random(state) << 32) | (rk_random(state));
 256 #endif
 257 }
 258
 259 unsigned long
 260 rk_interval(unsigned long max, rk_state *state)
 261 {
 262     unsigned long mask = max, value;
 263
 264     if (max == 0) {
 265         return 0;
 266     }
 267     /* Smallest bit mask >= max */
 268     mask |= mask >> 1;
 269     mask |= mask >> 2;
 270     mask |= mask >> 4;
 271     mask |= mask >> 8;
 272     mask |= mask >> 16;
 273 #if ULONG_MAX > 0xffffffffUL
 274     mask |= mask >> 32;
 275 #endif
 276
 277     /* Search a random value in [0..mask] <= max */
 278 #if ULONG_MAX > 0xffffffffUL
 279     if (max <= 0xffffffffUL) {
 280         while ((value = (rk_random(state) & mask)) > max);
 281     }
 282     else {
 283         while ((value = (rk_ulong(state) & mask)) > max);
 284     }
 285 #else
 286     while ((value = (rk_ulong(state) & mask)) > max);
 287 #endif
 288     return value;
 289 }
 290
 291 double
 292 rk_double(rk_state *state)
 293 {
 294     /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */
 295     long a = rk_random(state) >> 5, b = rk_random(state) >> 6;
 296     return (a * 67108864.0 + b) / 9007199254740992.0;
 297 }
 298
 299 void
 300 rk_fill(void *buffer, size_t size, rk_state *state)
 301 {
 302     unsigned long r;
 303     unsigned char *buf = buffer;
 304
 305     for (; size >= 4; size -= 4) {
 306         r = rk_random(state);
 307         *(buf++) = r & 0xFF;
 308         *(buf++) = (r >> 8) & 0xFF;
 309         *(buf++) = (r >> 16) & 0xFF;
 310         *(buf++) = (r >> 24) & 0xFF;
 311     }
 312
 313     if (!size) {
 314         return;
 315     }
 316     r = rk_random(state);
 317     for (; size; r >>= 8, size --) {
 318         *(buf++) = (unsigned char)(r & 0xFF);
 319     }
 320 }
 321
 322 rk_error
 323 rk_devfill(void *buffer, size_t size, int strong)
 324 {
 325 #ifndef _WIN32
 326     FILE *rfile;
 327     int done;
 328
 329     if (strong) {
 330         rfile = fopen(RK_DEV_RANDOM, "rb");
 331     }
 332     else {
 333         rfile = fopen(RK_DEV_URANDOM, "rb");
 334     }
 335     if (rfile == NULL) {
 336         return RK_ENODEV;
 337     }
 338     done = fread(buffer, size, 1, rfile);
 339     fclose(rfile);
 340     if (done) {
 341         return RK_NOERR;
 342     }
 343 #else
 344
 345 #ifndef RK_NO_WINCRYPT
 346     HCRYPTPROV hCryptProv;
 347     BOOL done;
 348
 349     if (!CryptAcquireContext(&hCryptProv, NULL, NULL, PROV_RSA_FULL,
 350             CRYPT_VERIFYCONTEXT) || !hCryptProv) {
 351         return RK_ENODEV;
 352     }
 353     done = CryptGenRandom(hCryptProv, size, (unsigned char *)buffer);
 354     CryptReleaseContext(hCryptProv, 0);
 355     if (done) {
 356         return RK_NOERR;
 357     }
 358 #endif
 359
 360 #endif
 361     return RK_ENODEV;
 362 }
 363
 364 rk_error
 365 rk_altfill(void *buffer, size_t size, int strong, rk_state *state)
 366 {
 367     rk_error err;
 368
 369     err = rk_devfill(buffer, size, strong);
 370     if (err) {
 371         rk_fill(buffer, size, state);
 372     }
 373     return err;
 374 }
 375
 376 double
 377 rk_gauss(rk_state *state)
 378 {
 379     if (state->has_gauss) {
 380         const double tmp = state->gauss;
 381         state->gauss = 0;
 382         state->has_gauss = 0;
 383         return tmp;
 384     }
 385     else {
 386         double f, x1, x2, r2;
 387
 388         do {
 389             x1 = 2.0*rk_double(state) - 1.0;
 390             x2 = 2.0*rk_double(state) - 1.0;
 391             r2 = x1*x1 + x2*x2;
 392         }
 393         while (r2 >= 1.0 || r2 == 0.0);
 394
 395         /* Box-Muller transform */
 396         f = sqrt(-2.0*log(r2)/r2);
 397         /* Keep for next call */
 398         state->gauss = f*x1;
 399         state->has_gauss = 1;
 400         return f*x2;
 401     }
 402 }