Single and Range wrappers for minimal PCG PRNG

General FreeBASIC programming questions.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

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 miss-reading 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.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

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
If you put two 64-bit 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 64-bit samples with multiplicative prediction resistance to build a 128-bit value, you end up with a random number with 128 bits of prediction resistance (2^128 * 2^128 = 2^256). Combine two of those 128-bit values together, and you get a 256-bit 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 Non-deterministic 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 800-90. RdRand complies with NIST SP 800-90A and RdSeed complies with 800-90B & 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 pre-Vista. 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 CBC-MAC 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^64-1 which could , of course, be a valid return albeit a very unlikely one. So, we could test for 2^64-1 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^64-1
 
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.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

Here is a Get64Bit() test.

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
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by jj2007 »

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=i-1
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.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

jj2007 wrote:What is the expected average here? I get strange values:
Expected value = (2^64-1)/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>
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by jj2007 »

deltarho[1859] wrote:Expected value = (2^64-1)/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^64-1)/2
print "sub: ";subtract

For i = 1 To 10000
  temp=Get64Bit()
  if i<=20 Then
	print temp
  endif
  Sum+=temp-subtract
Next
i=i-1
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
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

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>
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by jj2007 »

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.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

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.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Single and Range wrappers for minimal PCG PRNG

Post by dodicat »

You would expect the average RND to be less than .5.
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  
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

dodicat wrote:RND can be zero, but always falls short of one.
Do you know why it "falls short of one"?
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Single and Range wrappers for minimal PCG PRNG

Post by dodicat »

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 
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

Internally we are dealing with [0, 2^32-1] ie 2^32 values.

We normalize by dividing by 2^32 and get [0, (2^32-1)/2^32]

which gives [0, 1) since (2^32-1)/2^32 < 1

What we actually have is [0, 0.9999999997671694]
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Single and Range wrappers for minimal PCG PRNG

Post by dodicat »

Thanks Deltarho[].
I hadn't really put much thought into RND before.
Nice neat explanation.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Single and Range wrappers for minimal PCG PRNG

Post by deltarho[1859] »

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 32-bit mode. For 64-bit mode the increase is 25%; which is significant but the 32-bit mode increase is something else.
Post Reply