Single and Range wrappers for minimal PCG PRNG

General FreeBASIC programming questions.
deltarho[1859]
Posts: 4313
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

MyRandomize is fairly fast - comes in at about 5 microseconds. I needed to change just the sequence but 5 microseconds was not fast enough.

I added a procedure: Declare Sub SetSequence( Byval NewSequence As Ulongint ). If you use an even number it will be made odd - Or 1.

That is so fast it is coming in below the resolution of my Performance Counter at 10MHz - 100 nanoseconds.

Whilst I was at it I added 'Declare Sub SetState( Byval NewState As Ulongint ).

PCG32II.bas latest version (23 Nov 2020)

Code: Select all

'#Console On
Type pcg32
  Public:
  declare Constructor
  Declare Sub MyRandomize( byval seed As Ulongint = 0, byval sequence As Ulongint = 0 )
  Declare Function rand() As Ulong
	Declare Function randse() As Double
  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 function GetState() as ulongint
  Declare Sub SetState( Byval NewState As Ulongint )
  Declare Function GetSeed() As Ulongint
  Declare Function GetSequence() As Ulongint
  Declare Sub SetSequence( Byval NewSequence As Ulongint )
	declare sub GetSnapshot()
  Declare Sub SetSnapshot()
  declare Function Gauss() as double
  Private:
  state As Ulongint
  seed  As Ulongint
  sequence As Ulongint
  stateSnapshot As ULongInt
  sequenceSnapshot As ulongint
End Type

#ifdef __FB_WIN32__
  #Include Once "windows.bi"
  #Inclib "bcrypt"
  #Include Once "win/wincrypt.bi"
    
  Private Function Get64Bit() As Ulongint
  Dim As BCRYPT_ALG_HANDLE Ptr hRand, hAESAlg
  Dim As BCRYPT_KEY_HANDLE Ptr hKey
  Dim As String IV, AESKey, sRandomBytes, sReturn
  Dim As dword dwResult
    
    BCryptOpenAlgorithmProvider(@hRand, BCRYPT_RNG_ALGORITHM, 0, 0) ' Prepare For Random number generation
    IV = String(16,0) ' Generate a Random IV
    ' Generate a random IV
    BCryptGenRandom(hRand, Strptr( IV ), 16, 0)
     AESKey = String( 32, 0 ) ' Generate a Random encryption key
    ' Generate a random AES key
    BCryptGenRandom(hRand, Strptr( AESKey ), 32, 0) ' 256 bits
    sRandomBytes = String( 4096, 0 )
    ' Generate 256 * AES block size random numbers ' 4KB and thought to be the point at which the entropy pool is refreshed
    BCryptGenRandom(hRand, Strptr( sRandomBytes ), 4096, 0)
    BCryptCloseAlgorithmProvider(hRand, 0)
    BCryptOpenAlgorithmProvider( @hAESAlg, BCRYPT_AES_ALGORITHM, 0, 0 )
    ' Use default CBC algorithm. 6th parameter of 16 indicates AES128
    BCryptGenerateSymmetricKey( hAESAlg, @hKey, 0, 0, Strptr( AESKey ), 16, 0 ) ' We want hKey
    ' Perform the encryption in place; Output buffer = Input buffer
    ' Note that we do not need padding so last paramter is zero. 6th parameter of 16 is AES block size
    BCryptEnCrypt( hKey, Strptr( sRandomBytes ), Len( sRandomBytes ), 0, Strptr( IV ), _
      16, Strptr( sRandomBytes ), Len( sRandomBytes ), @dwResult, 0 )
    BCryptCloseAlgorithmProvider( hAESAlg, 0 )
    BCryptDestroyKey( hKey )
    sReturn += Right( sRandomBytes, 8 ) ' Half of Last block of ciphertext - 64 bits
    Function = Cast(Ulongint, Cvlongint(sReturn))
    If hRand <> 0 Then BCryptCloseAlgorithmProvider( hRand, 0 )
    If hAESAlg <> 0 Then BCryptCloseAlgorithmProvider( hAESAlg, 0 )
    If hKey <> 0 Then BCryptDestroyKey( hKey )
  End Function
#else
  Private Function Get64Bit() As UlongInt ' Will use 'Randomize , 5' by virtue of 'Sub pcg32.MyRandomize'
    return (Cast( Ulongint, Rnd*(2^32) ) Shl 32) Or Cast( Ulongint, Rnd*(2^32) )
  End Function
#endif

Private Function pcg32.rand() As Ulong
  Dim As Ulongint oldstate = this.state
  this.state = oldstate * 6364136223846793005ULL + this.sequence
  Dim As Ulong xorshifted = ((oldstate Shr 18u) xor oldstate) Shr 27u
  Dim As Ulong rot = oldstate Shr 59u
  Return (xorshifted Shr rot) Or (xorshifted Shl ((-rot) And 31))
End Function

Private Sub pcg32.MyRandomize( byval seed As Ulongint = 0, byval sequence As Ulongint = 0 )

  Randomize , 5
  If seed = 0 Then
    If sequence = 0 Then
      this.state = Get64Bit
      this.sequence = Get64Bit or 1
    Else
      this.state = Get64Bit
      this.sequence = sequence or 1
    End If
  Else
    If sequence = 0 Then
      this.state = seed
      this.sequence = Get64Bit or 1
    Else
      this.state = seed
      this.sequence = sequence or 1
    End If
  End If
  This.seed = This.state ' Save initial state
  ' Warm up generator
  For i As Ulong = 1 To 4096
    this.rand()
  Next
    
End Sub

Private Function pcg32.randse() As Double ' Uses 1 x 32-bit random number giving granularity of 32 bits
  dim TempVar As Ulong
  Dim As Ulongint oldstate = this.state
  this.state = oldstate * 6364136223846793005ULL + this.sequence
  Dim As Ulong xorshifted = ((oldstate Shr 18u) xor oldstate) Shr 27u
  Dim As Ulong rot = oldstate Shr 59u
  TempVar = (xorshifted Shr rot) Or (xorshifted Shl ((-rot) And 31))
  Return TempVar/2^32
End Function

Private Function pcg32.randd() As Double ' Uses 2 x 32-bit random numbers giving granularity of 53 bits
  Dim as ulong TempVar1, TempVar2
  Dim As Ulongint oldstate = this.state
  this.state = oldstate * 6364136223846793005ULL + this.sequence
  Dim As Ulong xorshifted = ((oldstate Shr 18u) xor oldstate) Shr 27u
  Dim As Ulong rot = oldstate Shr 59u
  TempVar1 =  (xorshifted Shr rot) Or (xorshifted Shl ((-rot) And 31))
  oldstate = this.state
  this.state = oldstate * 6364136223846793005ULL + this.sequence
  xorshifted = ((oldstate Shr 18u) xor oldstate) Shr 27u
  rot = oldstate Shr 59u
  TempVar2 =  (xorshifted Shr rot) Or (xorshifted Shl ((-rot) And 31))
  Asm
    mov eax, dword ptr [TempVar1]
    movd xmm0, eax
    mov eax, dword ptr [TempVar2]
    movd xmm1, eax
    punpckldq xmm0, xmm1
    psrlq xmm0, 12
    mov eax, 1
    cvtsi2sd xmm1, eax
    por xmm0, xmm1
    subsd xmm0, xmm1
    movq [Function], xmm0
  end asm
end Function

Private Function pcg32.range( Byval One As Long, Byval Two As Long ) As Long 
  Dim As ulong TempVar
  Dim As Ulongint oldstate = this.state
  this.state = oldstate * 6364136223846793005ULL + this.sequence
  Dim As Ulong xorshifted = ((oldstate Shr 18u) xor oldstate) Shr 27u
  Dim As Ulong rot = oldstate Shr 59u
  Tempvar =  (xorshifted Shr rot) Or (xorshifted Shl ((-rot) And 31))
  return clng(TempVar Mod (Two-One+1)) + One ' By dodicat
    
End Function

Private Function pcg32.range( Byval One As Double, Byval Two As Double ) As Double
  Dim TempVar As Ulong
  Dim As Ulongint oldstate = this.state
  this.state = oldstate * 6364136223846793005ULL + this.sequence
  Dim As Ulong xorshifted = ((oldstate Shr 18u) xor oldstate) Shr 27u
  Dim As Ulong rot = oldstate Shr 59u
  TempVar = (xorshifted Shr rot) Or (xorshifted Shl ((-rot) And 31))
  Return TempVar/2^32*( Two - One ) + One
End Function

Private Function pcg32.GetState() As Ulongint
	Return This.state
End Function

Private Sub pcg32.SetState( NewState As UlongInt )
  This.state = NewState
End Sub

Private Function pcg32.GetSeed() As Ulongint
	Return This.seed
End Function

Private Function pcg32.GetSequence() As Ulongint
	Return This.sequence
End Function

Private Sub pcg32.SetSequence( NewSequence as UlongInt )
  This.sequence = NewSequence Or 1
End Sub

Private sub pcg32.GetSnapshot()
  This.stateSnapshot = This.state
  This.sequenceSnapshot = This.sequence
end sub

Private sub pcg32.SetSnapshot()
  This.state = This.stateSnapshot
  This.sequence = This.sequenceSnapshot
end sub

constructor pcg32
  This.MyRandomize
  This.GetSnapshot
end constructor

Private Function pcg32.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.randse
      x2 = This.randse
      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

dim Shared as pcg32 pcg
#undef Rnd
#define Rnd pcg.randse
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by jj2007 »

A bit OT (but you will like it): Better default RNG in the future?
deltarho[1859]
Posts: 4313
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

Thanks Jochen.

It was good to see so many writing about random number generators and the consensus was pretty much where I am. There are many fans of PCG. I also see that xoshiro256** is getting a thumbs up. I wrote this in the Squares thread: "xoroshiro128** and xoshiro256** are worthy of note having periods of 2^128 and 2^256 respectively and 53-bit granularity." They are excellent in 64-bit.

I also noticed someone had settled on gcc 8.3 and -O2 which I have found to be the best combination for FreeBASIC and is my default after much reading and testing. The FB package has an 8.1 option. That should be updated to 8.3.
Post Reply