mirror of
https://github.com/ultimatepp/ultimatepp.git
synced 2026-08-25 14:22:54 -06:00
Core: Random refactored with xoshiro256**
git-svn-id: svn://ultimatepp.org/upp/trunk@13393 f0d560ea-af0d-0410-9eb7-867de7ffcac7
This commit is contained in:
cxl
parent
adc94a71c2
commit
1cbd73441c
3 changed files with 133 additions and 275 deletions
|
|
@ -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 <http://creativecommons.org/publicdomain/zero/1.0/>. */
|
||||
|
||||
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 <stdio.h>
|
||||
|
||||
/* 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<N; mti++) {
|
||||
mt[mti] = (1812433253UL * (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] &= 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<N-M;kk++) {
|
||||
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
|
||||
mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
|
||||
}
|
||||
for (;kk<N-1;kk++) {
|
||||
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
|
||||
mt[kk] = mt[kk+(M-N)] ^ (y >> 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 <objbase.h>
|
||||
#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<uint64 *> 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;
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
|
|
@ -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])&]
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue