Firstly, credit to cbruce for bringing this PRNG to our attention and, secondly, credit to counting_pine for porting the minimal C implementation to FreeBASIC.
The wrapper overhead was expensive so I have used the age old method of code repetition by using the same pcg code in the head of each of the wrappers.
If we employ no seeding the first two random 32 bit values will be zero. All random number generators should be 'warmed up'. In An experimental exploration of Marsaglia’s xorshift generators,scrambled by Sebastiano Vigna he produces a graph showing the number of samples which are required to warm up some PRNGs, including his own. It is quite clear that the the larger the period of a generator the greater the number of samples should be executed before using a generator. Not surprisingly the Mersenne Twister needs a lot of warming up. Since PCG shares some similarities with Vigna's PRNGs and we have a period of 2^64 with this minimal PCG I have opted for a warm up of 200 samples. This should prove more than adequate and only takes one micro-second to complete, so you will not have time to put the kettle on. <smile>
The warm up is done in the function RandomizePCG32( Seed )
Seed is optional. If no seed is given a cryptographic random Ulongint is used from FB's PRNG option 5. To repeat a sequence a seed is required. Quite frankly, anything will do: 123456, say, is fine.
For single precision we simply execute PCG32S. For range, 0 to 100 for example, we simply execute PCG32R(0,100).
The code at the end of this post has the code to add to your apps and an example usage.
I now use these compiler option switches at all times: -gen gcc -Wc -O1
With no compiler optimization the throughput collapses.
This is a typical output using RandomizePCG32 - that is an unknown seed.
Code: Select all
0.3077029
0.5957927
0.65415
0.7683138
0.6042935
0.1181773
0.07335913
0.6330329
0.7821538
0.4110711
116
204
86
7
211
37
128
147
240
157
Averge of 10^8 singles: 0.5000328068581712
Average of 10^8 range (0,100): 49.99446299
Estimate of 'raw' throughput (Millions per second)
Singles: 219.113
Range: 180.463
In 64 bit mode I get 222.410 and 206.462 for singles and range respectively.
A period of 2^64 may be regarded as small but just to put that into perspective it would take over 3,000 years on my machine before a sequence repetition and that is just executing random numbers and doing nothing with them. Also, if we used 2^32 random numbers then we would be using 2^32/2^64 = 1/2^32 of those available, that is 23*10^-9% so the chance of 'landing' on numbers previously used is very small. Of course,the chance decreases as the number used decreases.
srvaldez has published the full implementation of PCG in [1] but I have not had a chance to look at it yet.
Some folk may find the following code adequate for their purposes.
Code: Select all
'' *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org
'' Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)
'' FreeBASIC translation by Matthew Fearnley 2017-06-04
' Include the code between ' ******/' ****** in your apps
' **********************************************
Type pcg32_random_t
As Ulongint state, inc
End Type
Dim Shared pcg32 As pcg32_random_t
Function pcg32_random_r(Byval rng As pcg32_random_t Ptr) As Ulong
Dim As Ulongint 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 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 RandomizePCG32( Seed As Ulongint = 0 )
Dim i As Integer
If seed = 0 Then
Randomize , 5
pcg32.state = Rnd*(2^64-1)
Else
pcg32.state = seed
End If
For i As Integer = 1 To 200
pcg32_random_r(@pcg32)
Next i
End Sub
Function PCG32S() As Single
Dim tempvar As Ulong
Dim rng As pcg32_random_t Ptr = @pcg32
Dim As Ulongint 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 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 eax, dword Ptr [TempVar]
movd xmm0, eax
psrlq xmm0, 9
mov eax, 1
cvtsi2ss xmm1, eax
por xmm0, xmm1
subss xmm0, xmm1
movd [Function], xmm0
End Asm
End Function
Function PCG32R( Byval One As Long, Byval Two As Long ) As Long
Dim tempvar As Ulong
Dim rng As pcg32_random_t Ptr = @pcg32
Dim As Ulongint 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 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 Now1LT2
xchg eax, ecx
Now1LT2:
Sub eax, ecx
inc eax
jz doTheRnd
mul edx
Add edx, ecx
doTheRnd:
mov [Function], edx
End Asm
End Function
' *****************************************************
RandomizePCG32
Dim As Ulong I, N
' Knock out 10 singles
For i = 1 To 10
Print PCG32S
Next
Print
' Knock out 10 range
For i = 1 To 10
Print PCG32R(0,255)
Next
Print
Dim As Double tot
' Knock out average of 10^8 singles
For i = 1 To 10^8
tot += PCG32S
Next
tot = tot/10^8
Print "Averge of 10^8 singles: ";tot
Print
tot = 0
' Knock out average of 10^8 range (0,100)
For i = 1 To 10^8
tot += PCG32R(0,100)
Next
tot = tot/10^8
Print "Average of 10^8 range (0,100): ";tot
Print
Dim As Double t
Print "Estimate of 'raw' throughput (Millions per second)"
Print
t = Timer
For i = 1 To 10^8
PCG32S
Next
t = Timer - t
Print "Singles: ";Using "###.###";10^8/(t*1000000)
t = Timer
For i = 1 To 10^8
PCG32R(0,100)
Next
t = Timer - t
Print "Range: ";Using "###.###";10^8/(t*1000000)
Sleep