Fastest PRNG around, the XORSHIFT

[ShawnLG]

If you need a fast PRNG. XORSHIFTs are the fastest. Use them for generating random events in games or for Perlin noise for textures or terrain.
xorshift32() has a short random sequence of states. It's only 32bits. You can see this in the random walk test. Wait a few minutes(on a i5 3470) and you will see a fractal like pattern.

Code: Select all

'XorShift PRNG variations
'Translated from https://en.wikipedia.org/wiki/Xorshift
'seeds must be initialize other than zero or they will fail.
type u128longint
    a as ulongint
    b as ulongint
end type
type u128int
    a as ulong
    b as ulong
    c as ulong
    d as ulong
end type
type xorwow128
    a as ulong
    b as ulong
    c as ulong
    d as ulong
    e as ulong 'counter
end type


declare function xorshift32(byref seed as ulong) as ulong
declare function xorshift64(byref seed as ulongint) as ulongint
declare function xorshift64s(byref seed as ulongint) as ulongint
declare function xorshift128(byref seed as u128int) as ulong
declare function xorshift128p(byref seed as u128longint) as ulongint
declare function xorwow128(byref seed as xorwow128) as ulong

dim as ulong s32
dim as ulongint s64
dim as u128int s128
dim as u128longint p128
dim as xorwow128 w128
'initialize the seeds with some random integer (zero will fail). use timer for different seeds when progam starts.
s32 = 6765
s64 = 1456
s128.a = 4567
s128.b = 4635656
s128.c = 3485789
s128.d = 43588923
w128.a = 3405
p128.a = 4395879

const sx = 1280 'screeen resolution
const sy = 720
const zoomfactor = 7' sets the zoom out factor for random walk (2^zoomfactor)

dim as integer x, y
screenres sx, sy, 32, 1

'random point test
do
    x = xorshift32(s32) mod sx
    y = xorshift32(s32) mod sy
    pset (x,y), rgb(255,255,255)
loop until inkey <> ""
cls

'random walk test
x = (sx \ 2) shl zoomfactor 'start at center of screen
y = (sy \ 2) shl zoomfactor
do
    x = x + xorshift32(s32) mod 3 - 1
    y = y + xorshift32(s32) mod 3 - 1
    pset (x shr zoomfactor,y shr zoomfactor), rgb(255,255,255)
loop until inkey <> ""




function xorshift32(byref seed as ulong) as ulong
    seed = seed xor (seed shr 13)
    seed = seed xor (seed shl 17)
    seed = seed xor (seed shr 5)
	return seed
end function

function xorshift64(byref seed as ulongint) as ulongint
    seed = seed xor (seed shr 13)
    seed = seed xor (seed shl 7)
    seed = seed xor (seed shr 17)
    return seed
end function

function xorshift64s(byref seed as ulongint) as ulongint
    seed = seed xor (seed shr 12)
    seed = seed xor (seed shl 25) 
    seed = seed xor (seed shr 27)
    return seed * &h2545F4914F6CDD1D
end function

function xorshift128(byref seed as u128int) as ulong
    ' Algorithm "xor128" from p. 5 of Marsaglia, "Xorshift RNGs" */
	dim as ulong t = seed.d
	dim as ulong s = seed.a
	seed.d = seed.c
	seed.c = seed.b
	seed.b = s
    t = t xor (t shr 11)
    t = t xor (t shl 8)
	seed.a = (t xor s) xor (s shl 19)
    return seed.a
end function

function xorwow128(byref seed as xorwow128) as ulong
    ' Algorithm "xorwow" from p. 5 of Marsaglia, "Xorshift RNGs" */
	dim as ulong t = seed.d
	dim as ulong s = seed.a
	seed.d = seed.c
	seed.c = seed.b
	seed.b = s
    t = t xor (t shl 2)
    t = t xor (t shr 1)
    t = (t xor s) xor (s shr 4)
	seed.a = t
	seed.e = seed.e + 362437
	return t + seed.e
end function

function xorshift128p(byref seed as u128longint) as ulongint
    dim as ulongint t = seed.a
	dim as ulongint s = seed.b
	seed.a = s
    t = t xor (t shl 23)
    t = t xor (t shr 17)
    t = t xor s xor (s shr 26)
    seed.b = t
    return t + s
end function

[D.J.Peters]

You have to replace uInteger with ulong 32 vs 64 bit !

Joshy

[ShawnLG]

Can't do that. Xorshift32() and xorshift64() are optimized for their own bit depth. The original author did not choose the shift parameters for no reason.

[Imortis]

Can't do that. Xorshift32() and xorshift64() are optimized for their own bit depth. The original author did not choose the shift parameters for no reason.

FB's Integer and UInteger types fit the size of the native type of the processor. So on 32bit they are 32bit but on 64bit, they are 64bit. If you want it to be consistant, anything that needs a 32bit memory space should be changed to Long as that is 32bit no matter what the CPU type.

[MrSwiss]

XORshift may be fastest, however, at the cost of randomness (simply put: quality).

For a newer RNG (2014) with better quality at about a equal speed, see:
pcg32rr.bi - minimal PCG32_random_r implementation ( (c) Melissa E. O'Neill )
There are also links to external articles, provided by experts in the field.

[ShawnLG]

According to the FB manual, only long and ulong vary depending on platform. Everything else is static. Maybe the FB manual is wrong?

[MrSwiss]

According to the FB manual, only long and ulong vary depending on platform.

No, Long/ULong are static (32 bit always), U/Integer is the 'dynamic' one
(depending on pointer size).
Where in the manual (is it current?), did you get that from?

See: Integer Types (Standard Data Types page).

[ShawnLG]

Wow. I found that FBIDE I am using still has the old FB-manual-0.23.0.chm file that came with it. I updated it to 1.060 which is the current compiler version I am using.

[deltarho[1859]]

It is worth noting that PCG32 passes the PractRand test suite up to 16TB. Many of the Xorshift algorithms, especially those by Sebastian Vigna, fail PractRand at 64GB. It follows then that if a game, for example, uses less than 32GB then any quality deficiency is below the radar. Mersenne Twister fails PractRand v0.94 at 256GB. 128GB is an awful lot of random number requests. It took v0.94 to spot the weakness- all was well for v0.93 and earlier.

[paul doe]

@David: no mention of the Middle Square Weyl Sequence PRNG? ;)

Minimal implementation here:

Code: Select all

'' Convenience
const _
  ulongMax => culng( -1 ), _
  ulongIntMax => culngint( -1 )

/'
  Middle Square Weyl Sequence Pseudorandom Number Generator
'/
function _
  mswsrnd( _
    byval seed as ulongint => 0 ) _
  as ulong
  
  static as ulongint _
    x, _
    w = &h5851f42d4c957f2d, _
    s = &hb5ad4eceda1ce2a9
  
  if( seed <> 0 ) then
    w => seed
    x => 0
  end if
  
  x *= x : w += s : x += w
  x = ( x shr 32 ) or ( x shl 32 )
  
  return( x )
end function

'' Seed the rng
randomize( timer() )
mswsrnd( rnd() * ulongIntMax )

'' Optional: you can warm it up (it isn't really necessary, but still)
for _
  i as integer => 0 _
  to 100
  
  mswsrnd()
next

'' A simple example usage
for _
  i as integer => 0 _
  to 10
  
  '' To get a value in the 0 <= n < 1 range
  ? ( mswsrnd() / ulongMax )
next

sleep()

[deltarho[1859]]

Hi Paul

I mentioned PCG32 because it was mentioned earlier.

This is my latest version of MsWs. I don't use 's' because of the issues we had with it when discussing MsWs some time ago. Instead, I use the largest prime number < 2^64.

With my latest speed test method PCG32 is faster than MsWs but not by a large margin. If I remember correctly I went to 2TB with Practrand so nobody would get sacked for using MsWs.

Code: Select all

#define Get64Bit (Cast( Ulongint, Rnd*(2^32) ) Shl 32) Or Cast( Ulongint, Rnd*(2^32) )

Type MsWsParams
  Declare Constructor
  Declare Sub MyRandomize( Byval seed As Ulongint = 0, Byval sequence As Ulongint = 0 )
  Declare Function Rand() As Ulong
  Declare Function RandD() As Double
  declare function Range overload ( First as double, Last as double ) as double
  declare function Range overload ( First as longint, Last as longint ) as longint
  declare Function Gauss() as double
  Seed As Ulongint
  Sequence As Ulongint
End Type

Sub MsWsParams.MyRandomize( Byval seed As Ulongint = 0, Byval sequence As Ulongint = 0 )
  Randomize , 5
  If seed = 0 Then
    This.Seed = Get64Bit
    This.Sequence = Get64Bit
  Else
    If sequence = 0 Then
      This.Seed = seed
      This.Sequence = Get64Bit
    Else ' For a sequence repeat
      This.Seed = seed
      This.Sequence = sequence
    End If
  End If
  ' Warm up generator - essential for any PRNG
  For i As Ulong = 1 To 1000
    This.rand()
  Next
End Sub

' 18446744073709551557 is the largest prime number < 2^64

Function MsWsParams.Rand() As Ulong
  This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
  Return( This.Seed )
End Function

Function MsWsParams.RandD() As Double
  Dim As Ulong temp
  This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
  temp = This.Seed
  Return temp/2^32
End Function

Function MsWsParams.Range( First as double, Last as double ) as double
  Dim As Ulong temp
  This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
  temp = This.Seed
  Function = temp/2^32*( Last - First ) + First 
End Function

Function MsWsParams.Range(First As Longint, Last As Longint) As LongInt
  Dim As Ulong temp
  This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
  temp = This.Seed
  Function = temp/2^32 * (Last - First + 1) - .5 + First
End Function

Function MsWsParams.Gauss As double
Static As Long u2_cached
Static As double u1, u2, x1, x2, w
 
  If u2_cached = -1 Then
    u2_cached = 0
    Function = u2
  Else
    Do
      x1 = This.RandD
      x2 = This.RandD
      w = x1 * x1 + x2 * x2
    Loop While w >= 1
    w = Sqr( -2 * Log(w)/w )
    u1 = x1 * w
    u2 = x2 * w
    u2_cached = -1
    Function = u1
  End If
 
End Function

Constructor MsWsParams
  This.MyRandomize()
End Constructor

[paul doe]

...
This is my latest version of MsWs. I don't use 's' because of the issues we had with it when discussing MsWs some time ago. Instead, I use the largest prime number < 2^64.

I use the same above, based on your suggestion.

...
With my latest speed test method PCG32 is faster than MsWs but not by a large margin. If I remember correctly I went to 2TB with Practrand so nobody would get sacked for using MsWs.

I grew bored at 8TB and I think Bernard took it to 16TB. He was also kind enough to update its paper and inform me. Link is the same that was given in your other thread (sorry, don't have time to look for it right now XD)

[deltarho[1859]]

I use the same above, based on your suggestion.

So you do - sorry.

I grew bored at 8TB...

When you get to my age you stop at 2TB.

How long has he been dead?
I don't know, but look he got to 16TB with PractRand.
I'd like to think he actually saw that.
With that smile on his face he probably did.

[deltarho[1859]]

Comparing Paul's code with mine I spotted some faulty logic with mine but the results were still random 'enough' not to bother PractRand. Hmmm, I suppose that could mean my faulty logic was random. Corrected code posted above.

[MrSwiss]

Hmmm, I suppose that could mean my faulty logic was random. Corrected code posted above.

Your faulty logic starts off, with a wrong formula (calculating maximum value of ULong).
It isn't: 2^32 but: 2^32-1

Below code proofes, that your formula 'breaks' a ULong (the result needs 33 bit's):

Code: Select all

Const As ULong  ULmax = &hFFFFFFFF      ' known binary result (in HEX)

'' the below 'number crunching' is totally useless, because the binary
'' result is already known (see const above), regardless of knowing,
'' what it means in decimal (which is a: interpreted value)
Dim As ULongInt wResult = 2 ^ 32        ' blows a ULong (is 33 bit)
Dim As ULong    cResult = 2 ^ 32 - 1    ' max. ULong 'calculated'


Print "wrong formula (ULongInt) 33 bits!"   ' blows the ULong (33 bits)
Print "2 ^ 32 (decimal) = "; wResult
Print "2 ^ 32 (binary)  = "; Bin(wResult, 64)
Print "2 ^ 32 (hex)     = "; Hex(wResult, 16)
Print
Print "correct formula (ULong)"         ' correct ULong maximum
Print "2 ^ 32 - 1 (decimal) = "; cResult
Print "2 ^ 32 - 1 (binary)  = "; Bin(cResult, 32)
Print "2 ^ 32 - 1 (hex)     = "; Hex(cResult, 8)
Print
Print "ULong max. value, for comparison (Const ULong)"  ' dito (as above)
Print "Const (decimal) = "; ULmax
Print "Const (binary)  = "; Bin(ULmax, 32)
Print "Const (hex)     = "; Hex(ULmax, 8)
Print

Sleep

You are also inducing unnecessary 'number crunching' which causes speed loss.
Use a once defined Const instead ...

(Btw. if you feel the need to get angry at someone, taget yourself this time.)