#include <math.h>
#include "kmath.h"
-/**************************************
- *** Pseudo-random number generator ***
- **************************************/
-
-/*
- 64-bit Mersenne Twister pseudorandom number generator. Adapted from:
-
- http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c
-
- which was written by Takuji Nishimura and Makoto Matsumoto and released
- under the 3-clause BSD license.
-*/
-
-#define KR_NN 312
-#define KR_MM 156
-#define KR_UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */
-#define KR_LM 0x7FFFFFFFULL /* Least significant 31 bits */
-
-struct _krand_t {
- int mti, n_iset;
- double n_gset;
- krint64_t mt[KR_NN];
-};
-
-static void kr_srand0(krint64_t seed, krand_t *kr)
-{
- kr->mt[0] = seed;
- for (kr->mti = 1; kr->mti < KR_NN; ++kr->mti)
- kr->mt[kr->mti] = 6364136223846793005ULL * (kr->mt[kr->mti - 1] ^ (kr->mt[kr->mti - 1] >> 62)) + kr->mti;
-}
-
-krand_t *kr_srand(krint64_t seed)
-{
- krand_t *kr;
- kr = (krand_t*)calloc(1, sizeof(krand_t));
- kr_srand0(seed, kr);
- return kr;
-}
-
-krint64_t kr_rand(krand_t *kr)
-{
- krint64_t x;
- static const krint64_t mag01[2] = { 0, 0xB5026F5AA96619E9ULL };
- if (kr->mti >= KR_NN) {
- int i;
- if (kr->mti == KR_NN + 1) kr_srand0(5489ULL, kr);
- for (i = 0; i < KR_NN - KR_MM; ++i) {
- x = (kr->mt[i] & KR_UM) | (kr->mt[i+1] & KR_LM);
- kr->mt[i] = kr->mt[i + KR_MM] ^ (x>>1) ^ mag01[(int)(x&1)];
- }
- for (; i < KR_NN - 1; ++i) {
- x = (kr->mt[i] & KR_UM) | (kr->mt[i+1] & KR_LM);
- kr->mt[i] = kr->mt[i + (KR_MM - KR_NN)] ^ (x>>1) ^ mag01[(int)(x&1)];
- }
- x = (kr->mt[KR_NN - 1] & KR_UM) | (kr->mt[0] & KR_LM);
- kr->mt[KR_NN - 1] = kr->mt[KR_MM - 1] ^ (x>>1) ^ mag01[(int)(x&1)];
- kr->mti = 0;
- }
- x = kr->mt[kr->mti++];
- x ^= (x >> 29) & 0x5555555555555555ULL;
- x ^= (x << 17) & 0x71D67FFFEDA60000ULL;
- x ^= (x << 37) & 0xFFF7EEE000000000ULL;
- x ^= (x >> 43);
- return x;
-}
-
-double kr_normal(krand_t *kr)
-{
- if (kr->n_iset == 0) {
- double fac, rsq, v1, v2;
- do {
- v1 = 2.0 * kr_drand(kr) - 1.0;
- v2 = 2.0 * kr_drand(kr) - 1.0;
- rsq = v1 * v1 + v2 * v2;
- } while (rsq >= 1.0 || rsq == 0.0);
- fac = sqrt(-2.0 * log(rsq) / rsq);
- kr->n_gset = v1 * fac;
- kr->n_iset = 1;
- return v2 * fac;
- } else {
- kr->n_iset = 0;
- return kr->n_gset;
- }
-}
-
-#ifdef _KR_MAIN
-int main(int argc, char *argv[])
-{
- long i, N = 200000000;
- krand_t *kr;
- if (argc > 1) N = atol(argv[1]);
- kr = kr_srand(11);
- for (i = 0; i < N; ++i) kr_rand(kr);
-// for (i = 0; i < N; ++i) lrand48();
- free(kr);
- return 0;
-}
-#endif
-
/******************************
*** Non-linear programming ***
******************************/
extern "C" {
#endif
- /**********************************
- * Pseudo-random number generator *
- **********************************/
-
- typedef uint64_t krint64_t;
-
- struct _krand_t;
- typedef struct _krand_t krand_t;
-
- #define kr_drand(_kr) ((kr_rand(_kr) >> 11) * (1.0/9007199254740992.0))
- #define kr_sample(_kr, _k, _cnt) ((*(_cnt))++ < (_k)? *(_cnt) - 1 : kr_rand(_kr) % *(_cnt))
-
- krand_t *kr_srand(krint64_t seed);
- krint64_t kr_rand(krand_t *kr);
- double kr_normal(krand_t *kr);
-
/**************************
* Non-linear programming *
**************************/