FreeBASIC's PRNG #2

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

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

@dafhi

Apologies, I missed your post.

Have you compared 'i = (742938285 * i ) Mod (2^31-1)', from Cook's blog, with 'i = 742938285 * ( bit_rotate( i, 1 ) + 1 )'?

I didn't think that we could use a circular shift for such a modulus. [1]

If you did use 'i = (742938285 * i ) Mod (2^31-1)' then your code is fine.

Anyway, I will not test using PractRand because the algorithm is an LCG and will fail fairly quickly.

I used both methods above for computing state and the Chi-Squared for both were OK so the distributions are uniform. The other metrics were not bad either. So, maybe I am wrong with [1].

I hope my agent, paul doe, is not expecting a commision. <smile>
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Post by dodicat »

Could you give an example of a mean value for say 100000000 pcg32 values (.randse) in [0,1].
Your randomize ,5 seems to be chocking everything.

I am surprised that random ,4 fared so well in my (trivial) montecarlo/4 and mean things.
Also I am surprised how poorly randomize,5 fared, without the speed headache.
But my little tests are only in the freebasic world of course.

The rnd operations are in the rtl math_rnd.c
the first few lines
...............
#include "fb.h"

#if defined HOST_WIN32
#include <windows.h>
#include <wincrypt.h>
#elif defined HOST_LINUX
#include <fcntl.h>
#endif

#define RND_AUTO 0
#define RND_CRT 1
#define RND_FAST 2
#define RND_MTWIST 3
#define RND_QB 4
#define RND_REAL 5
..................

Also, if you didn't already know, #4 and #5 start seed is not confined to integer.

Code: Select all

 

for z as long=0 to 5
    print "randomize timer, ";z
for n as long=1 to 10
    randomize timer,z
    print rnd
next
print
next z

sleep   
I might rustle up a CHI squared test for #4.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Post by dodicat »

Here is the qb (#4) source used in fb.

Code: Select all


#include "crt.bi"
#define INITIAL_SEED	327680
dim shared as uint32_t iseed = INITIAL_SEED

function hRnd_QB cdecl ( n as single ) as double
	union ftoi
		f as single
		i as uint32_t
	end union

	dim _ftoi as ftoi
	
	if ( n <> 0.0 ) then
		if ( n < 0.0 ) then
			_ftoi.f = n
			dim as uint32_t s = _ftoi.i
			iseed = s + ( s shr 24 )
		end if
		iseed = ( ( iseed * &hFD43FD ) + &hC39EC3 ) and &hFFFFFF
	end if
	
	return cast(single, iseed) / cast(single, &h1000000)
end function


    
for n as long=1 to 20
    randomize timer
   print hRnd_QB(1)
next
sleep
 
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: FreeBASIC's PRNG #2

Post by jj2007 »

deltarho[1859] wrote:You know what, jj2007, in view of the fact that PCG32II knocks spots off anything that I have and there is a lot more to it than meets the eye perhaps I should write a Help file for it.
Indeed, that would be a good idea. In the meantime, I studied your code (excellent btw) and found that the real work is all done here:

Code: Select all

Private Function pcg32.randse() As Double
  Dim TempVar As Ulong
  Dim As Ulongint oldstate = this.state
  this.state = oldstate * 6364136223846793005ULL + (this.seq Or 1)
  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))
    Return TempVar/4294967296.0
End Function
Under the hood, in the disassembly, it looks messy, and I am tempted to rewrite it in assembly. It could be a lot shorter and faster. If I find the time...

What does the "u" stand for in oldstate Shr 18u?
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

@dodicat

10^8 mean value: 0.4999966466239317
Your randomize ,5 seems to be chocking everything.
What do you mean?
Also, if you didn't already know, #4 and #5 start seed is not confined to integer.
They will be converted to an integer. All the generators use integers internally.
I might rustle up a CHI squared test for #4.
Attached is a zip of chiSquare2g.exe written by John Gleason, in 2006, at the PowerBASIC forums. John has a sense of humour. If you execute by double clicking on it we get a messagebox saying "Oopsie, problem with file,try 'er again." It is designed to be run from a command line: chiSquare2g <filepath>. The output is a message box and is similar to ENT. I have just run it on a file and got a % of 96.319, ENT got 96.3%, ie strong evidence that the data is not random.

I will have a look at the qb (#4) source later. Thanks for that.

chiSquare2g.zip
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

@jj2007
I am tempted to rewrite it in assembly
That would be way beyond me. <smile>

"u" stands for UINTEGER. See Help file under "Standard Data Type Limits"
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

FB's #4 is from Microsoft Visual Basic (6 and earlier). It is an LCG so PractRand will tear it apart - I will try it though.

Added: I already have, according to my opening post. Blimey, I cannot even remember two days ago but I can remember in vivid detail things which happened years ago. You know what that is a sign of - yep, flaming old age. Came out of my local supermarket the other day and could I hell remember where I'd parked the car. A chap walked up to me and asked if I could spare some change. I described my car and said, "If you can find it before me there is a quid in it for you". I found it before he did but I gave him a quid anyway.
Last edited by deltarho[1859] on Sep 07, 2018 12:06, edited 1 time in total.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Post by dodicat »

Actually seems like #4,#5 are unaffected by the randomize seed.

Code: Select all


for z as long=0 to 3
  randomize 1,4
for n as long=1 to 10
    print rnd
next
print
next z

sleep  
deltarho said
They will be converted to an integer. All the generators use integers internally.

I know they use integers internally, but the conversions??

Sorry missed chocking.
I mean your rand is rendered slow by using #5.
If I wanted 1000000 random vectors say,.
I am only thinking fb practicalities, getting random numbers.
I normally use randomize ,2 for speed.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

dodicat wrote:Actually seems like #4,#5 are unaffected by the randomize seed.
#5 is cryptographic so it does not need a seed.

There is something strange going on with #4. Your last code should produce four identical blocks but it does not. Each execution of the code produces the same four blocks.

It is using a hard-wired seed. That is not intended. We have a bug, Houston.
I mean your rand is rendered slow by using #5.
I still do not understand. #5 is used in MyRandomize before calling Get64Bit. MyRandomize is completed in no time at all. #5 has nothing to with the random number generation.

#2 uses a LCG so will be fast. However, the above 'real life' tests show #2 and PCG32II are 'neck and neck' so why use #2?
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Post by dodicat »

That's why I posted.
#4, #5 are unaffected by a seed.
Note the help extract:

For any given seed, each algorithm will produce a specific, deterministic sequence of numbers for that seed. If you want each call to Randomize to produce a different sequence of numbers, a seed that is not quite predictable should be used - for example, the value returned from Timer. Omitting the seed parameter will use a value based on this.

And here was me thinking that randomize ,4 accepted a single as a seed.
Silly silly.
Everything must be tested in fb before use.
I have the opinion now than random generators are ten a penny.
(Random being the stickler word, as it has been since maths was invented)

I made one up.

Code: Select all


dim shared as long seed, k

function _rnd as double
k=k+1:if k>100019 then k=0
do
    seed=((seed*100+k) mod 100019)
    if seed<=99999 then return (seed/100000)
loop 
end function

sub _randomize(n as long=0)
    seed=n:k=n
    end sub

#define pseudo_range(f,l) _rnd*(l-f)+f

screen 19,32,,64
dim as ulong white=rgb(255,255,255)
_randomize 1

do
     pset(pseudo_range(50,750),pseudo_range(50,550)),_rnd*white
loop until len(inkey)
    
Could be a non randomness about (probability 1).
What the heck, it is more useful than #5 here.
There is a chance that #5 is crap on this machine, while others zip away at it.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Post by dodicat »

Example of using your prng

Code: Select all

 ' Generator PCG32II
' *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org
' Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)
' pcg32.rand is an adaptation of a FreeBASIC translation by Matthew Fearnley
' (aka FreeBASIC's counting_pine) 2017-06-04
' Object duplication method by FreeBASIC's St_W

Type pcg32
  Public:
  declare Constructor
  Declare Sub MyRandomize( seed As Ulongint = 0, seq As Ulongint = 0 )
  Declare Function rand() As Ulong
   Declare Function randse() As Double
  Declare Function range( Byval One As Long, Byval Two As Long ) As Long
   Declare Function GetSeed() As Ulongint
   Declare Function GetSeq() As Ulongint
  declare Function Gauss() as double
  state As Ulongint
  seq   As Ulongint
   seed  As Ulongint
End Type

Private Function Get64Bit() As UlongInt
return (Cast( Ulongint, Rnd*(2^32) ) Shl 32) Or Cast( Ulongint, Rnd*(2^32) )
End Function

Private Function pcg32.rand() As Ulong
  Dim As Ulongint oldstate = this.state
  this.state = oldstate * 6364136223846793005ULL + (this.seq Or 1)
  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

Private Sub pcg32.MyRandomize( seed As Ulongint = 0, seq As Ulongint = 0 )
  Dim i As Integer
  
  Randomize , 5
  If seed = 0 Then
    If seq = 0 Then
      this.state = Get64Bit
      this.seq = Get64Bit
    Else
      this.state = Get64Bit
      this.seq = seq
    End If
  Else
    If seq = 0 Then
      this.state = seed
      this.seq = Get64Bit
    Else
      this.state = seed
      this.seq = seq
    End If
  End If
  This.Seed = This.state ' Save initial state
  ' Warm up generator - essential
  For i As Integer = 1 To 200
    this.rand()
  Next i
End Sub

Private Function pcg32.randse() As Double
  Dim TempVar As Ulong
  Dim As Ulongint oldstate = this.state
  this.state = oldstate * 6364136223846793005ULL + (this.seq Or 1)
  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))
    Return TempVar/4294967296.0
End Function

Private Function pcg32.range( Byval One As Long, Byval Two As Long ) As Long
  Dim TempVar As Ulong
  Dim As Ulongint oldstate = this.state
  ' Advance internal state
  this.state = oldstate * 6364136223846793005ULL + (this.seq Or 1)
  ' rotate32((state ^ (state >> 18)) >> 27, state >> 59)
  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 0f
    xchg eax, ecx
    0:
    Sub eax, ecx
    inc eax
    jz 1f
    mul edx
    Add edx, ecx
    1:
    mov [Function], edx
  End Asm
End Function

Private Function pcg32.GetSeed() As Ulongint
   Return This.Seed
End Function

Private Function pcg32.GetSeq() As Ulongint
   Return This.Seq
End Function

constructor pcg32
  This.MyRandomize
end constructor

Private Function pcg32.Gauss As double
Static As Long u2_cached
Static As double u1, u2, x1, x2, w
 
  If u2_cached = -1 Then
    u2_cached = 0
    Function = u2
  Else
    Do
      x1 = This.randse
      x2 = This.randse
      w = x1 * x1 + x2 * x2
    Loop While w >= 1
    w = Sqr( -2 * Log(w)/w )
    u1 = x1 * w
    u2 = x2 * w
    u2_cached = -1
    Function = u1
  End If
 
End Function


type range
    as double _max,_min,sd
    as long _maxi,_mini
end type

function mean(a() as double,R as range) as double
    R=type(-1e10,1e10,0,0,0)
    dim as long lb=lbound(a),ub=ubound(a)
    dim as double acc,m
    for n as long=lb to ub
        acc+=a(n)
        if R._max<a(n) then R._max=a(n):R._maxi=n
        if R._min>a(n) then R._min=a(n):R._mini=n
    next
    m=acc/(ub-lb+1)
    acc=0
    for n as long=lb to ub
        acc+=(a(n)-m)*(a(n)-m)
    next
    R.sd =sqr(acc/(ub-lb))
    return m
end function
dim as double t
t=timer
redim as pcg32 b(1000000)
print "OK ";timer-t
redim as double a(ubound(b))
for n as long=0 to ubound(b)
    a(n)=b(n).randse
next

print "done"
print
dim as range r
var m=mean(a(),r)
print "mean ";m
print "smallest ";r._min
print "biggest ";r._max
print "Standard deviation ";r.sd
sleep

  
My result:

Code: Select all

OK  19.47298196237534
done

mean  0.4998160133220828
smallest  1.110835000872612e-006
biggest  0.9999997459817678
Standard deviation  0.2887414772820628
  
With #2 instead of #5 it takes 3 seconds.
(Is it my machine?)
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: FreeBASIC's PRNG #2

Post by srvaldez »

macOS High Sierra, FreeBASIC Compiler - Version 1.06.0 (03-01-2018), built for darwin-x86_64 (64bit)
fbc -w all -gen gcc -asm att -Wc -O2 "FreeBASIC's PRNG #2.bas"
my result of code as posted

Code: Select all

OK  2.277823209762573
done

mean  0.50002027459726
smallest  1.565553247928619e-06
biggest  0.9999944709707052
Standard deviation  0.288721683280281
changing line #38 from Randomize , 5 to Randomize , 2 the results are

Code: Select all

OK  0.3546240329742432
done

mean  0.5001081421162514
smallest  4.788860678672791e-06
biggest  0.9999831323511899
Standard deviation  0.2888171671143238
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Post by dodicat »

Thanks srvaldez.
Even with 64 bit -O3 it takes 18 seconds here.
It must be this computer.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: FreeBASIC's PRNG #2

Post by srvaldez »

hi dodicat
don't feel bad, the same PC in a Windows 10 VM times at about 24 seconds for the 32-bit version, tried various optimization options.
dafhi
Posts: 1640
Joined: Jun 04, 2005 9:51

Re: FreeBASIC's PRNG #2

Post by dafhi »

deltarho[1859] wrote:@dafhi

Apologies, I missed your post.

Have you compared 'i = (742938285 * i ) Mod (2^31-1)', from Cook's blog, with 'i = 742938285 * ( bit_rotate( i, 1 ) + 1 )'?
I have no testing experiencing other than staring at randomness pixels
deltarho[1859] wrote: If you did use 'i = (742938285 * i ) Mod (2^31-1)' then your code is fine.
crickets
deltarho[1859] wrote: Anyway, I will not test using PractRand because the algorithm is an LCG and will fail fairly quickly.
not an LCG :-)
deltarho[1859] wrote: I used both methods above for computing state and the Chi-Squared for both were OK so the distributions are uniform. The other metrics were not bad either. So, maybe I am wrong with [1].
*eats popcorn*
deltarho[1859] wrote: I hope my agent, paul doe, is not expecting a commision. <smile>
He's a good man
Post Reply