FreeBASIC's PRNG #2
-
- Posts: 4308
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
@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>
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>
Re: FreeBASIC's PRNG #2
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.
I might rustle up a CHI squared test for #4.
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
Re: FreeBASIC's PRNG #2
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
Re: FreeBASIC's PRNG #2
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: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.
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
What does the "u" stand for in oldstate Shr 18u?
-
- Posts: 4308
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
@dodicat
10^8 mean value: 0.4999966466239317
I will have a look at the qb (#4) source later. Thanks for that.
chiSquare2g.zip
10^8 mean value: 0.4999966466239317
What do you mean?Your randomize ,5 seems to be chocking everything.
They will be converted to an integer. All the generators use integers internally.Also, if you didn't already know, #4 and #5 start seed is not confined to integer.
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 might rustle up a CHI squared test for #4.
I will have a look at the qb (#4) source later. Thanks for that.
chiSquare2g.zip
-
- Posts: 4308
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
@jj2007
"u" stands for UINTEGER. See Help file under "Standard Data Type Limits"
That would be way beyond me. <smile>I am tempted to rewrite it in assembly
"u" stands for UINTEGER. See Help file under "Standard Data Type Limits"
-
- Posts: 4308
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
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.
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.
Re: FreeBASIC's PRNG #2
Actually seems like #4,#5 are unaffected by the randomize seed.
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.
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
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.
-
- Posts: 4308
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
#5 is cryptographic so it does not need a seed.dodicat wrote:Actually seems like #4,#5 are unaffected by the randomize 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 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.I mean your rand is rendered slow by using #5.
#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?
Re: FreeBASIC's PRNG #2
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.
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.
#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)
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.
Re: FreeBASIC's PRNG #2
Example of using your prng
My result:
With #2 instead of #5 it takes 3 seconds.
(Is it my machine?)
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
Code: Select all
OK 19.47298196237534
done
mean 0.4998160133220828
smallest 1.110835000872612e-006
biggest 0.9999997459817678
Standard deviation 0.2887414772820628
(Is it my machine?)
Re: FreeBASIC's PRNG #2
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
changing line #38 from Randomize , 5 to Randomize , 2 the results are
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
Code: Select all
OK 0.3546240329742432
done
mean 0.5001081421162514
smallest 4.788860678672791e-06
biggest 0.9999831323511899
Standard deviation 0.2888171671143238
Re: FreeBASIC's PRNG #2
Thanks srvaldez.
Even with 64 bit -O3 it takes 18 seconds here.
It must be this computer.
Even with 64 bit -O3 it takes 18 seconds here.
It must be this computer.
Re: FreeBASIC's PRNG #2
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.
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.
Re: FreeBASIC's PRNG #2
I have no testing experiencing other than staring at randomness pixelsdeltarho[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 )'?
cricketsdeltarho[1859] wrote: If you did use 'i = (742938285 * i ) Mod (2^31-1)' then your code is fine.
not an LCG :-)deltarho[1859] wrote: Anyway, I will not test using PractRand because the algorithm is an LCG and will fail fairly quickly.
*eats popcorn*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].
He's a good mandeltarho[1859] wrote: I hope my agent, paul doe, is not expecting a commision. <smile>