--- /dev/null
+/* Random kit 1.3 */
+
+/*
+ * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org)
+ *
+ * The rk_random and rk_seed functions algorithms and the original design of
+ * the Mersenne Twister RNG:
+ *
+ * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ *
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * 3. The names of its contributors may not be used to endorse or promote
+ * products derived from this software without specific prior written
+ * permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+ * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * Original algorithm for the implementation of rk_interval function from
+ * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by
+ * Magnus Jonsson.
+ *
+ * Constants used in the rk_double implementation by Isaku Wada.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+ * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+ * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+/* static char const rcsid[] =
+ "@(#) $Jeannot: randomkit.c,v 1.28 2005/07/21 22:14:09 js Exp $"; */
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <limits.h>
+#include <math.h>
+
+#ifdef _WIN32
+/*
+ * Windows
+ * XXX: we have to use this ugly defined(__GNUC__) because it is not easy to
+ * detect the compiler used in distutils itself
+ */
+#if (defined(__GNUC__) && defined(NPY_NEEDS_MINGW_TIME_WORKAROUND))
+
+/*
+ * FIXME: ideally, we should set this to the real version of MSVCRT. We need
+ * something higher than 0x601 to enable _ftime64 and co
+ */
+#define __MSVCRT_VERSION__ 0x0700
+#include <time.h>
+#include <sys/timeb.h>
+
+/*
+ * mingw msvcr lib import wrongly export _ftime, which does not exist in the
+ * actual msvc runtime for version >= 8; we make it an alias to _ftime64, which
+ * is available in those versions of the runtime
+ */
+#define _FTIME(x) _ftime64((x))
+#else
+#include <time.h>
+#include <sys/timeb.h>
+#define _FTIME(x) _ftime((x))
+#endif
+
+#ifndef RK_NO_WINCRYPT
+/* Windows crypto */
+#ifndef _WIN32_WINNT
+#define _WIN32_WINNT 0x0400
+#endif
+#include <windows.h>
+#include <wincrypt.h>
+#endif
+
+#else
+/* Unix */
+#include <time.h>
+#include <sys/time.h>
+#include <unistd.h>
+#endif
+
+#include "randomkit.h"
+
+#ifndef RK_DEV_URANDOM
+#define RK_DEV_URANDOM "/dev/urandom"
+#endif
+
+#ifndef RK_DEV_RANDOM
+#define RK_DEV_RANDOM "/dev/random"
+#endif
+
+char *rk_strerror[RK_ERR_MAX] =
+{
+ "no error",
+ "random device unvavailable"
+};
+
+void
+rk_seed(unsigned long seed, rk_state *state)
+{
+ int pos;
+ seed &= 0xffffffffUL;
+
+ /* Knuth's PRNG as used in the Mersenne Twister reference implementation */
+ for (pos = 0; pos < RK_STATE_LEN; pos++) {
+ state->key[pos] = seed;
+ seed = (1812433253UL * (seed ^ (seed >> 30)) + pos + 1) & 0xffffffffUL;
+ }
+ state->pos = RK_STATE_LEN;
+ state->gauss = 0;
+ state->has_gauss = 0;
+ state->has_binomial = 0;
+}
+
+/* Magic Mersenne Twister constants */
+#define N 624
+#define M 397
+#define MATRIX_A 0x9908b0dfUL
+#define UPPER_MASK 0x80000000UL
+#define LOWER_MASK 0x7fffffffUL
+
+/* Slightly optimised reference implementation of the Mersenne Twister */
+unsigned long
+rk_random(rk_state *state)
+{
+ unsigned long y;
+
+ if (state->pos == RK_STATE_LEN) {
+ int i;
+
+ for (i = 0; i < N - M; i++) {
+ y = (state->key[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK);
+ state->key[i] = state->key[i+M] ^ (y>>1) ^ (-(y & 1) & MATRIX_A);
+ }
+ for (; i < N - 1; i++) {
+ y = (state->key[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK);
+ state->key[i] = state->key[i+(M-N)] ^ (y>>1) ^ (-(y & 1) & MATRIX_A);
+ }
+ y = (state->key[N - 1] & UPPER_MASK) | (state->key[0] & LOWER_MASK);
+ state->key[N - 1] = state->key[M - 1] ^ (y >> 1) ^ (-(y & 1) & MATRIX_A);
+
+ state->pos = 0;
+ }
+ y = state->key[state->pos++];
+
+ /* Tempering */
+ y ^= (y >> 11);
+ y ^= (y << 7) & 0x9d2c5680UL;
+ y ^= (y << 15) & 0xefc60000UL;
+ y ^= (y >> 18);
+
+ return y;
+}
+
+long
+rk_long(rk_state *state)
+{
+ return rk_ulong(state) >> 1;
+}
+
+unsigned long
+rk_ulong(rk_state *state)
+{
+#if ULONG_MAX <= 0xffffffffUL
+ return rk_random(state);
+#else
+ return (rk_random(state) << 32) | (rk_random(state));
+#endif
+}
+
+unsigned long
+rk_interval(unsigned long max, rk_state *state)
+{
+ unsigned long mask = max, value;
+
+ if (max == 0) {
+ return 0;
+ }
+ /* Smallest bit mask >= max */
+ mask |= mask >> 1;
+ mask |= mask >> 2;
+ mask |= mask >> 4;
+ mask |= mask >> 8;
+ mask |= mask >> 16;
+#if ULONG_MAX > 0xffffffffUL
+ mask |= mask >> 32;
+#endif
+
+ /* Search a random value in [0..mask] <= max */
+#if ULONG_MAX > 0xffffffffUL
+ if (max <= 0xffffffffUL) {
+ while ((value = (rk_random(state) & mask)) > max);
+ }
+ else {
+ while ((value = (rk_ulong(state) & mask)) > max);
+ }
+#else
+ while ((value = (rk_ulong(state) & mask)) > max);
+#endif
+ return value;
+}
+
+double
+rk_double(rk_state *state)
+{
+ /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */
+ long a = rk_random(state) >> 5, b = rk_random(state) >> 6;
+ return (a * 67108864.0 + b) / 9007199254740992.0;
+}
--- /dev/null
+/* Random kit 1.3 */
+
+/*
+ * Copyright (c) 2003-2005, Jean-Sebastien Roy (js@jeannot.org)
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+ * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+ * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+ * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+/* @(#) $Jeannot: randomkit.h,v 1.24 2005/07/21 22:14:09 js Exp $ */
+
+/*
+ * Typical use:
+ *
+ * {
+ * rk_state state;
+ * unsigned long seed = 1, random_value;
+ *
+ * rk_seed(seed, &state); // Initialize the RNG
+ * ...
+ * random_value = rk_random(&state); // Generate random values in [0..RK_MAX]
+ * }
+ *
+ * Instead of rk_seed, you can use rk_randomseed which will get a random seed
+ * from /dev/urandom (or the clock, if /dev/urandom is unavailable):
+ *
+ * {
+ * rk_state state;
+ * unsigned long random_value;
+ *
+ * rk_randomseed(&state); // Initialize the RNG with a random seed
+ * ...
+ * random_value = rk_random(&state); // Generate random values in [0..RK_MAX]
+ * }
+ */
+
+/*
+ * Useful macro:
+ * RK_DEV_RANDOM: the device used for random seeding.
+ * defaults to "/dev/urandom"
+ */
+
+#include <stddef.h>
+
+#ifndef _RANDOMKIT_
+#define _RANDOMKIT_
+
+#define RK_STATE_LEN 624
+
+typedef struct rk_state_
+{
+ unsigned long key[RK_STATE_LEN];
+ int pos;
+ int has_gauss; /* !=0: gauss contains a gaussian deviate */
+ double gauss;
+
+ /* The rk_state structure has been extended to store the following
+ * information for the binomial generator. If the input values of n or p
+ * are different than nsave and psave, then the other parameters will be
+ * recomputed. RTK 2005-09-02 */
+
+ int has_binomial; /* !=0: following parameters initialized for
+ binomial */
+ double psave;
+ long nsave;
+ double r;
+ double q;
+ double fm;
+ long m;
+ double p1;
+ double xm;
+ double xl;
+ double xr;
+ double c;
+ double laml;
+ double lamr;
+ double p2;
+ double p3;
+ double p4;
+
+}
+rk_state;
+
+typedef enum {
+ RK_NOERR = 0, /* no error */
+ RK_ENODEV = 1, /* no RK_DEV_RANDOM device */
+ RK_ERR_MAX = 2
+} rk_error;
+
+/* error strings */
+extern char *rk_strerror[RK_ERR_MAX];
+
+/* Maximum generated random value */
+#define RK_MAX 0xFFFFFFFFUL
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*
+ * Initialize the RNG state using the given seed.
+ */
+extern void rk_seed(unsigned long seed, rk_state *state);
+
+/*
+ * Returns a random unsigned long between 0 and RK_MAX inclusive
+ */
+extern unsigned long rk_random(rk_state *state);
+
+/*
+ * Returns a random long between 0 and LONG_MAX inclusive
+ */
+extern long rk_long(rk_state *state);
+
+/*
+ * Returns a random unsigned long between 0 and ULONG_MAX inclusive
+ */
+extern unsigned long rk_ulong(rk_state *state);
+
+/*
+ * Returns a random unsigned long between 0 and max inclusive.
+ */
+extern unsigned long rk_interval(unsigned long max, rk_state *state);
+
+/*
+ * Returns a random double between 0.0 and 1.0, 1.0 excluded.
+ */
+extern double rk_double(rk_state *state);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* _RANDOMKIT_ */
projects / graphviz / commitdiff