Single and Range wrappers for minimal PCG PRNG

 Posts: 2408
 Joined: Jan 02, 2017 0:34
 Location: UK
Re: Single and Range wrappers for minimal PCG PRNG
When I read the pcg paper it was not clear to me whether the sequence number, referred to as increment by O'Neill, had an upper limit of < 2^63 or < 2^64 so I played it safe by using < 2^63; effectively halving the number of sequences from 2^64 to 2^63.
O'Neill's blog, last entry 20 August 2017, fleshes out some aspects of pcg and it appears that the upper limit is, in fact, < 2^64. There is a restriction in that the sequence number must be odd. This is ensured by 'this.seq Or 1' in the code. If this.seq is odd the 'this.seq Or 1' is redundant. If this.seq is even then 'this.seq Or 1' makes it odd. Disallowing even numbered, halves, obviously, the total number of sequences that we can use. So, my statement about a key length of 127 bits per round is still true.
In the Sub pcg32.MyRandomize there are two instances of Get64Bit\2. They can now be changed to Get64Bit. We now have 0 <= seq <= 2^64  1. In use seq will be >= 1 by virtue of 'this.seq Or 1'.
Strictly speaking when I refer to different sequences  I should be referring to different streams.
All that my missreading did was to put an earlier block on how many steams we could use.
By the way the blog is a good read  that is where I found out that O'Neill is heavily into using PractRand now.
O'Neill's blog, last entry 20 August 2017, fleshes out some aspects of pcg and it appears that the upper limit is, in fact, < 2^64. There is a restriction in that the sequence number must be odd. This is ensured by 'this.seq Or 1' in the code. If this.seq is odd the 'this.seq Or 1' is redundant. If this.seq is even then 'this.seq Or 1' makes it odd. Disallowing even numbered, halves, obviously, the total number of sequences that we can use. So, my statement about a key length of 127 bits per round is still true.
In the Sub pcg32.MyRandomize there are two instances of Get64Bit\2. They can now be changed to Get64Bit. We now have 0 <= seq <= 2^64  1. In use seq will be >= 1 by virtue of 'this.seq Or 1'.
Strictly speaking when I refer to different sequences  I should be referring to different streams.
All that my missreading did was to put an earlier block on how many steams we could use.
By the way the blog is a good read  that is where I found out that O'Neill is heavily into using PractRand now.

 Posts: 2408
 Joined: Jan 02, 2017 0:34
 Location: UK
Re: Single and Range wrappers for minimal PCG PRNG
I have decided to replace Get64Bit() for Windows users. A Linux version is totally beyond me and I should not think it a walk in the park for experienced Linux fans.
As most of you may know Intel introduced RdRand to it's instruction set with Ivy Bridge processors and later. RdSeed was introduced with Broadwell processors and later.
There intended use is simply:
RdSeed for seeding another pseudorandom number generator (PRNG).
RdRand for all other purposes.
From The Difference Between RDRAND and RDSEED
'together' was misspelt.
We should have seen 2^64 * 2^64 = 2^128.
Clearly, the author John M. (Intel) had a bad hair day. <smile>
RdRand has additive prediction resistance and RdSeed has multiplicative prediction resistance. RdRand is described as a Cryptographically secure pseudorandom number generator (CSPRNG) and RdSeed as a Nondeterministic random bit generator.
From Digital Random Number Generator (DRNG) Software Implementation Guide there is at 4.2.6 code to convert RdRand to have fully forward and backward prediction resistance for those who do not have RdSeed. I have RdRand but not RdSeed.
Microsoft's BCryptGenRandom complies with NIST SP 80090. RdRand complies with NIST SP 80090A and RdSeed complies with 80090B & C. I am fairly confident then that BCryptGenRandom is only additive prediction resistant. It does not make sense to me for BCryptGenRandom to be multiplicative prediction resistant. The current GetBit64() uses FB's algorithm 5 which is in Windows Crypto API. I cannot remember what that complied with being preVista. We don't need 'Randomize , 5' now in MyRandomize.
It follows that we could use the Intel code to make BCryptGenRandom multiplicative prediction resistant. Using BCryptGenRandom means not just Windows but Windows Vista and later.
The code uses 512 * 128 bit random numbers ( 512 * AES blocks ), encrypts them with the CBCMAC mode of AES and uses the last 128 bits. Of course, the name of our game is not authentication so we use a random IV (Initial Values) and a random AES128 key; as does the Intel code.. There is no explanation why 512 was chosen. I have used 256 and use the last 64 bits ie half a block. 256 * 128 / 8 = 4KB which is where I think BCryptGenRandom refreshes the entropy pool.
Looking at the code it may seem that it will take a while to generate our 64 bits but, in fact, on my machine it takes about 180 microseconds so no chance to put the kettle on.
The error trapping was for development  the second parameter of BCryptGenerateSymmetricKey was hKey and it should have been @hKey. There should be no run time errors and it is not easy to report one. I am pushing out a messagebox but we will not be able to recover from that. In this case the return value is set at 2^641 which could , of course, be a valid return albeit a very unlikely one. So, we could test for 2^641 and recover from that even though we could be rejecting a valid return.
PCG32II now has a very posh seeding algorithm for MyRandomize for Windows Vista and later users.
New GetBit64()
As most of you may know Intel introduced RdRand to it's instruction set with Ivy Bridge processors and later. RdSeed was introduced with Broadwell processors and later.
There intended use is simply:
RdSeed for seeding another pseudorandom number generator (PRNG).
RdRand for all other purposes.
From The Difference Between RDRAND and RDSEED
If you put two 64bit values with additive prediction resistance togehter, the prediction resistance of the resulting value is only 65 bits (2^64 + 2^64 = 2^65).
'together' was misspelt.
If you use two 64bit samples with multiplicative prediction resistance to build a 128bit value, you end up with a random number with 128 bits of prediction resistance (2^128 * 2^128 = 2^256). Combine two of those 128bit values together, and you get a 256bit number with 256 bits of prediction resistance.
We should have seen 2^64 * 2^64 = 2^128.
Clearly, the author John M. (Intel) had a bad hair day. <smile>
RdRand has additive prediction resistance and RdSeed has multiplicative prediction resistance. RdRand is described as a Cryptographically secure pseudorandom number generator (CSPRNG) and RdSeed as a Nondeterministic random bit generator.
From Digital Random Number Generator (DRNG) Software Implementation Guide there is at 4.2.6 code to convert RdRand to have fully forward and backward prediction resistance for those who do not have RdSeed. I have RdRand but not RdSeed.
Microsoft's BCryptGenRandom complies with NIST SP 80090. RdRand complies with NIST SP 80090A and RdSeed complies with 80090B & C. I am fairly confident then that BCryptGenRandom is only additive prediction resistant. It does not make sense to me for BCryptGenRandom to be multiplicative prediction resistant. The current GetBit64() uses FB's algorithm 5 which is in Windows Crypto API. I cannot remember what that complied with being preVista. We don't need 'Randomize , 5' now in MyRandomize.
It follows that we could use the Intel code to make BCryptGenRandom multiplicative prediction resistant. Using BCryptGenRandom means not just Windows but Windows Vista and later.
The code uses 512 * 128 bit random numbers ( 512 * AES blocks ), encrypts them with the CBCMAC mode of AES and uses the last 128 bits. Of course, the name of our game is not authentication so we use a random IV (Initial Values) and a random AES128 key; as does the Intel code.. There is no explanation why 512 was chosen. I have used 256 and use the last 64 bits ie half a block. 256 * 128 / 8 = 4KB which is where I think BCryptGenRandom refreshes the entropy pool.
Looking at the code it may seem that it will take a while to generate our 64 bits but, in fact, on my machine it takes about 180 microseconds so no chance to put the kettle on.
The error trapping was for development  the second parameter of BCryptGenerateSymmetricKey was hKey and it should have been @hKey. There should be no run time errors and it is not easy to report one. I am pushing out a messagebox but we will not be able to recover from that. In this case the return value is set at 2^641 which could , of course, be a valid return albeit a very unlikely one. So, we could test for 2^641 and recover from that even though we could be rejecting a valid return.
PCG32II now has a very posh seeding algorithm for MyRandomize for Windows Vista and later users.
New GetBit64()
Code: Select all
#Include Once "windows.bi"
#Include Once "win\bcrypt.bi"
#Inclib "bcrypt"
#Include Once "win/wincrypt.bi"
#Inclib "crypt32"
Function Get64Bit() As Ulongint
Dim dwStatus As Ulong
Dim As BCRYPT_ALG_HANDLE Ptr hRand, hAESAlg
Dim As BCRYPT_KEY_HANDLE Ptr hKey
Dim As String IV, AESKey, sRandomBytes, sReturn, sTemp
Dim As Long lError, lSize
Dim As dword dwResult
dwStatus = BCryptOpenAlgorithmProvider(@hRand, BCRYPT_RNG_ALGORITHM, 0, 0) ' Prepare For Random number generation
If dwStatus <> 0 Then lError = 10 : Goto ErrorTrap
IV = String(16,0) ' Generate a Random IV
' Generate a random IV
dwStatus = BCryptGenRandom(hRand, Strptr( IV ), 16, 0)
If dwStatus <> 0 Then lError = 20 : Goto ErrorTrap
AESKey = String( 32, 0 ) ' Generate a Random encryption key
' Generate a random AES key
dwStatus = BCryptGenRandom(hRand, Strptr( AESKey ), 32, 0) ' 256 bits
If dwStatus <> 0 Then lError = 30 : Goto ErrorTrap
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
dwStatus = BCryptGenRandom(hRand, Strptr( sRandomBytes ), 4096, 0)
If dwStatus <> 0 Then lError = 40 : Goto ErrorTrap
dwStatus = BCryptCloseAlgorithmProvider(hRand, 0)
If dwStatus <> 0 Then lError = 50 : Goto ErrorTrap
dwStatus = BCryptOpenAlgorithmProvider( @hAESAlg, BCRYPT_AES_ALGORITHM, 0, 0 )
If dwStatus <> 0 Then lError = 60 : Goto ErrorTrap
' Use default CBC algorithm. 6th parameter of 16 indicates AES128
dwStatus = BCryptGenerateSymmetricKey( hAESAlg, @hKey, 0, 0, Strptr( AESKey ), 16, 0 ) ' We want hKey
If dwStatus <> 0 Then lError = 70 : Goto ErrorTrap
' 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
dwStatus = BCryptEnCrypt( hKey, Strptr( sRandomBytes ), Len( sRandomBytes ), 0, Strptr( IV ), _
16, Strptr( sRandomBytes ), Len( sRandomBytes ), @dwResult, 0 )
If dwStatus <> 0 Then lError = 80 : Goto ErrorTrap
BCryptCloseAlgorithmProvider( hAESAlg, 0 )
BCryptDestroyKey( hKey )
sReturn += Right( sRandomBytes, 8 ) ' Half of Last block of ciphertext  64 bits
Function = Cast(Ulongint, Cvlongint(sReturn))
Goto TidyUp
ErrorTrap:
MessageBox 0, "Error at Position " + Str(lError) + " " + Hex(dwStatus), "GetBit64() function", MB_ICONERROR
Function = 2^641
TidyUp:
If hRand <> 0 Then BCryptCloseAlgorithmProvider( hRand, 0 )
If hAESAlg <> 0 Then BCryptCloseAlgorithmProvider( hAESAlg, 0 )
If hKey <> 0 Then BCryptDestroyKey( hKey )
End Function
Last edited by deltarho[1859] on Feb 23, 2018 13:29, edited 1 time in total.

 Posts: 2408
 Joined: Jan 02, 2017 0:34
 Location: UK
Re: Single and Range wrappers for minimal PCG PRNG
Here is a Get64Bit() test.
Typical output:
Code: Select all
Dim As long i
Dim As Ulongint temp
For i = 1 To 20
temp = Get64Bit()
Print temp, Hex(temp, 16)
Next
Sleep
Typical output:
Code: Select all
1179217869368033519 105D6C67887980EF
8061447756788657677 6FE00405FEC1EE0D
11532571073257096558 A00BEC49F86F456E
5244478776237625113 48C821A919B80F19
1643232616617765285 16CDEF5BCCD669A5
1124810037282964205 0F9C20C4F80B82ED
1264657024004987329 118CF6DCE68519C1
15163745403503081578 D270721A8043786A
17587994337889448690 F4151BAA457B72F2
5955891729638608324 52A793F8F1CBCDC4
8450419287641279037 7545EB91E2F7CA3D
4912272513105289905 442BE5D3078E62B1
17856521655769314901 F7CF1BD688333655
15128184744814889080 D1F21BDF732CC878
16795082945703182868 E9141EF477DC0214
9378041748964242630 82257F489FCBF4C6
7892165538176175302 6D869ABE13920CC6
6908801614862817025 5FE0FE76DD8B1701
14347557251325062244 C71CC34DAFCF6C64
11575454837111675515 A0A446D84F9E5A7B
Re: Single and Range wrappers for minimal PCG PRNG
deltarho[1859] wrote:Here is a Get64Bit() test.
What is the expected average here? I get strange values:
Code: Select all
Dim As long i
Dim As Ulongint temp
Dim As Double Sum=0
For i = 1 To 100000
temp=Get64Bit()
' print temp
Sum+=temp
Next
i=i1
print " count=",i
Sum=Sum/i
print "average=",Sum
Sleep
Output:
Code: Select all
count= 100000
average= 9.22706565734578e+018
I have tried
Code: Select all
For i = 1 To 1000
temp=Get64Bit()2^63
if i<=20 Then
print temp
endif
Sum+=temp
Next
but no success  positive numbers only. With a similar algo I get e.g.
Code: Select all
8334395501155713472.
1329069443539533820.
2521605505986592764.
2513648907272060932.
6348034446003798020.
3307815087883943940.
4066175165531160572.
3922227352330829828.
428641309553590268.0
8910786327096590340.
100000000 iterations, average: [b]7.702250e+14[/b] generated in 980 ms
Last edited by jj2007 on Feb 23, 2018 17:35, edited 1 time in total.

 Posts: 2408
 Joined: Jan 02, 2017 0:34
 Location: UK
Re: Single and Range wrappers for minimal PCG PRNG
jj2007 wrote:What is the expected average here? I get strange values:
Expected value = (2^641)/2 ie Max value of Ulongint/2 ===> 9.223372036854776e+018
From docs:
ULONGINT 64 unsigned integer 0 +18446744073709551615 ull 19+
The code uses: Function = Cast(Ulongint, CVLongInt(sReturn)) where sReturn is the last 8 bytes of ciphertext as a String.
I have just run your code with 10,000,000 samples and got
Code: Select all
count= 10000000
average= 9.224077230537447e+018
It is getting there. I don't want to push the boat out too much  at my age I may never see the result. <laugh>
Re: Single and Range wrappers for minimal PCG PRNG
deltarho[1859] wrote:Expected value = (2^641)/2 ie Max value of Ulongint/2 ===> 9.223372036854776e+018
Code: Select all
Dim As long i
Dim As Ulongint temp
Dim As Double Sum=0, subtract=(2^641)/2
print "sub: ";subtract
For i = 1 To 10000
temp=Get64Bit()
if i<=20 Then
print temp
endif
Sum+=tempsubtract
Next
i=i1
print " count=",i
Sum=Sum/i
print "average=",Sum
Sleep
Should be close to zero now, right?
Code: Select all
sub: 9.223372036854776e+018
17915931214730098112
16981617639080706205
6851483647008877885
17667523011585827458
...
11465691595791434169
count= 10000
average= 4.997765434240366e+016
With 100,000 iterations:
Code: Select all
average= 3.48358037525256e+016
One Million (warning, takes ages):
Code: Select all
average= 218437685626371.6

 Posts: 2408
 Joined: Jan 02, 2017 0:34
 Location: UK
Re: Single and Range wrappers for minimal PCG PRNG
average= 218437685626371.6
Looks large but it isn't compared to Ulongint interval.
To get a feel of the convergence we should normalise that figure 218437685626371.6/2^64 = 0.00000118415.
The Sun is 93 million miles away but we are almost touching it compared with the 'width' of the universe. <smile>
Re: Single and Range wrappers for minimal PCG PRNG
You are right! In fact, I had calculated the inverse of the average (corrected above): I get 7.7e+14, which is a bit higher than yours (in absolute terms) but the same order of magnitude.

 Posts: 2408
 Joined: Jan 02, 2017 0:34
 Location: UK
Re: Single and Range wrappers for minimal PCG PRNG
the same order of magnitude
... or, similar in relative terms as opposed absolute terms.
Another way of looking at it.
I have just calculated the average of 100 million RND and got 0.5000039778499213. If we subtract 0.5 and multiply by 10^9 we get 3977.8499213. Should we be closer to zero? No, because we should also multiply 0.5 by 10^9. In relative terms they are the same.
Re: Single and Range wrappers for minimal PCG PRNG
You would expect the average RND to be less than .5.
RND can be zero, but always falls short of one.
RND can be zero, but always falls short of one.
Code: Select all
randomize 10
dim as double a,n,z,x
do
a=0
randomize n
for n as long=1 to 100000000
a+=rnd
next
x= a/100000000
n+=1
print x,n;" of 10"
z+=x
loop until n=10
print
print
print "Average average"
print z/10
sleep

 Posts: 2408
 Joined: Jan 02, 2017 0:34
 Location: UK
Re: Single and Range wrappers for minimal PCG PRNG
dodicat wrote:RND can be zero, but always falls short of one.
Do you know why it "falls short of one"?
Re: Single and Range wrappers for minimal PCG PRNG
I can only conclude
Code: Select all
#include "crt.bi"
printf(!"%22.18f = First step below one \n",nextafter(1,0))
printf(!"%22.18f = First step above zero \n",nextafter(0,1))
sleep

 Posts: 2408
 Joined: Jan 02, 2017 0:34
 Location: UK
Re: Single and Range wrappers for minimal PCG PRNG
Internally we are dealing with [0, 2^321] ie 2^32 values.
We normalize by dividing by 2^32 and get [0, (2^321)/2^32]
which gives [0, 1) since (2^321)/2^32 < 1
What we actually have is [0, 0.9999999997671694]
We normalize by dividing by 2^32 and get [0, (2^321)/2^32]
which gives [0, 1) since (2^321)/2^32 < 1
What we actually have is [0, 0.9999999997671694]
Re: Single and Range wrappers for minimal PCG PRNG
Thanks Deltarho[].
I hadn't really put much thought into RND before.
Nice neat explanation.
I hadn't really put much thought into RND before.
Nice neat explanation.

 Posts: 2408
 Joined: Jan 02, 2017 0:34
 Location: UK
Re: Single and Range wrappers for minimal PCG PRNG
The asm code in the integral range function has been replaced by a single FreeBASIC statement devised by dodicat. The throughput has increased by 60%  an astonishing increase.
Rather than post the source code the following link is a zip including the updated PCG32II.bas and PCG32II.chm. The Help hasn't changed other than to give dodicat the credit for the new integral range code.
Thanks, dodicat.
PCG32II.zip
PS The 60% refers to 32bit mode. For 64bit mode the increase is 25%; which is significant but the 32bit mode increase is something else.
Rather than post the source code the following link is a zip including the updated PCG32II.bas and PCG32II.chm. The Help hasn't changed other than to give dodicat the credit for the new integral range code.
Thanks, dodicat.
PCG32II.zip
PS The 60% refers to 32bit mode. For 64bit mode the increase is 25%; which is significant but the 32bit mode increase is something else.
Who is online
Users browsing this forum: No registered users and 5 guests