The 32 bit minimal C code from Melissa E. O'Neill translated to FreeBASIC.
Only the really necessary procedures added, to deal with issues like:
- Initialization (random initial values, Mersenne Twister based, FB's CRT lib.)
- Initialization (user supplied initial values), user controlled init
- returning a [0, 1) interval, as Double (as the 'built in' randomizers do)
- returning a positive range (can easily be converted to other ranges)
Main aim is a simple way, to use it ... preferably, with FBC 64, due to GCC's
optimization options, without which, even Mersenne Twister is faster!
(Above by comparing 'Apple with Apples', the inverval [0, 1) procedure! And
NOT the raw ULong output!)
By contrast, to try and implement all the maybe needed, possible variants.
IMPORTANT:
GCC optimisation options -O 2, should be used as minimum, for speed!
(I've NOT tested the -O 1 option! Nor any FBC 32 bit tests made!)
Optimal speed results are with -O 3 (and above, like -O max).
Above -O 3 there are only insignificant changes seen ...
Additional adjustments like -fpu sse, did not seem to have a impact.
For modularisation purposes, the namespace 'pcg' is used.
Include guards (FB style) in place.
A header only implementation (all the code is in the 'pcg32rr.bi' file).
I'm not interested, in making a lot of comments about it, since others have
done so, who are in a far better position, to do so.
If you want to know more about it, use the provided links below.
Random Numbers Generation Basics (M.E. O'Neill)
PCG32 Minimal C Implementation (M.E. O'Neill) the actually translated C-code
Testing the PCG random number generator (John D. Cook)
Information about *true-random* vs. *pseudo-random* generated values
PCG Passes PractRand (M.E. O'Neill)
PCG seems to have passed, just about all the existing tests, such as:
- TestU01 (by original author), initial testing
- DIEHARDER (extended from DIEHARD, Marsaglia et.al.)
- NIST (not certain, their site is currently down, due to 'lockdown')
- PractRand (by original author), up to 32 TB
The 'pcg32rr.bi':
Code: Select all
/'
pcg32rr.bi -- Copyright (c) 2019, MrSwiss (except, as noted in comments)
--------------------------------------------------------------------------
Licensed under Apache License 2.0 (NO WARRANTY, etc. see: pcg-random.org)
The same, as the original C-code by Melissa E. O'Neill (for simplicity)
--------------------------------------------------------------------------
Instead of translating every used variable, I've used type alias(es).
The only other translation is therefore, the translation of operators.
Apart from that, the addition of: 'Dim As' which isn't needed in 'C'.
Removed all semi-colon's (C's line end marker).
--------------------------------------------------------------------------
'/
#Pragma Once ' FB style include guard
Namespace pcg ' for modularity
Type uint16_t As UShort ' substitute the FB specific data-types
Type uint32_t As ULong ' make type alias(es)
Type uint64_t As ULongInt
' internally used constants
Const As ULong ULmax = &hFFFFFFFFul ' used by Mersenne Twister (initializers)
Const As Double DmulF = 1.0 / (ULmax + 1ull) ' ULong-range multiplier factor _
' used by pcg32dbl() Function
'// *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org
'// Licensed under Apache License 2.0 (NO WARRANTY, etc. see ebsite)
'// translated C-code - start
Type pcg32_random_t
As uint64_t state
As uint64_t inc
End Type
Function pcg32rr( _ ' pcg32_random_r renamed (shorter)
ByVal rng As pcg32_random_t Ptr _
) As uint32_t ' aka: ULong
Dim As uint64_t oldstate = rng->state
' Advance internal state
rng->state = oldstate * 6364136223846793005ULL + (rng->inc Or 1)
' Calculate output function (XSH RR), uses old state for max ILP
Dim As uint32_t xorshifted = ((oldstate Shr 18) Xor oldstate) Shr 27
Dim As uint32_t rot = oldstate Shr 59
Return (xorshifted Shr rot) Or (xorshifted Shl ((-rot) And 31))
End Function
'// translated C-code - end
' just for demo purposes, because unlikely to be used, otherwise ...
' probably better explained as: pcg32x2rr (calling the same sequence
' twice or at your otion, using two different sequences), 64 bit result
Function pcg64rr( _ ' 64 bit, using: 2 x 32 bit
ByVal pt1 As pcg32_random_t Ptr, _ ' mandatory
ByVal pt2 As pcg32_random_t Ptr = 0 _ ' optional
) As uint64_t ' aka: ULongInt
#Ifdef __FB_64BIT__ ' FBC 64 (no type-cast's required)
If pt1 = pt2 OrElse pt2 = 0 Then ' depending on: provided parameter's
Return (pcg32rr(pt1) Shl 32) + pcg32rr(pt1) ' same sequence (twice)
Else
Return (pcg32rr(pt1) Shl 32) + pcg32rr(pt2) ' diff. sequence's
End If
#Else ' FBC 32 (with required type-cast's)
If pt1 = pt2 OrElse pt2 = 0 Then ' as above
Return (CULngInt(pcg32rr(pt1)) Shl 32) + pcg32rr(pt1)
Else
Return (CULngInt(pcg32rr(pt1)) Shl 32) + pcg32rr(pt2)
End If
#EndIf
End Function
' fastest method of returning a [0, 1) interval as double (32bit granularity)
' (result of performed speed test vs. the below 'commented out' implementations,
' with various GCC optimisation options applied, averaged = 5 executed tests
' a single testrun is 1 Billion call's, aka: 1'000 * 1'000'000, see testcode)
Function pcg32dbl( _
ByVal pt1 As pcg32_random_t Ptr _
) As Double
Return pcg32rr(pt1) * DmulF
End Function
'' for spped tests only, as well as FBC 64 only!!! (just delete leading ')
'Function pcg32dbl2( _
' ByVal pt1 As pcg32_random_t Ptr _
' ) As Double
' Dim As uint64_t u64 = (&h3FF Shl 52) + (pcg32rr(pt1) Shl 20)
' Return (*cptr(Double Ptr, @u64) - 1)
'End Function
'
'Function pcg32dbl3( _
' ByVal pt1 As pcg32_random_t Ptr _
' ) As Double
' Union ULL2DBL
' As uint64_t u
' As Double d
' End Union
' Dim As ULL2DBL conv
' conv.u = (&h3FF Shl 52) + (pcg32rr(pt1) Shl 20)
' Return (conv.d - 1)
'End Function
'' surprisingly, with most GCC optimisation options: pcg32dbl() is fastest!
'' equally, with unoptimized GCC, Mersenne Twister is fastest! (FBC's CRT)
Function pcg32range( _
ByVal pt1 As pcg32_random_t Ptr, _
ByVal l_lim As uint32_t, _
ByVal h_lim As uint32_t _
) As uint32_t
Return Int(pcg32dbl(pt1) * ((h_lim + 1) - l_lim) + l_lim) ' round to floor
End Function
' used internally by init_pcg_rnd() and/or init_pcg_user()
Private Sub init_MT( _ ' initializer for Mersenne Twister
ByVal size As UShort = 8000 _ ' 1000 x 8 call's
)
Randomize(Timer, 3)
For i As UInteger = 1 To (size Shr 3) ' size \ 8 (fast)
Rnd() : Rnd() : Rnd() : Rnd() : _
Rnd() : Rnd() : Rnd() : Rnd()
Next
End Sub
Private Function MT64() As uint64_t ' 2 x 32 bits combined (MT)
#Ifdef __FB_64BIT__ ' FBC 64 (no type-cast's required)
Return ((Rnd * ULmax) Shl 32) + (Rnd * ULmax)
#Else ' FBC 32 (with required type-cast's)
Return (CULngInt(Rnd * ULmax) Shl 32) + CULng(Rnd * ULmax)
#EndIf
End Function
Private Sub init_pcg32_t( _
ByVal pcg_t As pcg32_random_t Ptr, _
ByVal size As UShort = 8000 _
)
For x As UInteger = 0 To (size Shr 3) ' size \ 8 (fast)
pcg32rr(pcg_t) : pcg32rr(pcg_t) : _
pcg32rr(pcg_t) : pcg32rr(pcg_t) : _
pcg32rr(pcg_t) : pcg32rr(pcg_t) : _
pcg32rr(pcg_t) : pcg32rr(pcg_t)
Next
End Sub
' end: used internally by init_pcg_rnd() and/or init_pcg_user()
Sub init_pcg_rnd( _ ' initializer (uses random values)
ByVal pcg_t As pcg32_random_t Ptr, _
ByVal init As Boolean = TRUE, _ ' manual switch-off option (initializer)
ByVal isize As uint16_t = 8000 _
)
Static As Boolean _once_
' init_MT() is run once only
If _once_ = FALSE Then
init_MT(isize) : _once_ = TRUE
End If
' get random uint64_t's (ULongInt)
pcg_t->state = MT64()
pcg_t->inc = MT64()
If init Then init_pcg32_t(pcg_t, isize) ' run initializer (default: yes)
End Sub
Sub init_pcg_user( _ ' initializer (uses user defined values)
ByVal pcg_t As pcg32_random_t Ptr, _
ByVal state As uint64_t, _
ByVal inc As uint64_t, _
ByVal init As Boolean = FALSE, _ ' manual switch-on option (initializer)
ByVal isize As uint16_t = 8000 _
)
pcg_t->state = state
pcg_t->inc = inc
If init Then init_pcg32_t(pcg_t, isize) ' run initializer (default: no)
End Sub
End Namespace ' pcg -- end
' ----- EOF -----Code: Select all
/'
pcg32rr.bi -- copyright (c) 2019, MrSwiss (except as noted in comments)
--------------------------------------------------------------------------
For additional information and/or details read: pcg32rr.bi
'/
' testing code - start - pcg32rr_test1.bas (with 'using' directive)
#Include "pcg32rr.bi"
Using pcg ' include pcg-namespace into 'global-namespace'
Dim As pcg32_random_t pcgt1, pcgt2 ' use two instances (for 64/2 x 32 bit only needed!)
init_pcg_rnd(@pcgt1) ' initialize the type (random values, default init = TRUE)
init_pcg_user(@pcgt2, 123456789, 987654321, TRUE, 16000) ' user suppled init values (incl. init)
' ----- user code start
For i As UInteger = 0 To 59
Select Case As Const i Mod 4
Case 0
Print Tab(11);
Print pcg32rr(@pcgt1) ' @ = type ptr = MANDATORY!
Case 1
Print Tab(11);
Print Str(pcg32dbl(@pcgt2)) ' returns: ULong converted to Double [0, 1)
Case 2
Print Tab(11); ' return range: a signed byte (modified from ubyte)
Print pcg32range(@pcgt1, 0, 255) - 128 ' any positive range (incl. 0) _
' however, easy to convert to +/- or -/0 etc.
Case 3
Print Tab(11);
Print pcg64rr(@pcgt1, @pcgt2) ' using both sequences (2 x 32 bit)
Print
End Select
Next
' ----- user code end
Sleep
' testing code - end ' ----- EOF -----Code: Select all
/'
pcg32rr.bi -- copyright (c) 2019, MrSwiss (except as noted in comments)
--------------------------------------------------------------------------
For additional information and/or details read: pcg32rr.bi
'/
' testing code - start - pcg32rr_test2.bas (without 'using' directive)
#Include "pcg32rr.bi"
Dim As pcg.pcg32_random_t pcgt1, pcgt2 ' use two instances (for 64/2 x 32 bit only needed!)
pcg.init_pcg_rnd(@pcgt1) ' initialize the type (random values, default init = TRUE)
pcg.init_pcg_user(@pcgt2, 123456789, 987654321, TRUE, 16000) ' user suppled init values (incl. init)
' ----- user code start
For i As UInteger = 0 To 59
Select Case As Const i Mod 4
Case 0
Print Tab(11);
Print pcg.pcg32rr(@pcgt1) ' @ = type ptr = MANDATORY!
Case 1
Print Tab(11);
Print Str(pcg.pcg32dbl(@pcgt2)) ' returns: ULong converted to Double [0, 1)
Case 2
Print Tab(11); ' return range: a signed byte (modified from ubyte)
Print pcg.pcg32range(@pcgt1, 0, 255) - 128 ' any positive range (incl. 0) _
' however, easy to convert to +/- or -/0 etc.
Case 3
Print Tab(11);
Print pcg.pcg64rr(@pcgt1, @pcgt2) ' using both sequences (2 x 32 bit)
Print
End Select
Next
' ----- user code end
Sleep
' testing code - end ' ----- EOF -----