standard libraries

[deltarho[1859]]

And a bit slower with 32 bits.

pcg32.randse() leaves FB's Rnd #2 or #3 standing in 32-bit mode.

Show me your benchmarking code, sir.

[dodicat]

I used this.

Code: Select all

#cmdline "-Wc -O3"
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 Function GetSeed() As Ulongint
  Declare Function Getsequence() 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"
    
  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
  Function Get64Bit() As UlongInt
    return (Cast( Ulongint, Rnd*(2^32) ) Shl 32) Or Cast( Ulongint, Rnd*(2^32) )
  End Function
#endif

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

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 500
    this.rand()
  Next
  
End Sub

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

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

Function pcg32.range( Byval One As Long, Byval Two As Long ) As Long 
  Dim As ulong TempVar
  dim as long map 
  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))
'  Asm
'    mov edx, dword Ptr [TempVar]
'    mov ecx, [One]
'    mov eax, [Two]
'    cmp ecx, eax
'    jl 0f
'    xchg eax, ecx
'  0:
'    Sub eax, ecx
'    inc eax
'    jz 1f
'    mul edx
'    Add edx, ecx
'  1:
'    mov [Function], edx
'  End Asm
' map = Abs(One) + 5
' One = One + map
' Two = Two + map
 return clng(TempVar Mod (Two-One+1)) + One
  
End Function

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

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

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

Function pcg32.Getsequence() As Ulongint
	Return This.sequence
End Function

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

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

constructor pcg32
  This.MyRandomize
  This.GetSnapshot
end constructor

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 as  pcg32 x
randomize 2,0

dim as double d
for k as long=1 to 10
dim as long lim=10000
var t=timer
for n as long=1 to lim
      d=x.randse
next
print timer-t,"randse"

t=timer
for n as long=1 to lim
      d=rnd
next
print timer-t,"rnd"
print
next k
sleep


 

and my 32 bit result

Code: Select all

 0.000579399994037999       randse
 0.0001568000096199285      rnd

 0.0005907999960896859      randse
 0.0001762999999073145      rnd

 0.0005170000045069401      randse
 0.0001572999954930765      rnd

 0.0004809999877863902      randse
 0.0001555999860585189      rnd

 0.0004788000093327582      randse
 0.0001553999964158948      rnd

 0.0005156000019752582      randse
 0.0001562000066286373      rnd

 0.0005314000086400483      randse
 0.0001569000129535425      rnd

 0.0004792000129469898      randse
 0.0002046999862699295      rnd

 0.0004836999940494025      randse
 0.0002048999999431089      rnd

 0.0005509000137493558      randse
 0.0001558000065813303      rnd

 

[deltarho[1859]]

@dodicat

I have been at this all morning, and it was doing my head in, so I went for a walk in my local park. I was only there for a few minutes and the penny dropped.

#cmdline "-Wc -O3"

You haven't specified a backend code emitter, so gas will be used.

Use #cmdline "-gen gcc -Wc -O3

[dodicat]

I am just back from the shore with my dog.
But I didn't twig, too busy chucking twigs into the water for the dog to fetch.
But you are correct, my apologies, I forgot the -gen gcc.
I wonder why I didn't get a compiler warning saying you can use -Wc -O3 till you are blue in the face, but it won't make a hoot of a difference because you are on the gas backend, or something along those lines.
It would have saved you doing your head in and my acute embarrassment.

[deltarho[1859]]

@dodicat

Good point about the lack of some comment when using -Wc and gas.

Now all you have to do is confirm that PCG32II's randse is faster than FB's #2 and #3 generators when using gcc; even if it is embarrassing.

By the way, -O3 almost invariably generates larger binaries, and every so often we get a performance hit compared with -O2. If -O3 is faster, it is usually negligible. My default is -O2 and srvaldez favours -O2. Most pundits on the internet favour -O2. With your code above -O2 is a little faster on average.

I look at -O3, but it gets rejected most of the time.

[deltarho[1859]]

I am off-topic now, but my latest effort, MsWsII, has a faster randse than PCG32II ( RIP )

[dodicat]

I haven't compiled the fb compiler for years.
I wonder what optimisations are given to compile the rtlib .c functions. (math_rnd.c included)

[srvaldez]

hi dodicat
from the FB makefile

FBC := fbc
CFLAGS := -Wfatal-errors -O2
# Avoid gcc exception handling bloat
CFLAGS += -fno-exceptions -fno-unwind-tables -fno-asynchronous-unwind-tables

I didn't find a separate makefile for the rtlib so I am assuming that those options apply to the rtlib as well

[dodicat]

Thanks srvaldez.
O2 option the best option I suppose.