Oh no, not another PRNG

General FreeBASIC programming questions.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Oh no, not another PRNG

Post by dodicat »

Dafhi.
Make up a pipe file (compile one)
Use the command as per deltarho[] or just simply as given below.
Pipefile example

Code: Select all


'pipefilecrt2.exe

'command >>  pipefilecrt2.exe | rng_test stdin32 -multithreaded

'======
'my own generater (skew)
'======================
#include "crt.bi"

randomize 0,1 


#define _r_(l,f) Rnd*(l-f)+f
#define rnd2 log10(_r_(1,9.99999999999999))
#define Irange(f,l) Int(Rnd2*(((l)+1)-(f)))+(f)
'===============


Dim Shared S As String * 1048576
Dim As Ulong Ptr SPtr, BasePtr
Dim As Long j

SPtr = Cptr(Ulong Ptr, StrPtr( S ))
BasePtr = SPtr
 
Do
  For j = 1 to 262144
    *SPtr = Irange(0,65534)'<--- your random function
    SPtr += 1
  Next
  print S;
  SPtr = BasePtr
Loop
  
My skew generator of course is crap and doesn't get far.

Code: Select all

Microsoft Windows [Version 10.0.18362.900]
(c) 2019 Microsoft Corporation. All rights reserved.

C:\Users\User\Desktop\bin\random\PractRand_0.94\PractRand_094\bin\msvc12_64bit>pipefilecrt2.exe | rng_test stdin32 -mult
ithreaded
RNG_test using PractRand version 0.94
RNG = RNG_stdin32, seed = unknown
test set = core, folding = standard (32 bit)

rng=RNG_stdin32, seed=unknown
length= 64 megabytes (2^26 bytes), time= 2.4 seconds
  Test Name                         Raw       Processed     Evaluation
  BCFN(2+0,13-3,T)                  R>+99999  p = 0           FAIL !!!!!!!!
  BCFN(2+1,13-3,T)                  R>+99999  p = 0           FAIL !!!!!!!!
  BCFN(2+2,13-4,T)                  R>+99999  p = 0           FAIL !!!!!!!!
  BCFN(2+3,13-4,T)                  R>+99999  p = 0           FAIL !!!!!!!!
  BCFN(2+4,13-5,T)                  R>+99999  p = 0           FAIL !!!!!!!!
  BCFN(2+5,13-5,T)                  R>+99999  p = 0           FAIL !!!!!!!!
  BCFN(2+6,13-6,T)                  R>+99999  p = 0           FAIL !!!!!!!!
  BCFN(2+7,13-6,T)                  R>+99999  p = 0           FAIL !!!!!!!!
  BCFN(2+8,13-7,T)                  R=+69534  p = 0           FAIL !!!!!!!!
  BCFN(2+9,13-8,T)                  R=+41141  p = 0           FAIL !!!!!!!!
  BCFN(2+10,13-8,T)                 R=+20506  p =  3e-5205    FAIL !!!!!!!!
  BCFN(2+11,13-9,T)                 R=+11761  p =  1e-2643    FAIL !!!!!!!!
  BCFN(2+12,13-9,T)                 R= +5862  p =  2e-1318    FAIL !!!!!!!!
  DC6-9x1Bytes-1                    R>+99999  p = 0           FAIL !!!!!!!!
  Gap-16:A                          R>+99999  p = 0           FAIL !!!!!!!!
  Gap-16:B                          R>+99999  p = 0           FAIL !!!!!!!!
  FPF-14+6/16:(0,14-0)              R=+40620  p = 0           FAIL !!!!!!!!
  FPF-14+6/16:(1,14-0)              R=+48067  p = 0           FAIL !!!!!!!!
  FPF-14+6/16:(2,14-1)              R=+33939  p = 0           FAIL !!!!!!!!
  FPF-14+6/16:(3,14-2)              R=+23971  p = 0           FAIL !!!!!!!!
  FPF-14+6/16:(4,14-2)              R=+23375  p = 0           FAIL !!!!!!!!
  FPF-14+6/16:(5,14-3)              R=+16849  p = 0           FAIL !!!!!!!!
  FPF-14+6/16:(6,14-4)              R=+11662  p =  6e-9529    FAIL !!!!!!!!
  FPF-14+6/16:(7,14-5)              R= +8359  p =  6e-6930    FAIL !!!!!!!!
  FPF-14+6/16:(8,14-5)              R= +6538  p =  8e-5420    FAIL !!!!!!!!
  FPF-14+6/16:(9,14-6)              R= +3938  p =  5e-3014    FAIL !!!!!!!!
  FPF-14+6/16:(10,14-7)             R= +3379  p =  9e-2689    FAIL !!!!!!!!
  FPF-14+6/16:(11,14-8)             R= +1983  p =  7e-1427    FAIL !!!!!!!!
  FPF-14+6/16:(12,14-8)             R= +2000  p =  1e-1439    FAIL !!!!!!!!
  FPF-14+6/16:(13,14-9)             R= +1005  p =  4.0e-633   FAIL !!!!!!!
  FPF-14+6/16:(14,14-10)            R=+709.1  p =  7.7e-378   FAIL !!!!!!!
  FPF-14+6/16:(15,14-11)            R=+797.5  p =  2.0e-348   FAIL !!!!!!!
  FPF-14+6/16:(16,14-11)            R=+546.9  p =  3.1e-239   FAIL !!!!!!
  FPF-14+6/16:all                   R>+99999  p = 0           FAIL !!!!!!!!
  FPF-14+6/16:cross                 R>+99999  p = 0           FAIL !!!!!!!!
  BRank(12):128(4)                  R= +2544  p~=  4e-1354    FAIL !!!!!!!!
  BRank(12):256(2)                  R= +3748  p~=  3e-1129    FAIL !!!!!!!!
  BRank(12):384(1)                  R= +4028  p~=  1e-1213    FAIL !!!!!!!!
  BRank(12):512(2)                  R= +7644  p~=  3e-2302    FAIL !!!!!!!!
  BRank(12):768(1)                  R= +8161  p~=  1e-2457    FAIL !!!!!!!!
  BRank(12):1K(1)                   R=+10916  p~=  3e-3287    FAIL !!!!!!!!
  mod3n(5):(0,9-2)                  R>+99999  p = 0           FAIL !!!!!!!!
  mod3n(5):(1,9-2)                  R>+99999  p = 0           FAIL !!!!!!!!
  TMFn(2+0):wl                      R=+12030  p~=  0          FAIL !!!!!!!!
  TMFn(2+1):wl                      R= +5985  p~=  0          FAIL !!!!!!!!
  [Low8/32]DC6-9x1Bytes-1           R= +24.8  p =  4.4e-14    FAIL !
  [Low8/32]Gap-16:B                 R=  +6.4  p =  1.7e-5   mildly suspicious
  [Low8/32]FPF-14+6/16:(0,14-1)     R=+152.8  p =  3.9e-135   FAIL !!!!!
  [Low8/32]FPF-14+6/16:(1,14-2)     R=+121.2  p =  9.2e-106   FAIL !!!!!
  [Low8/32]FPF-14+6/16:(2,14-2)     R= +60.0  p =  2.8e-52    FAIL !!!!
  [Low8/32]FPF-14+6/16:(3,14-3)     R= +42.4  p =  6.7e-37    FAIL !!!
  [Low8/32]FPF-14+6/16:(4,14-4)     R= +37.8  p =  8.4e-31    FAIL !!!
  [Low8/32]FPF-14+6/16:(5,14-5)     R= +24.9  p =  1.5e-20    FAIL !
  [Low8/32]FPF-14+6/16:all          R=+208.1  p =  8.9e-195   FAIL !!!!!!
  [Low8/32]FPF-14+6/16:cross        R=+130.8  p =  3.6e-103   FAIL !!!!!
  [Low8/32]mod3n(5):(0,9-3)         R= +14.8  p =  8.1e-8   very suspicious
  ...and 84 test result(s) without anomalies

 
You notice that I work directly in the
PractRand_0.94\PractRand_094\bin\msvc12_64bit folder.
I have a startshell.exe file in there
startshell.bas

Code: Select all

  shell "cmd" 
Note deltarho wrote the pilefile format a while back.
dafhi
Posts: 1640
Joined: Jun 04, 2005 9:51

Re: Oh no, not another PRNG

Post by dafhi »

thanks dodicat. x128xx (xorshiro128**) mylit as ULong passes to 32 Gigs. Total Fail @ UShort
I suspect it has something to do with multiplier granularity

my generator-in-development doesn't trigger as many fails at 16 bits, but doesn't succeed @ 32

Code: Select all

#include "crt.bi"

randomize 0,1


type mylit as ulong

const bit_size = len(mylit) * 8

#define rotl(in, amount)  (in shl amount or in shr (bit_size - amount))

function mul_conv(m as integer) as integer
  return (bit_size * m + 32) \ 64 '' add "0.5" then truncate divide (FB)
end function

function x128xx as mylit
  static As mylit s0, s1, result
 
  static as mylit s(1) = {1,0} '' seed
 
  s0 = s(0)
  s1 = s(1)
  result = rotl(s0 * mul_conv(5), mul_conv(7)) * mul_conv(9)

  s1 xor= s0
  s(0) = rotl(s0, mul_conv(24)) xor s1 xor (s1 shl mul_conv(16)) ' a, b
  s(1) = rotl(s1, mul_conv(37)) ' c

  return result
end function


function csg_hp_2 as mylit
  static as mylit  a
  static as mylit  w
 
  const rot_amount = bit_size / 2 - 1
 
  a *= 3
  a xor= w
  a = rotl(a, rot_amount)
  w += 1
  return a
End Function


#define _r_(l,f) Rnd*(l-f)+f
#define rnd2 log10(_r_(1,9.99999999999999))
#define Irange(f,l) Int(Rnd2*(((l)+1)-(f)))+(f)
'===============


Dim Shared S As String * 1048576
Dim As mylit Ptr SPtr, BasePtr
Dim As Long j

SPtr = Cptr(Ulong Ptr, StrPtr( S ))
BasePtr = SPtr
 
Do
  For j = 1 to 262144 * 4 \ len(mylit)
    *SPtr = x128xx
    '*SPtr = csg_hp_2
    SPtr += 1
  Next
  print S;
  SPtr = BasePtr
Loop
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Oh no, not another PRNG

Post by deltarho[1859] »

xoshiro256** Image

I had to have a look.

The most important question is will it get to 1TB with PractRand? It did with only two small anomalies on the way. So, after 'tanking' with 64GB failures a couple of years ago Vigna and Blackman have come up trumps with these two gems.

This has a period of 2^256. It has a state vector of four Ulongints, twice as large as xoroshiro128**. The engine is a little 'heavier' with much more activity with the state vector and this is reflected in being about 22% slower than xoroshiro128**. Further, tests on PCG32II and xoroshiro128** show that they are actually pretty much 'neck and neck' on speed in 32-bit mode. In 64-bit mode the Vigna generators are close in speed, with xoshiro256** being a little faster, and both are over 21% faster than PCG32II. Of course, PCG32II's granularity is 32-bit whereas the two Vigna generators have a granularity of 53-bit. On my machine, i7 3.9GHz, seven years old now, the two Vigna generators are GHz generators. I thought that I would have to wait until my next PC upgrade before I saw that.

So what is the argument for using the 256 as opposed to the 128. I don't have one and will use the 128.

From a parallelization perspective the 128 could be configured to have a period of 2^64 and 2^64 sequences. Someone who favours 2^128 generators would prefer the 256 because that could be configured to have a period of 2^128 and 2^128 sequences.

I wrote an edited version of the usage program and it behaved in the same way as with xoroshiro128**; that is it did what it was supposed to do.

One thing I have noticed when running PractRand is that on getting to 1TB I close the command window but the exe is still running forcing me to 'End task' in Task Manager. I spotted that when I realized that the CPU Usage was at 13%, indicating a program 'hang' on my machine, and the CPU temperature was higher than normal.

There is an interesting aspect to periods.

Suppose we wrote some code to do nothing else but generate random floats and ran it for 28 days. Remember earlier I likened the sequence to a circle. Well, at the end of the 28 days we will, obviously, have used a percentage of the circumference. If we repeated the exercise the chance of starting within the arc covered with the first run would be 3.6 x 10^30 to 1 against with the Vigna generators. With PCG32II the chance would be about 3.0 x 10^9 to 1 against; highly unlikely. With the Vigna generators we can say "It ain't going to happen". That is in 32-bit mode. With 64-bit mode we have 1.5 x 10^28 to 1 against for the Vigna generators and 9.8 x 10^8 to 1 against for PCG32II; which is still highly unlikely. [Added: PCG32II's likelihood will be greatly reduced if we assume random sequences and more in line with the Vigna generators.] With PowerBASIC's generator its period, 2^32, would be expired in less than one hour so another run would guarantee an overlap. With PractRand 4GB, and we would be done, not that it can get to 4GB because it is a LCG and fails at 4KB. PowerBASIC's generator should have been killed years ago. However, to replace it with anything half decent would require a 64-bit unsigned integer which was on the 'ToDo list' in 2001 and never got implemented.

PCG32II has not been put into the long grass. Having another generator in a secondary thread of execution is easy to set up. All we do is create another instance and use a different sequence number to the instance in the primary thread; giving no collisions and different sequences. Do I have any applications using that? Not one, but I am prepared if ever I do. PowerBASIC's generator is not thread safe so only one instance is allowed.

This is my seventh generator, and I am still loving it. There is no accounting for taste, is there?

So, there you have, yet another PRNG.

xoshiro256.bas

Code: Select all

' http://prng.di.unimi.it/xoshiro256starstar.c
 
/'
static inline uint64_t rotl(const uint64_t x, int k) {
	return (x << k) | (x >> (64 - k));
}
 
static uint64_t s[4];
 
uint64_t next(void) {
	const uint64_t result = rotl(s[1] * 5, 7) * 9;
 
	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] = rotl(s[3], 45);
 
	return result;
}
'/
 
#define Make64Bit ( Cast( Ulongint, Rnd*(2^32) ) Shl 32) Or Cast( Ulongint, Rnd*(2^32) )
#define rotl(x,k) ( (x Shl k) Or (x Shr(64-k)) ) ' Note the extra parentheses
 
#macro Engine
  result = rotl(This.s(1)*5,7)*9
  t = This.s(1) shl 17
  This.s(2) Xor= This.s(0)
  This.s(3) Xor= This.s(1)
  This.s(1) Xor= This.s(2)
  This.s(0) Xor= This.s(3)
 
  This.s(2) Xor= t
 
  This.s(3) =rotl(This.s(3),45)
#endmacro
 
' Generator is Sebastian Vigna's xoshiro256**
 
Type xoshiro256
  Public:
  Declare Constructor
  Declare Sub MyRandomize( seed as string )
  Declare Function rand() As Ulongint
  Declare Function randD() As Double
  Declare Function range overload( Byval One As Long, Byval Two As Long ) As Long
  Declare Function range overload ( byval One as double, Byval Two as double ) as double
  Declare Sub GetSnapshot()
  Declare Sub SetSnapshot()
  Private:
  As Ulongint s(0 To 3)
  As Ulongint snapshot(0 to 3)
End Type
 
Private Function Get64Bit() As Ulongint
  Dim As Ulongint dummy = Make64Bit
  Return rotl(dummy, 16)
End Function
 
Private Sub xoshiro256.MyRandomize( seed as string )
Dim as ulong i
  ' If seed = "" then use cryptographic random numbers else use Mersenne Twister
  If seed = "" Then Randomize , 5 Else randomize val(seed), 3
  This.s(0) = Get64Bit
  This.s(1) = Get64Bit
  ' Warm up generator - essential
  For i = 1 To 8192
    this.rand()
  Next i
End Sub
 
Private Function xoshiro256.rand() As ulongint
  Dim As Ulongint t, result
  Engine
  Return result
End Function
 
Private Function xoshiro256.randD() As Double
Dim As Ulongint t, result
  Engine
  Return result/2^64
End Function
 
Private Function xoshiro256.range( Byval One As Long, Byval Two As Long ) As Long
  Dim As Ulongint t, result
  Engine
  Return Clng(result Mod ( Two-One+1 )) + One ' By dodicat
End Function
 
 Private Function xoshiro256.range( Byval One As Double, Byval Two As Double ) As Double
   Dim As Ulongint t, result
   Engine
   Return result/2^64*( Two - One ) + One
 End Function
 
 Private Sub xoshiro256.GetSnapshot()
  This.snapshot(0) = This.s(0)
  This.snapshot(1) = This.s(1)
  This.snapshot(2) = This.s(2)
  This.snapshot(3) = This.s(3)
end sub
 
Private Sub xoshiro256.SetSnapshot()
  This.s(0) = This.snapshot(0)
  This.s(1) = This.snapshot(1)
  This.s(2) = This.snapshot(2)
  This.s(3) = This.snapshot(3)
end sub
 
Constructor xoshiro256
  This.MyRandomize( "" )
  This.GetSnapshot
End Constructor
 
Dim Shared As xoshiro256 x256
#undef Rnd
#define Rnd x256.randD
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: Oh no, not another PRNG

Post by srvaldez »

Hi deltarho[1859]
just for fun with the help of fbfrog I translated xoroshiro1024starstar.c, I hope you don't mind me posting it here, if you do object then I will delete the code, it's probably not good for any practical use anyway.

Code: Select all

private function rotl(byval x as const ulongint, byval k as long) as ulongint
	return (x shl k) or (x shr (64 - k))
end function

dim shared p as long
dim shared s(0 to 15) as ulongint

private function next_() as ulongint
	dim q as const long = p
	p = (p + 1) and 15
	dim as const ulongint s0 = s(p)
	dim s15 as ulongint = s(q)
	dim result as const ulongint = rotl(s0 * 5, 7) * 9
	s15 xor= s0
	s(q) = (rotl(s0, 25) xor s15) xor (s15 shl 27)
	s(p) = rotl(s15, 36)
	return result
end function

private sub jump()
	static JUMP_(0 to ...) as const ulongint = {&h931197d8e3177f17, &hb59422e0b9138c5f, &hf06a6afb49d668bb, &hacb8a6412c8a1401, &h12304ec85f0b3468, &hb7dfe7079209891e, &h405b7eec77d9eb14, &h34ead68280c44e4a, &he0e4ba3e0ac9e366, &h8f46eda8348905b7, &h328bf4dbad90d6ff, &hc8fd6fb31c9effc3, &he899d452d4b67652, &h45f387286ade3205, &h03864f454a8920bd, &ha68fa28725b1b384}
	dim t(0 to ubound(s)) as ulongint
	'memset(t, 0, sizeof t);
	' not needed in FB
	for i as long = 0 to ubound(JUMP_)
		for b as long = 0 to 63
			if (JUMP_(i) and 1ull shl b) then
				for j as long = 0 to ubound(s)
					t(j) xor= s((j + p) and ubound(s))
				next
			end if
			next_()
		next
	next

	for i as long = 0 to ubound(s)
		s((i + p) and ubound(s)) = t(i)
	next
end sub

private sub long_jump()
	static LONG_JUMP_(0 to ...) as const ulongint = {&h7374156360bbf00f, &h4630c2efa3b3c1f6, &h6654183a892786b1, &h94f7bfcbfb0f1661, &h27d8243d3d13eb2d, &h9701730f3dfb300f, &h2f293baae6f604ad, &ha661831cb60cd8b6, &h68280c77d9fe008c, &h50554160f5ba9459, &h2fc20b17ec7b2a9a, &h49189bbdc8ec9f8f, &h92a65bca41852cc1, &hf46820dd0509c12a, &h52b00c35fbf92185, &h1e5b3b7f589e03c1}
	dim t(0 to ubound(s)) as ulongint
	'memset(t, 0, sizeof t);
	' not needed in FB
	for i as long = 0 to ubound(LONG_JUMP_)
		for b as long = 0 to 63
			if (LONG_JUMP_(i) and 1ull shl b) then
				for j as long = 0 to ubound(s)
					t(j) xor= s((j + p) and ubound(s))
				next
			end if
			next_()
		next
	next

	for i as long = 0 to ubound(s)
		s((i + p) and ubound(s)) = t(i)
	next
end sub

for i as long =0 to 15
	s(i)=i+1
next
p=3
jump()
for i as long =0 to 15
	? next_()
next
Last edited by srvaldez on Jul 12, 2020 16:25, edited 1 time in total.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Oh no, not another PRNG

Post by deltarho[1859] »

@srvaldez

Why should I object? It doesn't contradict the thread's title.

I didn't go that far because even Vigna wrote: "Its state however is too large--in general, the xoshiro256 family should be preferred."

My default is now xoroshiro128** but as mentioned above PCG32II has not been put into the long grass.
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: Oh no, not another PRNG

Post by paul doe »

Xoshiro256 aka 'Vigna's Revenge' XD

Quite nice. Here's a version that caters to my own tastes:

Code: Select all

type _
  Xoshiro256Snapshot
  
  declare constructor( _
    byval as ulongint, _
    byval as ulongint, _
    byval as ulongint, _
    byval as ulongint )
  
  as ulongint _
    s0, s1, s2, s3
end type

constructor _
  Xoshiro256Snapshot( _
    byval aS0 as ulongint, _
    byval aS1 as ulongint, _
    byval aS2 as ulongint, _
    byval aS3 as ulongint )
  
  s0 => aS0 : s1 => aS1 : s2 => aS2 : s3 => aS3
end constructor

type _
  Xoshiro256
  
  public:
    declare constructor()
    declare destructor()
    
    declare function _
      randomized( _
        byval as double => 0.0 ) _
      byref as Xoshiro256
    declare function _
      rand() as ulongint
    declare function _
      randD() as double
    declare function _
      range( _
        byval as long, _
        byval as long ) _
      as long
    declare function _
      range( _
        byval as double, _
        byval as double ) _
      as double
    declare function _
      getSnapshot() as Xoshiro256Snapshot
    declare function _
      setSnapshot( _
        byref as const Xoshiro256Snapshot ) _
      byref as Xoshiro256
    
  private:
    declare function _
      get64Bit() as ulongint
    
    static as const ulongint _
      LONG_MAX
    static as const ulongint _
      LONGINT_MAX
    static as const double _
      INV_LONGINT_MAX
    
    as ulongint _
      _s( 0 to 3 )
    
  #define make64Bit _
    ( culngint( rnd() * LONG_MAX ) shl 32 ) or _
      culngint( rnd() * LONG_MAX )
  #define rotl( x, k ) _
    ( ( x shl k ) or ( x shr ( 64 - k ) ) )
  
  #macro shuffle( _r_ )
    _r_ => rotl( _s( 1 ) * 5, 7 ) * 9
    
    dim as ulongint _
      t => _s( 1 ) shl 17
    
    _s( 2 ) xor=> _s( 0 )
    _s( 3 ) xor=> _s( 1 )
    _s( 1 ) xor=> _s( 2 )
    _s( 0 ) xor=> _s( 3 )
   
    _s( 2 ) xor=> t
   
    _s( 3 ) => rotl( _s( 3 ), 45 )
  #endmacro
end type

dim as const ulongint _
  Xoshiro256.LONG_MAX => 2 ^ 32
dim as const ulongint _
  Xoshiro256.LONGINT_MAX => 2 ^ 64
dim as const double _
  Xoshiro256.INV_LONGINT_MAX => 1.0 / ( 2 ^ 64 )

constructor _
  Xoshiro256()
  
  randomized()
end constructor

destructor _
  Xoshiro256()
end destructor

function _
  Xoshiro256.get64bit() as ulongint
  
  dim as ulongint _
    n => make64bit
  
  return( rotl( n, 16 ) )
end function

function _
  Xoshiro256.randomized( _
    byval seed as double => 0.0 ) _
  byref as Xoshiro256
  
  if( seed ) then
    randomize( seed, 3 )
  else _
    randomize( , 5 )
  end if
  
  _s( 0 ) => get64bit()
  _s( 1 ) => get64bit()
  _s( 2 ) => get64bit()
  _s( 3 ) => get64bit()
  
  '' Warm up generator - essential
  for i as integer => 1 to 8192
    rand()
  next
  
  return( this )
end function

function _
  Xoshiro256.rand() as ulongint
  
  dim as ulongint _
    result
  
  shuffle( result )
  
  return( result )
end function

function _
  Xoshiro256.randD() as double
  
  dim as ulongint _
    result
  
  shuffle( result )
  
  return( result * INV_LONGINT_MAX )
end function

function _
  Xoshiro256.range( _
    byval aMin as long, _
    byval aMax as long ) _
  as long
  
  dim as ulongint _
    result
  
  shuffle( result )
  
  return( clng( result mod ( aMax - aMin + 1 ) ) + aMin )
end function

function _
  Xoshiro256.range( _
    byval aMin as double, _
    byval aMax as double ) _
  as double
  
  dim as ulongint _
    result
  
  shuffle( result )
  
  return( result * INV_LONGINT_MAX * ( aMax - aMin ) + aMin )
end function

function _
  Xoshiro256.getSnapshot() as Xoshiro256Snapshot
  
  return( Xoshiro256Snapshot( _
    _s( 0 ), _s( 1 ), _s( 2 ), _s( 3 ) ) )
end function

function _
  Xoshiro256.setSnapshot( _
    byref ss as const Xoshiro256Snapshot ) _
  byref as Xoshiro256
  
  _s( 0 ) => ss.s0 : _s( 1 ) => ss.s1
  _s( 2 ) => ss.s2 : _s( 3 ) => ss.s3
  
  return( this )
end function
I like the idea of 'snapshots'. This way, you can persist as many states as you want, and restore them at will. I can see the applications for procedurally generated content. Nice work!
Last edited by paul doe on Jul 13, 2020 1:31, edited 2 times in total.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Oh no, not another PRNG

Post by deltarho[1859] »

I don't know about gas but gcc will convert '1/2^64' to a multiplication for us. Similarly, if we wrote 'x^2' gcc would change that to 'x*x'. I should imagine that a lot of 'tricks' we used in the old days are now done for us. I wish that I knew what they all were because I'd write code faster. Image
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Oh no, not another PRNG

Post by deltarho[1859] »

As is Snapshots are saved to memory - they could be saved to a file. You could be watching a simulation evolve on screen, save the state, power off, go to bed, get up, load the state and continue where you left off.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Oh no, not another PRNG

Post by deltarho[1859] »

Benchmark tests are useful but if generator A is faster than generator B it does not follow that will be true for a particular application.

On running some Monte Carlo code I found that the relationship between 128 and 256 mirrored what I got with the benchmark tests. However, on this occasion PCG32II beat both the Vigna generators on speed. If we wanted 53-bit granularity then we would have to call PCG32II twice and the two Vigna generators would push PCG32II to one side. With regard accuracy of the results there was nothing in it so a 32-bit granularity was good enough allowing PCG32II to take first place.

We have a similar situation with gcc compilers. I use 8.3 as my default as that has the greatest likelihood of producing the fastest binaries. However, if I want the nth degree on speed I have to give 7.5, 9.3 and 10.1 a chance as they may beat 8.3.

I will default to xoroshiro128** but will give PCG32II a try if I am desperate for speed and 32-bit granularity is god enough.

Fortunately, many applications are perfectly happy to be given gcc 5.2 and Mersenne Twister and do not scream "I need more speed". Of course, if they do scream "I need more speed" then the more options we have to achieve that the better. The first port of call, of course, should be to examine the program design as very often a change in design can give us a significant performance boost.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Oh no, not another PRNG

Post by deltarho[1859] »

@paul doe

I would recommend executing 'randomized' in the Constructor because if we forget to initialize the state vector it will be zero everywhere and our random number sequence will be nothing but zeros. We can, of course, override the state vector by calling 'randomized' with a fixed seed, for example, before requesting any numbers.

I also noticed that you are only initializing _s( 0 ) and _s( 1 ) in 'randomized'. That is my fault, I do the same, and now corrected. I ported the function from 128 which only has two Ulongints for the state vector. It is not the end of the world because the warm up code will change _s( 2 ) and _s( 3 ). It is, however, better to initialize all four of the _s().
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: Oh no, not another PRNG

Post by paul doe »

I see. Well, I'll update the code then (I just refactored as it was). Thanks!
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Oh no, not another PRNG

Post by deltarho[1859] »

@paul doe

Code: Select all

constructor _
  Xoshiro256()
 
  randomized()
end constructor
That will use Mersenne Twister(MT) with a zero seed. Most folk will want a random default seed.

Come to think of it I think that you need to look at 'randomized'. As is a positive seed will use cryptographic seeding and zero will use MT with a zero seed. I think the other way around would be better: Zero for cryptographic seeding and a positive seed will use whatever that value is for MT.

My method is a bit 'quaint': An empty string for crypto and Val(<string>) for MT.
paul doe
Moderator
Posts: 1730
Joined: Jul 25, 2017 17:22
Location: Argentina

Re: Oh no, not another PRNG

Post by paul doe »

Indeed, I coded it the other way around XD

There.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Oh no, not another PRNG

Post by dodicat »

Looks like the rnd function is pre processed.
I try to capture rnd 1, 5 (Win32 Crypto API, Linux /dev/urandom) by creating a pointer to the function.
I know that rnd 5 cannot be re seeded, so I expect the totals to be different, but suddenly my pointer seems to point to the twister without my permission.

Code: Select all



randomize 1,5
dim as function(as single=0) as double fn= @rnd()

dim as double t=timer,value=0

for n as long=1 to 100000
   value+= fn(n)
next
print timer-t,value

'===============================
value=0
randomize 1,3

 t=timer
for n as long=1 to 100000
   value+= fn(n)
next
print timer-t,value
sleep

 
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Oh no, not another PRNG

Post by deltarho[1859] »

@dodicat

I have just run your code and got this:

Code: Select all

0.424914799998362           50080.62991477177
 0.001207800009552784        49971.21065035835
How is that confirming Rnd is pre-processed or "my pointer seems to point to the twister without my permission."?

PS: Is this related to the thread's title or should it be in a new thread?
Post Reply