ref: b09398e98563c23b6988e77b7b5ab84c0e03ceb7
parent: 0f11cb9cc11baa2960aa47eeb5e2d8a28ce271b3
author: Sigrid Solveig Haflínudóttir <[email protected]>
date: Sat Apr 15 20:08:08 EDT 2023
replace mt19937 with Cocoa with Love modified 64-bit version
--- /dev/null
+++ b/3rd/mt19937-64.c
@@ -1,0 +1,178 @@
+/*
+ A C-program for MT19937-64 (2004/9/29 version).
+ Coded by Takuji Nishimura and Makoto Matsumoto.
+
+ This is a 64-bit version of Mersenne Twister pseudorandom number
+ generator.
+
+ Before using, initialize the state by using init_genrand64(seed)
+ or init_by_array64(init_key, key_length).
+
+ Copyright (C) 2004, 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.
+
+ References:
+ T. Nishimura, ``Tables of 64-bit Mersenne Twisters''
+ ACM Transactions on Modeling and
+ Computer Simulation 10. (2000) 348--357.
+ M. Matsumoto and T. Nishimura,
+ ``Mersenne Twister: a 623-dimensionally equidistributed
+ uniform pseudorandom number generator''
+ ACM Transactions on Modeling and
+ Computer Simulation 8. (Jan. 1998) 3--30.
+
+ Any feedback is very welcome.
+ http://www.math.hiroshima-u.ac.jp/~m-mat/MT/emt.html
+ email: m-mat @ math.sci.hiroshima-u.ac.jp (remove spaces)
+*/
+
+#include "mt19937-64.h"
+
+#define NN 312
+#define MM 156
+#define MATRIX_A 0xB5026F5AA96619E9ULL
+#define UM 0xFFFFFFFF80000000ULL /* Most significant 33 bits */
+#define LM 0x7FFFFFFFULL /* Least significant 31 bits */
+
+/* initializes mt[NN] with a seed */
+void
+init_genrand64(mt19937_64 *context, unsigned long long seed)
+{+ context->mt[0] = seed;
+ for(context->mti=1; context->mti<NN; context->mti++)
+ context->mt[context->mti] = (6364136223846793005ULL * (context->mt[context->mti-1] ^ (context->mt[context->mti-1] >> 62)) + context->mti);
+}
+
+/* initialize by an array with array-length */
+/* init_key is the array for initializing keys */
+/* key_length is its length */
+void
+init_by_array64(mt19937_64 *context, unsigned long long init_key[], unsigned long long key_length)
+{+ unsigned long long i, j, k;
+
+ init_genrand64(context, 19650218ULL);
+ i = 1;
+ j = 0;
+ k = NN > key_length ? NN : key_length;
+ for(; k; k--){+ context->mt[i] = context->mt[i] ^ ((context->mt[i-1] ^ (context->mt[i-1] >> 62)) * 3935559000370003845ULL) + init_key[j] + j; /* non linear */
+ i++;
+ j++;
+ if(i >= NN){+ context->mt[0] = context->mt[NN-1];
+ i = 1;
+ }
+ if(j >= key_length)
+ j = 0;
+ }
+ for(k = NN-1; k; k--){+ context->mt[i] = context->mt[i] ^ ((context->mt[i-1] ^ (context->mt[i-1] >> 62)) * 2862933555777941757ULL) - i; /* non linear */
+ i++;
+ if(i >= NN){+ context->mt[0] = context->mt[NN-1];
+ i = 1;
+ }
+ }
+
+ context->mt[0] = 1ULL << 63; /* MSB is 1; assuring non-zero initial array */
+}
+
+/* generates a random number on [0, 2^64-1]-interval */
+unsigned long long
+genrand64_int64(mt19937_64 *context)
+{+ /* This is the altered Cocoa with Love implementation. */
+ int i, j;
+ unsigned long long result;
+
+ if(context->mti >= NN){ /* generate NN words at one time */+ int mid = NN / 2;
+ unsigned long long stateMid = context->mt[mid];
+ unsigned long long x;
+ unsigned long long y;
+
+ /* NOTE: this "untwist" code is modified from the original to improve
+ * performance, as described here:
+ * http://www.cocoawithlove.com/blog/2016/05/19/random-numbers.html
+ * These modifications are offered for use under the original icense at
+ * the top of this file.
+ */
+ for(i = 0, j = mid; i != mid - 1; i++, j++){+ x = (context->mt[i] & UM) | (context->mt[i + 1] & LM);
+ context->mt[i] = context->mt[i + mid] ^ (x >> 1) ^ ((context->mt[i + 1] & 1) * MATRIX_A);
+ y = (context->mt[j] & UM) | (context->mt[j + 1] & LM);
+ context->mt[j] = context->mt[j - mid] ^ (y >> 1) ^ ((context->mt[j + 1] & 1) * MATRIX_A);
+ }
+ x = (context->mt[mid - 1] & UM) | (stateMid & LM);
+ context->mt[mid - 1] = context->mt[NN - 1] ^ (x >> 1) ^ ((stateMid & 1) * MATRIX_A);
+ y = (context->mt[NN - 1] & UM) | (context->mt[0] & LM);
+ context->mt[NN - 1] = context->mt[mid - 1] ^ (y >> 1) ^ ((context->mt[0] & 1) * MATRIX_A);
+ context->mti = 0;
+ }
+
+ result = context->mt[context->mti];
+ context->mti = context->mti + 1;
+
+ result ^= (result >> 29) & 0x5555555555555555ULL;
+ result ^= (result << 17) & 0x71D67FFFEDA60000ULL;
+ result ^= (result << 37) & 0xFFF7EEE000000000ULL;
+ result ^= (result >> 43);
+
+ return result;
+}
+
+/* generates a random number on [0, 2^63-1]-interval */
+long long
+genrand64_int63(mt19937_64 *context)
+{+ return (long long)(genrand64_int64(context) >> 1);
+}
+
+/* generates a random number on [0,1]-real-interval */
+double
+genrand64_real1(mt19937_64 *context)
+{+ return (genrand64_int64(context) >> 11) * (1.0/9007199254740991.0);
+}
+
+/* generates a random number on [0,1)-real-interval */
+double
+genrand64_real2(mt19937_64 *context)
+{+ return (genrand64_int64(context) >> 11) * (1.0/9007199254740992.0);
+}
+
+/* generates a random number on (0,1)-real-interval */
+double
+genrand64_real3(mt19937_64 *context)
+{+ return ((genrand64_int64(context) >> 12) + 0.5) * (1.0/4503599627370496.0);
+}
--- /dev/null
+++ b/3rd/mt19937-64.h
@@ -1,0 +1,14 @@
+typedef struct mt19937_64 mt19937_64;
+
+struct mt19937_64 {+ unsigned long long mt[312];
+ int mti;
+};
+
+void init_genrand64(mt19937_64 *context, unsigned long long seed);
+void init_by_array64(mt19937_64 *context, unsigned long long init_key[], unsigned long long key_length);
+unsigned long long genrand64_int64(mt19937_64 *context);
+long long genrand64_int63(mt19937_64 *context);
+double genrand64_real1(mt19937_64 *context);
+double genrand64_real2(mt19937_64 *context);
+double genrand64_real3(mt19937_64 *context);
--- a/3rd/mt19937ar.c
+++ /dev/null
@@ -1,171 +1,0 @@
-/*
- A C-program for MT19937, with initialization improved 2002/1/26.
- Coded by Takuji Nishimura and Makoto Matsumoto.
-
- Before using, initialize the state by using init_genrand(seed)
- or init_by_array(init_key, key_length).
-
- 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.
-
-
- Any feedback is very welcome.
- http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
- email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
-*/
-
-/* Period parameters */
-#define mtN 624
-#define mtM 397
-#define MATRIX_A 0x9908b0dfU /* constant vector a */
-#define UPPER_MASK 0x80000000U /* most significant w-r bits */
-#define LOWER_MASK 0x7fffffffU /* least significant r bits */
-
-static uint32_t mt[mtN]; /* the array for the state vector */
-static int mti=mtN+1; /* mti==mtN+1 means mt[mtN] is not initialized */
-
-/* initializes mt[mtN] with a seed */
-void init_genrand(uint32_t s)
-{- mt[0]= s & 0xffffffffU;
- for (mti=1; mti<mtN; mti++) {- mt[mti] =
- (1812433253U * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
- /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
- /* In the previous versions, MSBs of the seed affect */
- /* only MSBs of the array mt[]. */
- /* 2002/01/09 modified by Makoto Matsumoto */
- mt[mti] &= 0xffffffffU;
- /* for >32 bit machines */
- }
-}
-
-/* initialize by an array with array-length */
-/* init_key is the array for initializing keys */
-/* key_length is its length */
-/* slight change for C++, 2004/2/26 */
-void init_by_array(uint32_t init_key[], int key_length)
-{- int i, j, k;
- init_genrand(19650218U);
- i=1; j=0;
- k = (mtN>key_length ? mtN : key_length);
- for (; k; k--) {- mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525U))
- + init_key[j] + j; /* non linear */
- mt[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
- i++; j++;
- if (i>=mtN) { mt[0] = mt[mtN-1]; i=1; }- if (j>=key_length) j=0;
- }
- for (k=mtN-1; k; k--) {- mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941U))
- - i; /* non linear */
- mt[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
- i++;
- if (i>=mtN) { mt[0] = mt[mtN-1]; i=1; }- }
-
- mt[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
-}
-
-/* generates a random number on [0,0xffffffff]-interval */
-uint32_t genrand_int32(void)
-{- uint32_t y;
- static uint32_t mag01[2]={0x0U, MATRIX_A};- /* mag01[x] = x * MATRIX_A for x=0,1 */
-
- if (mti >= mtN) { /* generate mtN words at one time */- int kk;
-
- if (mti == mtN+1) /* if init_genrand() has not been called, */
- init_genrand(5489U); /* a default initial seed is used */
-
- for (kk=0;kk<mtN-mtM;kk++) {- y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
- mt[kk] = mt[kk+mtM] ^ (y >> 1) ^ mag01[y & 0x1U];
- }
- for (;kk<mtN-1;kk++) {- y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
- mt[kk] = mt[kk+(mtM-mtN)] ^ (y >> 1) ^ mag01[y & 0x1U];
- }
- y = (mt[mtN-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
- mt[mtN-1] = mt[mtM-1] ^ (y >> 1) ^ mag01[y & 0x1U];
-
- mti = 0;
- }
-
- y = mt[mti++];
-
- /* Tempering */
- y ^= (y >> 11);
- y ^= (y << 7) & 0x9d2c5680U;
- y ^= (y << 15) & 0xefc60000U;
- y ^= (y >> 18);
-
- return y;
-}
-
-#if 0
-/* generates a random number on [0,0x7fffffff]-interval */
-long genrand_int31(void)
-{- return (long)(genrand_int32()>>1);
-}
-
-/* generates a random number on [0,1]-real-interval */
-double genrand_real1(void)
-{- return genrand_int32()*(1.0/4294967295.0);
- /* divided by 2^32-1 */
-}
-
-/* generates a random number on [0,1)-real-interval */
-double genrand_real2(void)
-{- return genrand_int32()*(1.0/4294967296.0);
- /* divided by 2^32 */
-}
-
-/* generates a random number on (0,1)-real-interval */
-double genrand_real3(void)
-{- return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
- /* divided by 2^32 */
-}
-
-/* generates a random number on [0,1) with 53-bit resolution*/
-double genrand_res53(void)
-{- uint32_t a=genrand_int32()>>5, b=genrand_int32()>>6;
- return(a*67108864.0+b)*(1.0/9007199254740992.0);
-}
-#endif
--- a/Makefile
+++ b/Makefile
@@ -54,6 +54,7 @@
3rd/mp/u16.o\
3rd/mp/u32.o\
3rd/mp/u64.o\
+ 3rd/mt19937-64.o\
.PHONY: all default test clean
--- a/builtins.c
+++ b/builtins.c
@@ -372,11 +372,12 @@
BUILTIN("rand", rand) {USED(args); USED(nargs);
+ uint64_t x = genrand_int63();
fixnum_t r;
#ifdef BITS64
- r = ((((uint64_t)random())<<32) | random()) & 0x1fffffffffffffffLL;
+ r = x >> 3;
#else
- r = random() & 0x1fffffff;
+ r = x >> (32+3);
#endif
return fixnum(r);
}
@@ -384,31 +385,25 @@
BUILTIN("rand.uint32", rand_uint32) {USED(args); USED(nargs);
- uint32_t r = random();
-#ifdef BITS64
- return fixnum(r);
-#else
- return mk_uint32(r);
-#endif
+ return mk_uint32(genrand_uint32());
}
BUILTIN("rand.uint64", rand_uint64) {USED(args); USED(nargs);
- uint64_t r = (((uint64_t)random())<<32) | random();
- return mk_uint64(r);
+ return mk_uint64(genrand_uint64());
}
BUILTIN("rand.double", rand_double) {USED(args); USED(nargs);
- return mk_double(rand_double());
+ return mk_double(genrand_double());
}
BUILTIN("rand.float", rand_float) {USED(args); USED(nargs);
- return mk_float(rand_float());
+ return mk_float(genrand_double());
}
#define BUILTIN_(lname, cname) \
--- a/llt/llt.h
+++ b/llt/llt.h
@@ -81,15 +81,11 @@
uint32_t memhash32(const char* buf, size_t n);
/* random.c */
-#define random() genrand_int32()
-#define srandom(n) init_genrand(n)
-double rand_double(void);
-float rand_float(void);
-double randn(void);
void randomize(void);
-uint32_t genrand_int32(void);
-void init_genrand(uint32_t s);
-uint64_t i64time(void);
+double genrand_double(void);
+uint64_t genrand_uint64(void);
+uint32_t genrand_uint32(void);
+int64_t genrand_int63(void);
/* utils.c */
char *uint2str(char *dest, size_t len, uint64_t num, uint32_t base);
--- a/llt/random.c
+++ b/llt/random.c
@@ -2,58 +2,37 @@
random numbers
*/
#include "llt.h"
-#include "ieee754.h"
-#include "mt19937ar.c"
+#include "mt19937-64.h"
-double
-rand_double()
-{- union ieee754_double d;
+static mt19937_64 ctx;
- d.ieee.mantissa0 = genrand_int32();
- d.ieee.mantissa1 = genrand_int32();
- d.ieee.negative = 0;
- d.ieee.exponent = IEEE754_DOUBLE_BIAS + 0; /* 2^0 */
- return d.d - 1.0;
+uint64_t
+genrand_uint64(void)
+{+ return genrand64_int64(&ctx);
}
-float
-rand_float()
+uint32_t
+genrand_uint32(void)
{- union ieee754_float f;
+ return genrand64_int64(&ctx) >> 32;
+}
- f.ieee.mantissa = genrand_int32();
- f.ieee.negative = 0;
- f.ieee.exponent = IEEE754_FLOAT_BIAS + 0; /* 2^0 */
- return f.f - 1.0;
+int64_t
+genrand_int63(void)
+{+ return genrand64_int63(&ctx);
}
double
-randn()
+genrand_double(void)
{- double s, vre, vim, ure, uim;
- static double next = -42;
-
- if(next != -42){- s = next;
- next = -42;
- return s;
- }
- do{- ure = rand_double();
- uim = rand_double();
- vre = 2*ure - 1;
- vim = 2*uim - 1;
- s = vre*vre + vim*vim;
- }while(s >= 1);
- s = sqrt(-2*log(s)/s);
- next = s * vre;
- return s * vim;
+ return genrand64_real1(&ctx);
}
void
randomize(void)
{- uint64_t tm = i64time();
- init_by_array((uint32_t*)&tm, 2);
+ unsigned long long tm = i64time();
+ init_by_array64(&ctx, &tm, 1);
}
--- a/mkfile
+++ b/mkfile
@@ -25,6 +25,7 @@
iostream.$O\
string.$O\
table.$O\
+ 3rd/mt19937-64.$O\
3rd/wcwidth.$O\
llt/bitvector-ops.$O\
llt/bitvector.$O\