pcg32rr.bi - minimal PCG32_random_r implementation ( (c) Melissa E. O'Neill )

User projects written in or related to FreeBASIC.
MrSwiss
Posts: 3190
Joined: Jun 02, 2013 9:27
Location: Switzerland

pcg32rr.bi - minimal PCG32_random_r implementation ( (c) Melissa E. O'Neill )

Postby MrSwiss » Jan 12, 2019 20:58

Hi all,

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 -----
Two little tests (showing use in code):

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 -----
MrSwiss
Posts: 3190
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: pcg32rr.bi - minimal PCG32_random_r implementation ( (c) Melissa E. O'Neill )

Postby MrSwiss » Jan 12, 2019 21:19

Here is the code I've used for the above mentioned tests.

IMPORTANT:
- only for FBC 64 bit (otherwise, undefined results)
you'll have to FIRST uncomment the following two procedures in 'pcg32rr.bi':
- pcg32dbl2() and
- pcg32dbl3()
It's best to comment them after speed test(s) is/are run.

Code:

Code: Select all

/'
    pcg32rr.bi -- copyright (c) 2019, MrSwiss (except as noted in comments)
'/

' testing code - start - pcg32rr_speed-test.bas
#Include "pcg32rr.bi"

Using pcg

Dim As pcg32_random_t   pcgt1, pcgt2, pcgt3 ' use three instances
Dim As Double t, tt1, tt2, tt3, tt4

init_pcg_rnd(@pcgt1)                    ' initialize the type's
init_pcg_rnd(@pcgt2)
init_pcg_rnd(@pcgt3)

' ----- user code start
Print "Speed test running, please wait ..."
Print

For j As UInteger = 1 To 1e3            ' 1 to 1'000
    ' using multiplier const
    t = Timer
    For i As UInteger = 1 To 1e6        ' 1 to 1'000'000
        pcg32dbl(@pcgt1)
    Next
    tt1 += Timer - t
    ' using bit-shift & ptr cast
    t = Timer
    For i As UInteger = 1 To 1e6
        pcg32dbl2(@pcgt2)
    Next
    tt2 += Timer - t
    ' using bit-shift & union
    t = Timer
    For i As UInteger = 1 To 1e6
        pcg32dbl3(@pcgt3)
    Next
    tt3 += Timer - t
    ' Mersenne Twister (built in)
    t = Timer
    For i As UInteger = 1 To 1e6
        Rnd()
    Next
    tt4 += Timer - t
Next

Print "Times taken ... 1'000 x 1'000'000 run's each"
Print "pcg32dbl(@pcgt1)  "; tt1; Tab(40); "using multiplier const"
Print "pcg32dbl2(@pcgt2) "; tt2; Tab(40); "using shift's & ptr cast"
Print "pcg32dbl3(@pcgt3) "; tt3; Tab(40); "using shift's & union"
Print
Print "Mersenne Twister  "; tt4; Tab(40); "(FBC built in)"
Print
' ----- user code end

Sleep
' testing code - end    ' ----- EOF -----
a typical result wrote:Speed test running, please wait ...

Times taken ... 1'000 x 1'000'000 run's each
pcg32dbl(@pcgt1) 1.161500201094896 using multiplier const
pcg32dbl2(@pcgt2) 1.162808296037838 using shift's & ptr cast
pcg32dbl3(@pcgt3) 1.165347297443077 using shift's & union

Mersenne Twister 5.043453700141981 (FBC built in)

(GCC options: -O max [difference to -O 3 is insignificant])
MrSwiss
Posts: 3190
Joined: Jun 02, 2013 9:27
Location: Switzerland

Re: pcg32rr.bi - minimal PCG32_random_r implementation ( (c) Melissa E. O'Neill )

Postby MrSwiss » Mar 01, 2019 13:18

Another example, using four instances of pcg32_random_t, with the DrawRhomb() Sub
See: Tips & Tricks/DrawRhomb() Sub ..., using Mersenne Twister (built in).
To demonstrate the use of multiple (differently seeded) random generation.

Code: Select all

' DrawRhomb_Sub-test1.bas -- (c) 2019-02-26, MrSwiss
'
' compile: -s gui
'
#Include "pcg32rr.bi"   ' https://www.freebasic.net/forum/viewtopic.php?f=8&t=27343
Using pcg               ' include pcg namespace

Sub DrawRhomb( _                        ' draw a rhomb(us) (aka: diamond, lozenge)
    ByVal xc    As Single, _            ' X-axis rhomb center
    ByVal yc    As Single, _            ' y-axis rhomb center
    ByVal wid   As Single, _            ' width of rhomb
    ByVal hei   As Single, _            ' height of rhomb
    ByVal ang   As Single=0.0, _        ' rhomb angle of rotation (optional)
    ByVal bcl   As ULong=&hFFFFFFFF, _  ' border color (optional)
    ByVal fcl   As ULong=0, _           ' fill color (optional)
    ByVal fll   As Boolean=FALSE _      ' fill flag (default: none)
    )
    #Ifndef Pi                          ' only if not externally defined
    Static As Double    Pi  = 4d * Atn(1d)      ' semi circle (rad)
    Static As Double    d_r = Atn(1d) / 45.0    ' deg to rad mul. factor
    #EndIf
    Dim As Single   SinAng, CosAng, Theta, D_Theta, _   ' angles
                    x, y, rx, ry        ' position points
    ' initialize local variables
    ang *= d_r : SinAng = Sin(ang) : CosAng = Cos(ang)
    Theta = 0.0 : D_Theta = Pi * .5     ' 90° in RAD
    ' calculate the first point (2D vector)
    x = wid * Cos(Theta) : y = hei * Sin(Theta)
    ' rotate it
    rx = xc + x * CosAng + y * SinAng
    ry = yc - x * SinAng + y * CosAng
    ' set start point (aka: graphics cursor)
    PSet (rx, ry)                       ' doesn't change anything visible!
    ' draw the rhomb
    For i As UInteger = 0 To 3          ' four points
        Theta += D_Theta                ' add delta, 90° RAD
        x = wid * Cos(Theta)            ' new pos.
        y = hei * Sin(Theta)
        rx = xc + x * CosAng + y * SinAng   ' rotate it
        ry = yc - x * SinAng + y * CosAng
        Line -(rx, ry), bcl             ' draw line (relative, border color)
    Next
    If fll Then Paint (xc, yc), fcl, bcl' fill if requested (fill color)
End Sub


' demo code (uses: PCG32rr.bi & it's namespace, see: at top)
ScreenRes(1024, 800, 32)

Dim As pcg32_random_t   p1, p2, p3, p4  ' 4 instances of pcg32 (type)
' user defined initializers (instead of: randomly picked)
init_pcg_user(@p1, &h0, &h00FF00FF00FF00FF, TRUE)   ' TRUE = run it 8000 times
init_pcg_user(@p2, &h1000, &h0000FFFF0000FFFF, TRUE)
init_pcg_user(@p3, &h1000000, &h7F007F007F007F00, TRUE)
init_pcg_user(@p4, &h1000000000, &h007F007F007F007F, TRUE)

Dim As ULong    angle = 180, x = 511, y = 399, _    ' max. angle | center/center
                w = pcg32range(@p3, 180, 300), _    ' random width
                h = pcg32range(@p2, 100, 200), _    ' random height
                bc = &hFFFF3F00, fc = &hFF007FFF    ' preset colors
Dim As Boolean  clsf = FALSE, fllf = TRUE           ' Cls-/fill-flags

' MAIN-Loop
Do
    ScreenLock
    If clsf Then Cls : clsf = FALSE
    DrawRhomb(x, y, w, h, angle, bc, fc, fllf)
    ScreenUnLock
    angle -= 5                          ' smaller angle = CW rotation!
    If angle = 0 Then                   ' reset condition and new init
        angle = 180 : clsf = TRUE : fllf = Not fllf
        bc = pcg32range(@p1, &hFF000000, &hFFFFFFFF)    ' 24 bit RGB only
        fc = pcg32range(@p4, &hFF3F3F3F, &hFFBFBFBF)
        x  = pcg32range(@p2, 199, 823)
        y  = pcg32range(@p3, 199, 599)
        w  = pcg32range(@p2, 180, 320)
        h  = pcg32range(@p3, 180, 320)
    End If
    Sleep(16, 1)                        ' approx. 60 Hz (aka: FPS)
Loop Until InKey() = Chr(255, 107)      ' [Alt]+[F4] press or a [X] click
' ----- EOF -----
Tested: FBC 64 - WIN, 1.06.0

Return to “Projects”

Who is online

Users browsing this forum: No registered users and 5 guests