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] »

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

Re: FreeBASIC's PRNG #2

Post by dodicat »

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.

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

  
Just end by closing the console (saves loop time looking for a key)
srvaldez
Posts: 3374
Joined: Sep 25, 2005 21:54

Re: FreeBASIC's PRNG #2

Post by srvaldez »

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' "
hello deltarho[1859]
use

Code: Select all

Asm "nop" ' Get off my loop!
and compile with -asm att
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"
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

@dodicat

Code: Select all

function rnd2() as double
dim as ulongint rand=Get64Bit
rand = (6364136223846793005ull * rand + 1442695040888963407ull) And 18446744073709551615ull
return  rand/2^64
end function
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.
dafhi
Posts: 1641
Joined: Jun 04, 2005 9:51

Re: FreeBASIC's PRNG #2

Post by dafhi »

deltarho[1859] wrote:
dafhi wrote:you can drop the and :D
Your second post lacking detail to the extent that I do not know what you are talking about. Perhaps it is just me.

Code: Select all

comment mark -->' And 18446744073709551615ull
jj2007
Posts: 2326
Joined: Oct 23, 2016 15:28
Location: Roma, Italia
Contact:

Re: FreeBASIC's PRNG #2

Post by jj2007 »

The QWORD 18446744073709551615 is -1, therefore the AND doesn't do anything useful.
dodicat
Posts: 7979
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Post by dodicat »

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

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

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

Re: FreeBASIC's PRNG #2

Post by jj2007 »

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
If you get different results with / without the And 2^64-1, please show. I am curious.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

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

Re: FreeBASIC's PRNG #2

Post by dodicat »

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

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

       
   
And the code I used

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



  
NOTE not only is random 0,5 not very good, it is desperately slow.
dafhi
Posts: 1641
Joined: Jun 04, 2005 9:51

Re: FreeBASIC's PRNG #2

Post by dafhi »

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

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

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.

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

Re: FreeBASIC's PRNG #2

Post by jj2007 »

deltarho[1859] wrote:Here are some outputs from Fourmilab's ENT.exe employing a bit analysis of 100MB of random data.
But you wrote earlier that ENT.exe is not good enough, right?

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

Re: FreeBASIC's PRNG #2

Post by deltarho[1859] »

jj2007 wrote:But you wrote earlier that ENT.exe is not good enough, right?
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.
Is the PCG32 code posted here ("Here is what I am using now") still valid?
Yes.
It builds fine but I have not understood how to use it.
Maybe you did not follow me when I wrote the 'main' PCG thread - you did call there a few times.

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
We use Gauss as in 'Mean + pcg32A.Gauss * Standard_deviation'

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
PCG32II is thread safe.

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
There are two other functions, GetSeed and Getseq, which simply tell us what the current seed and sequence are.

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.
Post Reply