FreeBASIC's PRNG #2
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
BTW, the 'Asm nop' used above gives a 'compiling C FAILED' with gcc 8.1 64-bit: " internal compiler error: output_operand: invalid use of register 'frame' "
Re: FreeBASIC's PRNG #2
Going for 1/4* montecarlo (pi/4)
Then randomize 4 seems the best for accuracy.
It reaches a steady 5 decimal places after about 20 seconds.
The others never seem to get 5 places even after a coffee run.
Just end by closing the console (saves loop time looking for a key)
Then randomize 4 seems the best for accuracy.
It reaches a steady 5 decimal places after about 20 seconds.
The others never seem to get 5 places even after a coffee run.
Code: Select all
'deltarho[]
Function Get64Bit() As Ulongint
Return (Cast( Ulongint, Rnd*(2^32) ) Shl 32) Or Cast( Ulongint, Rnd*(2^32) )
End Function
function rnd2() as double
dim as ulongint rand=Get64Bit
rand = (6364136223846793005ull * rand + 1442695040888963407ull) And 18446744073709551615ull
return rand/2^64
end function
#define inbox(x,y,r) x>0 and x <r and y>0 and y<r
#define incircle(x,y,r) x*x + y*y < r*r
randomize ,4
dim as long r=1
dim as double x,y
dim as double countbox,countcircle,counter
dim as ulongint u=2^16-1'2^24-1
do
counter+=1
x=rnd
y=rnd
if inbox(x,y,r) then countbox+=1
if incircle(x,y,r ) then countcircle+=1
if (counter and (u)) =0 then 'print regularily
counter=0
locate 3,,0
print countcircle/countbox
print atn(1)
end if
loop
sleep
Re: FreeBASIC's PRNG #2
hello deltarho[1859]deltarho[1859] wrote:BTW, the 'Asm nop' used above gives a 'compiling C FAILED' with gcc 8.1 64-bit: " internal compiler error: output_operand: invalid use of register 'frame' "
use
Code: Select all
Asm "nop" ' Get off my loop!
without the asm "nop" I get 1.#INF MHz, with asm "nop" I get 1455 MHz
changing dim As Ulong i to dim shared As Ulong i and using asm "nop" I get about 470 MHz
FreeBASIC Compiler - Version 1.06.0 (09-03-2018), built for win64 (64bit)
WinFBE_Suite\FreeBASIC-1.06.0\fbc64.exe -w all -asm att -O 3 -v "%f"
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
@dodicat
is tantamount to executing Randomize each and every time Rnd is executed. rand should be a shared variable and Get64Bit executed only once. Alternatively, rand could be static within rnd2 and Get64bit used once at rnd2's head. That would slow down rnd2 but I'm not sure, without checking, by how much.
@srvaldez
Thanks for that. Nice to see someone else using the 'latest and greatest' builds. The more of us who did that the more likely a stable 1.06 will get released.
Code: Select all
function rnd2() as double
dim as ulongint rand=Get64Bit
rand = (6364136223846793005ull * rand + 1442695040888963407ull) And 18446744073709551615ull
return rand/2^64
end function
@srvaldez
Thanks for that. Nice to see someone else using the 'latest and greatest' builds. The more of us who did that the more likely a stable 1.06 will get released.
Re: FreeBASIC's PRNG #2
deltarho[1859] wrote:Your second post lacking detail to the extent that I do not know what you are talking about. Perhaps it is just me.dafhi wrote:you can drop the and :D
Code: Select all
comment mark -->' And 18446744073709551615ull
Re: FreeBASIC's PRNG #2
The QWORD 18446744073709551615 is -1, therefore the AND doesn't do anything useful.
Re: FreeBASIC's PRNG #2
deltarho
Sorry about not using your rand correctly.
But my main theme was
randomize ,4
Advertised as
4 - Uses a function that is designed to give the same random number sequences as QBASIC. This should be stable across all platforms, and provides 24-bit precision, with a low degree of randomness.
But it seems to produce the best randomness (5 places in pi/4 within a few seconds).
FreeBASIC 1.05.0 (The last official build)
Only a handful of users will use 1.06.0 I would guess.
Maybe with it now included in the Winfbe suite another handful or two.
Sorry about not using your rand correctly.
But my main theme was
randomize ,4
Advertised as
4 - Uses a function that is designed to give the same random number sequences as QBASIC. This should be stable across all platforms, and provides 24-bit precision, with a low degree of randomness.
But it seems to produce the best randomness (5 places in pi/4 within a few seconds).
FreeBASIC 1.05.0 (The last official build)
Only a handful of users will use 1.06.0 I would guess.
Maybe with it now included in the Winfbe suite another handful or two.
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
@dafhi and jj2007
For 'Mod 2^n' we can use 'And 2^n-1' when n is a multiple of 2.
So, for 'Mod 2^64' we can use 'And 2^64-1'
18446744073709551615ull = 2^64-1
@dodicat
Randomness does not come into it. What matters is the uniformity of the scattering. Theoretically, any algorithm which produces a uniform distribution will do. Whether the distribution is built randomly is neither here nor there. 24-bit is Single precision. The precision does not come into it either- uniformity is the key. The total number of points scattered is also important. The larger the total the larger will be both countcircle and countbox in your code resulting in more significant figures in the approximation.
For 'Mod 2^n' we can use 'And 2^n-1' when n is a multiple of 2.
So, for 'Mod 2^64' we can use 'And 2^64-1'
18446744073709551615ull = 2^64-1
@dodicat
Randomness does not come into it. What matters is the uniformity of the scattering. Theoretically, any algorithm which produces a uniform distribution will do. Whether the distribution is built randomly is neither here nor there. 24-bit is Single precision. The precision does not come into it either- uniformity is the key. The total number of points scattered is also important. The larger the total the larger will be both countcircle and countbox in your code resulting in more significant figures in the approximation.
Re: FreeBASIC's PRNG #2
If you get different results with / without the And 2^64-1, please show. I am curious.deltarho[1859] wrote:For 'Mod 2^n' we can use 'And 2^n-1' when n is a multiple of 2.
So, for 'Mod 2^64' we can use 'And 2^64-1'
18446744073709551615ull = 2^64-1
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
In theory, I am correct.
However, (2^64-1)/2^64 = 0.9999999999999999999458 ie 19 significant figures.
UlongInt has 19+ significant figures and Double has 15+ significant figures.
So, with our highly inaccurate PCs we can 'drop the and'.
So, dafhi and jj2007, I stand corrected.
I thought that I may get a reasonable performance boost with less work to do but I cannot measure any improvement.
Anyway, thanks gents.
However, (2^64-1)/2^64 = 0.9999999999999999999458 ie 19 significant figures.
UlongInt has 19+ significant figures and Double has 15+ significant figures.
So, with our highly inaccurate PCs we can 'drop the and'.
So, dafhi and jj2007, I stand corrected.
I thought that I may get a reasonable performance boost with less work to do but I cannot measure any improvement.
Anyway, thanks gents.
Re: FreeBASIC's PRNG #2
100000000 instances of a random number.
expectation -- a uniform distribution
looking at offset from the mean (expected to be .9999999997671694/2 in most cases)
As in my previous code randomize , 4 seems the best.
used randomize through 1 to 5 then c runtime (crudish mean) then deltarho's rnd
results from 64 bit fb 1.050
And the code I used
NOTE not only is random 0,5 not very good, it is desperately slow.
expectation -- a uniform distribution
looking at offset from the mean (expected to be .9999999997671694/2 in most cases)
As in my previous code randomize , 4 seems the best.
used randomize through 1 to 5 then c runtime (crudish mean) then deltarho's rnd
results from 64 bit fb 1.050
Code: Select all
randomize 0,0
offset from expected (mean) 3.977254664055963e-006
mean 0.5000039771382487 expected mean 0.4999999998835847
smallest 9.313225746154785e-010
biggest 0.9999999965075404
Standard deviation 0.2886589088592121
randomize 0,1
offset from expected (mean) 5.493224234931482e-005
mean 0.4999450676412354 expected mean 0.4999999998835847
smallest 0
biggest 0.999969482421875
Standard deviation 0.288673736185414
randomize 0,2
offset from expected (mean) 2.284705608379678e-005
mean 0.4999771528275009 expected mean 0.4999999998835847
smallest 4.656612873077393e-010
biggest 0.9999999802093953
Standard deviation 0.2886803978142094
randomize 0,3
offset from expected (mean) 3.977254664055963e-006
mean 0.5000039771382487 expected mean 0.4999999998835847
smallest 9.313225746154785e-010
biggest 0.9999999965075404
Standard deviation 0.2886589088592121
randomize 0,4
offset from expected (mean) 3.453724622459742e-006
mean 0.4999965461589622 expected mean 0.4999999998835847
smallest 0
biggest 0.9999999403953552
Standard deviation 0.2886751456277132
randomize 0,5
offset from expected (mean) 2.871379171459632e-005
mean 0.4999712860918701 expected mean 0.4999999998835847
smallest 4.656612873077393e-010
biggest 0.9999999990686774
Standard deviation 0.2886579723166759
C runtime
offset from expected (mean) 5.306510963098976e-005
mean 0.500053065109631 expected mean 0.5
smallest 0
biggest 1
Standard deviation 0.2886787037708659
deltarho[]
offset from expected (mean) 2.220751048742642e-005
mean 0.5000222073940721 expected mean 0.4999999998835847 (not sure of deltarho's expected value)
smallest 7.194511839633908e-009
biggest 0.999999990834915
Standard deviation 0.288669628358582
done
Code: Select all
#include "crt.bi"
const c = RAND_MAX
function rnd2() as double
return rand/c
end function
'deltarho[]
Function Get64Bit() As Ulongint
Return (Cast( Ulongint, Rnd*(2^32) ) Shl 32) Or Cast( Ulongint, Rnd*(2^32) )
End Function
dim shared as ulongint _rand:_rand=Get64Bit
function rnd3() as double
_rand = (6364136223846793005ull * _rand + 1442695040888963407ull) And 18446744073709551615ull
return _rand/2^64
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 expected
for z as long=0 to 7
if z=6 then expected=.5 else expected=0.9999999997671694/2
if z<=5 then randomize 0 ,z else randomize 0
dim as ulongint lim=100000000
redim as double a(lim)
for n as long=0 to lim
if z<=5 then a(n)=rnd
if z=6 then a(n)=rnd2
if z=7 then a(n)=rnd3
next
dim as range r
var m=mean(a(),r)
if z<= 5 then print ,"randomize 0," +str(z)
if z=6 then print "C runtime"
if z=7 then print "deltarho[]"
print
print "offset from expected (mean) ";abs(m-expected)
print "mean ";m,"expected mean ";expected
print "smallest ";r._min
print "biggest ";r._max
print "Standard deviation ";r.sd
print
print
next z
print "done"
sleep
Re: FreeBASIC's PRNG #2
paul doe suggested i run my algorithm by you, see how it fares
Code: Select all
function bit_rotate(x as ulong, r as ulong) as ulong
return (x shr r) or (x shl (-r and 31))
end function
sub RCSG_32f(byref i as ulong)
/'
https://www.johndcook.com/blog/2017/07/05/simple-random-number-generator/
# Linear congruence generator (LCG) constants
z = 20170705 # seed
a = 742938285 # multiplier
'/
i = 742938285 * ( bit_rotate( i, 1 ) + 1 )
End sub
dim as ulong state = 20170705
for i as integer = 1 to 1e1
RCSG_32f state
? state
Next
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
Here are some outputs from Fourmilab's ENT.exe employing a bit analysis of 100MB of random data.
Chi Square is quite good at determining a distribution's uniformity. A percentage value in the range 10 < x% < 90 is deemed to provide no evidence that data is not uniform. We cannot say, for example, that 60% is better than 70% - we are either within the envelope or we are not.
I have not been able to establish the algorithm used by FB #4 but, from what I have read, the QBasic algorithm is pretty 'ropey'. This is borne out by the following figures.
With FB's #4 Chi Square coming in at 0.01% then there is very strong evidence that the data is not uniform. The QBasic algorithm, from what I have read, does poorly with Marsaglia's Die Hard Suite. The figures for Mersenne Twister and Knuth64 are OK.
Chi Square is quite good at determining a distribution's uniformity. A percentage value in the range 10 < x% < 90 is deemed to provide no evidence that data is not uniform. We cannot say, for example, that 60% is better than 70% - we are either within the envelope or we are not.
I have not been able to establish the algorithm used by FB #4 but, from what I have read, the QBasic algorithm is pretty 'ropey'. This is borne out by the following figures.
With FB's #4 Chi Square coming in at 0.01% then there is very strong evidence that the data is not uniform. The QBasic algorithm, from what I have read, does poorly with Marsaglia's Die Hard Suite. The figures for Mersenne Twister and Knuth64 are OK.
Code: Select all
F:\ENT>ent -b 100MB.txt
FB #3 Mersenne Twister
Chi square distribution for 838860800 samples is 0.04, and randomly
would exceed this value 83.37 percent of the times.
Arithmetic mean value of data bits is 0.5000 (0.5 = random).
Monte Carlo value for Pi is 3.142015806 (error 0.01 percent).
Serial correlation coefficient is 0.000012 (totally uncorrelated = 0.0).
FB #4 QBasic
Chi square distribution for 838860800 samples is 52430357.01, and randomly
would exceed this value less than 0.01 percent of the times.
Arithmetic mean value of data bits is 0.3750 (0.5 = random).
Monte Carlo value for Pi is 3.570999205 (error 13.67 percent).
Serial correlation coefficient is 0.166660 (totally uncorrelated = 0.0).
Knuth64
Chi square distribution for 838860800 samples is 0.72, and randomly
would exceed this value 39.45 percent of the times.
Arithmetic mean value of data bits is 0.5000 (0.5 = random).
Monte Carlo value for Pi is 3.141754194 (error 0.01 percent).
Serial correlation coefficient is 0.000013 (totally uncorrelated = 0.0).
Re: FreeBASIC's PRNG #2
But you wrote earlier that ENT.exe is not good enough, right?deltarho[1859] wrote:Here are some outputs from Fourmilab's ENT.exe employing a bit analysis of 100MB of random data.
Is the PCG32 code posted here ("Here is what I am using now") still valid? It builds fine but I have not understood how to use it. A simple for loop showing the first 100 numbers would be helpful to understand...
-
- Posts: 4292
- Joined: Jan 02, 2017 0:34
- Location: UK
- Contact:
Re: FreeBASIC's PRNG #2
Correct. I used ENT to specifically examine Chi Square because that is a good indicator of a distribution's uniformity. From a randomness perspective, ENT is not very good compared with PractRand.jj2007 wrote:But you wrote earlier that ENT.exe is not good enough, right?
Yes.Is the PCG32 code posted here ("Here is what I am using now") still valid?
Maybe you did not follow me when I wrote the 'main' PCG thread - you did call there a few times.It builds fine but I have not understood how to use it.
Anyway, here is an example usage illustrating Ulongs, [0,1), range(a,b) and Gauss.
Code: Select all
#Include "PCG32II.bas"
Dim Shared pcg32A As pcg32
Dim As Ulong i
' Knock Out ULongs
Print "Ulongs" : Print
For i = 1 To 6
Print pcg32A.rand
Next
Print
' Knock Out [0,1) 32-Bit granularity
Print "[0,1)" : Print
For i = 1 To 6
Print pcg32A.randse
Next
Print
' Knock Out a range
Print "range(a,b) eg (0,255)" : Print
For i = 1 To 6
Print pcg32A.range(0,255)
Next
Print
' Knock Out normal distribution
Print "Normal distribution" : Print
For i = 1 To 6
Print pcg32A.Gauss
Next
Sleep
Anyone not familiar with the normal distribution will need to do a spot of reading to understand that.
There is a randomize function similar to FreeBASIC: MyRandomize( Seed, Sequence )
Most generators provide a single sequence. Take Mersenne Twister (MT) for example and imagine the sequence as a circle. For a given seed we enter the circle at the same position and then move around the circle for each random number request. Once the sequence is exhausted by using all of the generator's period then we find ourselves back at the starting position and move around the circle again. With MT that will take a while because it has a very large period. A generator with a period of 2^32 will see us back at the starting position after 4 Gigs of random numbers have been requested. With PCG32II we have available 2^63 streams. For a given seed with two generators with differing streams, we will have two different sequences ie two different circles. For a given sequence with two generators with differing seeds they will be using the same circle. This is similar to using FreeBASIC's Randomize more than once and should be avoided - we could end up with a territorial overlap. If we want to repeat a sequence we simply chose the same seed and the same sequence. Both the seed and sequence need to be a UlongInt. MyRandomize with a random seed and random sequence is invoked in a Constructor so we do not need to use it if we want a random 'starting' position. The generator is also 'warmed up' for us.
Here is an example using streams:
Code: Select all
#Include "PCG32II.bas"
Dim Shared pcg32A As pcg32
Dim as Ulong i
pcg32A.MyRandomize()
Print "Random seed and seq"
For i = 1 to 5
Print pcg32A.randse
Next
Print
pcg32A.MyRandomize(,99)
Print "Random seed, seq=99"
For i = 1 to 5
Print pcg32A.randse
Next
Print
pcg32A.MyRandomize(,99)
For i = 1 to 5
Print pcg32A.randse
Next
Print
pcg32A.MyRandomize(5,)
Print "Seed=5, random seq"
For i = 1 to 5
Print pcg32A.randse
Next
Print
pcg32A.MyRandomize(5,)
For i = 1 to 5
Print pcg32A.randse
Next
Print
pcg32A.MyRandomize(5,99)
Print "Seed=5 and seq=99"
For i = 1 to 5
Print pcg32A.randse
Next
Print
pcg32A.MyRandomize(5,99)
For i = 1 to 5
Print pcg32A.randse
Next
Print
Sleep
Here is an example using a single thread ( Twin = False ) or two threads ( Twin = True ). This shows, with 'Twin = True', that there are no collisions between the thread's generators.
Code: Select all
#Include "PCG32II.bas"
#define Twin false
Dim Shared pcg32A As pcg32
Dim Shared pcg32B As pcg32
pcg32A.MyRandomize
pcg32B.MyRandomize
Dim As Ulong i
Dim As Double t1, y
Dim Shared As Double t2
Dim shared As Ulong counter
Dim As Any Ptr hThread
Sub SecondThread( byval dummy as any ptr)
Dim As Ulong i
Dim As Double y
t2 = Timer
For i = 1 To counter
y = pcg32A.randse()
Next
t2 = Timer - t2
End Sub
counter = 10^8 ' 100 million
#If Twin
hThread = Threadcreate( @SecondThread, 0 )
#endif
t1 = Timer
For i = 1 To counter
y = pcg32B.randse()
Next
t1 = Timer - t1
#if Twin
Threadwait( hThread )
#endif
#if Twin
Print Int(100/t1);" miil/sec", Int(100/t2);" mill/sec"
#else
Print Int(100/t1);" mill/sec"
#endif
Sleep
All the above are covered in the 'main' PCG32II thread but it is a long thread.
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. Anyone visiting the 'main' PCG thread for the first time will find it heavy going; not heavy technically but it really is a long read and new ideas keep popping up throughout.
I will write a Help file for it. <smile> I have just got an update to HelpNDoc so this is an opportunity to give it a whirl.