diff --git a/uppsrc/Core/Random.cpp b/uppsrc/Core/Random.cpp index d1d782d26..b1679f4f6 100644 --- a/uppsrc/Core/Random.cpp +++ b/uppsrc/Core/Random.cpp @@ -2,277 +2,141 @@ namespace Upp { -/* - A C-program for MT19937, with initialization improved 2002/1/26. - Coded by Takuji Nishimura and Makoto Matsumoto. +/* Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) - Before using, initialize the state by using init_genrand(seed) - or init_by_array(init_key, key_length). +To the extent possible under law, the author has dedicated all copyright +and related and neighboring rights to this software to the public domain +worldwide. This software is distributed without any warranty. - Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, - All rights reserved. +See . */ - Redistribution and use in source and binary forms, with or without - modification, are permitted provided that the following conditions - are met: +/* This is xoshiro256** 1.0, our all-purpose, rock-solid generator. It has + excellent (sub-ns) speed, a state (256 bits) that is large enough for + any parallel application, and it passes all tests we are aware of. - 1. Redistributions of source code must retain the above copyright - notice, this list of conditions and the following disclaimer. + For generating just floating-point numbers, xoshiro256+ is even faster. - 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. + The state must be seeded so that it is not everywhere zero. If you have + a 64-bit seed, we suggest to seed a splitmix64 generator and use its + output to fill s. */ - 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.keio.ac.jp/matumoto/emt.html - email: matumoto@math.keio.ac.jp -*/ - -#include - -/* Period parameters */ -#define N 624 -#define M 397 -#define MATRIX_A 0x9908b0dfUL /* constant vector a */ -#define UPPER_MASK 0x80000000UL /* most significant w-r bits */ -#define LOWER_MASK 0x7fffffffUL /* least significant r bits */ - -struct MTrand { - dword mt[N]; - int mti; /* mti==N+1 means mt[N] is not initialized */ - dword mag01[2]; -#ifdef PLATFORM_POSIX - int pid; -#endif - - void seed(); - void init_genrand(dword s); - void init_by_array(dword *init_key, int key_length); - dword genrand(); - - MTrand(); -}; - -/* initializes mt[N] with a seed */ -void MTrand::init_genrand(dword s) -{ - mt[0]= s & 0xffffffffUL; - for (mti=1; mti> 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] &= 0xffffffffUL; - /* for >32 bit machines */ - } +static force_inline uint64_t s_rotl(const uint64_t x, int k) { + return (x << k) | (x >> (64 - k)); // GCC/CLANG/MSC happily optimize this } -/* 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 MTrand::init_by_array(dword *init_key, int key_length) +force_inline +static uint64 sNext(uint64 *s) { - int i, j, k; - init_genrand(19650218UL); - i=1; j=0; - k = (N>key_length ? N : key_length); - for (; k; k--) { - mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) - + init_key[j] + j; /* non linear */ - mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ - i++; j++; - if (i>=N) { mt[0] = mt[N-1]; i=1; } - if (j>=key_length) j=0; - } - for (k=N-1; k; k--) { - mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - - i; /* non linear */ - mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ - i++; - if (i>=N) { mt[0] = mt[N-1]; i=1; } - } + const uint64_t result_starstar = s_rotl(s[1] * 5, 7) * 9; - mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ + const uint64_t t = s[1] << 17; + + s[2] ^= s[0]; + s[3] ^= s[1]; + s[1] ^= s[2]; + s[0] ^= s[3]; + + s[2] ^= t; + + s[3] = s_rotl(s[3], 45); + + return result_starstar; } -/* generates a random number on [0,0xffffffff]-interval */ -dword MTrand::genrand() +never_inline +static void sSeed(uint64 *s) { - dword y; - - if (mti >= N) { /* generate N words at one time */ - int kk; - - for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; - } - for (;kk> 1) ^ mag01[y & 0x1UL]; - } - y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); - mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; - - mti = 0; - } - - y = mt[mti++]; - - /* Tempering */ - y ^= (y >> 11); - y ^= (y << 7) & 0x9d2c5680UL; - y ^= (y << 15) & 0xefc60000UL; - y ^= (y >> 18); - - return y; -} - -#ifdef PLATFORM_WIN32 - -#define Ptr Ptr_ -#include -#undef Ptr - -#pragma comment(lib, "ole32.lib") - -#endif - -MTrand::MTrand() -{ - mti = N + 1; - mag01[0] = 0; - mag01[1] = MATRIX_A; - seed(); -#ifdef PLATFORM_POSIX - pid = getpid(); -#endif -} - -void MTrand::seed() -{ - dword seed[1024]; #ifdef PLATFORM_POSIX int fd = open("/dev/urandom", O_RDONLY); if(fd != -1) { IGNORE_RESULT( - read(fd, seed, sizeof(seed)) - ); + read(fd, s, 4 * sizeof(uint64)) + ); close(fd); } #else - for(int i = 0; i < 1024; i++) { - Uuid uuid; - CoCreateGuid((GUID *)&uuid); - seed[i] = GetHashValue(uuid); + for(int pass = 0; pass < 4; pass++) { + for(int i = 0; i < 4; i++) { + uint64 h = 0; + for(int p = 0; p < 2; p++) { + Uuid uuid; + CoCreateGuid((GUID *)&uuid); + h = (h << 32) | GetHashValue(uuid); + } + s[i] ^= h; + } } #endif - init_by_array(seed, 1024); } -#ifdef CPU_BLACKFIN -//on blackfin toolchain we dont have TLS support, but we can live here without -MTrand *sRng; -byte sRb[sizeof(MTrand)]; +force_inline +static uint64 *sState() +{ +#ifdef COMPILER_MINGW + FastMingwTls s; #else -thread_local MTrand *sRng; -thread_local byte sRb[sizeof(MTrand)]; + thread_local uint64 *s; #endif + if(!s) { + thread_local uint64 state[4]; + s = state; + sSeed(s); +#ifdef PLATFORM_POSIX + if(Thread::IsMain()) // non-main threads do now work with fork anyway + pthread_atfork(NULL, NULL, [] { sSeed(s); }); // reseed random generator after fork +#endif + } + return s; +} + +never_inline +dword Random(void) { + return (dword)sNext(sState()); +} + +never_inline +uint64 Random64(void) { + return sNext(sState()); +} + +never_inline +dword Random(dword n) { + ASSERT(n); + uint64 *s = sState(); + dword mask = (1 << SignificantBits(n)) - 1; + dword r; + do + r = (dword)sNext(s) & mask; + while(r >= n); + return r; +} + +never_inline +uint64 Random64(uint64 n) { + ASSERT(n); + uint64 *s = sState(); + uint64 mask = ((uint64)1 << SignificantBits64(n)) - 1; + uint64 r; + do + r = sNext(s) & mask; + while(r >= n); + return r; +} + +never_inline +double Randomf() +{ + return (sNext(sState()) >> 11) * (1. / (UINT64_C(1) << 53)); +} void SeedRandom() { - if(!sRng) { - sRng = new(sRb) MTrand; - } - sRng->seed(); -} - -void SeedRandom(dword *seed, int len){ - if(!sRng) { - sRng = new(sRb) MTrand; - } - sRng->init_by_array(seed, len); + sSeed(sState()); } void SeedRandom(dword seed) { - if(!sRng) { - sRng = new(sRb) MTrand; - } - sRng->init_genrand(seed); -} - -dword Random() -{ - if(!sRng) { - sRng = new(sRb) MTrand; - } -#ifdef PLATFORM_POSIX // Be fork safe... - int pid = getpid(); - if(sRng->pid != pid) { - sRng->seed(); - sRng->pid = pid; - } -#endif - return sRng->genrand(); -} - -dword Random(dword n) -{ - ASSERT(n); - dword mask = n; - mask |= mask >> 1; - mask |= mask >> 2; - mask |= mask >> 4; - mask |= mask >> 8; - mask |= mask >> 16; - - dword r; - do - r = Random() & mask; - while(r >= n); - return r; -} - -qword Random64() -{ - return MAKEQWORD(Random(), Random()); -} - -qword Random64(qword n) -{ - qword mask = n, r; - mask |= mask >> 1; mask |= mask >> 2; - mask |= mask >> 4; mask |= mask >> 8; - mask |= mask >> 16; mask |= mask >> 32; - - do - r = Random64() & mask; - while(r >= n); - return r; -} - -double Randomf() -{ - return Random64(I64(4503599627370496)) / 4503599627370496.0; + uint64 *s = sState(); + s[0] = s[1] = s[2] = s[3] = seed + 12345678; } } diff --git a/uppsrc/Core/Util.h b/uppsrc/Core/Util.h index 3a4d401ed..8add7da46 100644 --- a/uppsrc/Core/Util.h +++ b/uppsrc/Core/Util.h @@ -265,11 +265,11 @@ inline void Dbl_Self(T *x) #define ZeroArray(x) memset((x), 0, sizeof(x)) -dword Random(); -dword Random(dword n); -qword Random64(); -qword Random64(qword n); -double Randomf(); +dword Random(); +dword Random(dword n); +uint64 Random64(); +uint64 Random64(uint64 n); +double Randomf(); void SeedRandom(dword *seed, int len); void SeedRandom(dword seed); diff --git a/uppsrc/Core/src.tpp/MathUtil_en-us.tpp b/uppsrc/Core/src.tpp/MathUtil_en-us.tpp index 780cc8edf..b6fe0c9ac 100644 --- a/uppsrc/Core/src.tpp/MathUtil_en-us.tpp +++ b/uppsrc/Core/src.tpp/MathUtil_en-us.tpp @@ -1,5 +1,4 @@ topic "Math utility"; -[2 $$0,0#00000000000000000000000000000000:Default] [i448;a25;kKO9;2 $$1,0#37138531426314131252341829483380:class] [l288;2 $$2,0#27521748481378242620020725143825:desc] [0 $$3,0#96390100711032703541132217272105:end] @@ -9,6 +8,7 @@ topic "Math utility"; [l288;i1121;b17;O9;~~~.1408;2 $$7,0#10431211400427159095818037425705:param] [i448;b42;O9;2 $$8,8#61672508125594000341940100500538:tparam] [b42;2 $$9,9#13035079074754324216151401829390:normal] +[2 $$0,0#00000000000000000000000000000000:Default] [{_} [ {{10000@(113.42.0) [s0;%% [*@7;4 Math utility]]}}&] [s3; &] @@ -81,22 +81,37 @@ topic "Math utility"; [s3; &] [s5;:Random`(`): [_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:dword`:`:typedef^ d word]_[* Random]()&] -[s2;%% Returns the next random generated number. MT safe and efficient -variant of Merssene twister MT19937 by Takuji Nishimura and Makoto -Matsumoto is used to get random numbers, with the seed obtained -from host platform specific resources (`'/dev/urandom`' in Posix -systems, GUID generator in Win32).&] +[s2;%% Returns the next random generated number. Algorithm used is +xoshiro256`*`*, with seed obtained from host platform specific +resources (`'/dev/urandom`' in Posix systems, GUID generator +in Win32).&] [s3;%% &] [s4; &] -[s5;:SeedRandom`(dword`*`,int`): [@(0.0.255) void]_[* SeedRandom]([_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:dword`:`:typedef^ d -word]_`*[*@3 seed],[@(0.0.255) int]_[*@3 len])&] -[s2;%% Seeds random with buffer.&] +[s5;:Random`(dword`): [_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:dword`:`:typedef^ d +word]_[* Random]([_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:dword`:`:typedef^ d +word]_[*@3 n])&] +[s2; Returns random generated number smaller than [*@3 n].&] +[s3;%% &] +[s4; &] +[s5;:Random64`(`): [_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:qword`:`:typedef^ q +word]_[* Random64]()&] +[s2;%% Returns 64bit random number.&] +[s3; &] +[s4; &] +[s5;:Upp`:`:Random64`(Upp`:`:uint64`): [_^Upp`:`:uint64^ uint64]_[* Random64]([_^Upp`:`:uint64^ u +int64]_[*@3 n])&] +[s2; [%% Returns the 64bit random number] smaller than [*@3 n].&] +[s3; &] +[s4; &] +[s5;:Randomf`(`): [@(0.0.255) double]_[* Randomf]()&] +[s2;%% Returns the floating point number <0, 1).&] [s3;%% &] [s4; &] [s5;:SeedRandom`(dword`): [@(0.0.255) void]_[* SeedRandom]([_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:dword`:`:typedef^ d word]_[*@3 seed])&] [s2;%% Seeds random with single value [%-*@3 seed]. This is good to -get always the same sequence of numbers (for the same seed).&] +get always the same sequence of numbers (for the same seed). +Usually used for testing.&] [s3;%% &] [s4; &] [s5;:Upp`:`:SeedRandom`(`): [@(0.0.255) void]_[* SeedRandom]()&] @@ -106,27 +121,6 @@ Random after seeding it with fixed value.&] [s4; &] [s5;:ZeroArray`(x`): [* ZeroArray]([*@3 x])&] [s2;%% Fills C array [%-*@3 x] with zeros.&] -[s3;%% &] -[s4; &] -[s5;:Random`(dword`): [_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:dword`:`:typedef^ d -word]_[* Random]([_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:dword`:`:typedef^ d -word]_[*@3 n])&] -[s2; Returns the next random generated number smaller than [*@3 n].&] -[s3;%% &] -[s4; &] -[s5;:Random64`(`): [_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:qword`:`:typedef^ q -word]_[* Random64]()&] -[s2;%% Returns the 64bit random number.&] -[s3; &] -[s4; &] -[s5;:Random64`(qword`): [_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:qword`:`:typedef^ q -word]_[* Random64]([_^topic`:`/`/Core`/src`/PrimitiveDataTypes`$en`-us`#Upp`:`:qword`:`:typedef^ q -word]_[*@3 n])&] -[s2; [%% Returns the 64bit random number] smaller than [*@3 n].&] -[s3;%% &] -[s4; &] -[s5;:Randomf`(`): [@(0.0.255) double]_[* Randomf]()&] -[s2;%% Returns the floating point number 0..1.&] [s3; &] [s4; &] [s5;:sqr`(double`): [@(0.0.255) double]_[* sqr]_([@(0.0.255) double]_[*@3 a])&]