FreeBASIC's PRNG #2

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

Re: FreeBASIC's PRNG #2

Postby deltarho[1859] » Sep 07, 2018 18:40

I cannot believe what I am looking at.

redim as pcg32 b(1000000) will create one million generators!!!

Each generator is asked for one random number and then an average is determined.

It is no bloody wonder the code takes a long time. Each new generator created goes past MyRandomize which uses FB's #5.

Above I show how to use PCG32II.

I will definitely write a Help file.

If the intention is to determine the average of 1000000 random numbers in [0,1) got from ONE generator then this is how it is done.

Code: Select all

#Include Once "PCG32II.bas"
 
Dim As Ulong i
Dim As Double tot
Dim As pcg32 pcg
 
For i = 1 To 1000000
  tot += pcg.randse
Next
 
Print tot/1000000
 
Sleep
dodicat
Posts: 5913
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Postby dodicat » Sep 07, 2018 20:23

Thanks deltarho.
Much better.
Sorry I missed your method of use, but there are many posts and I still don't see it.
Anyway

Code: Select all

number of elements  100000001

mean  0.4999848280108567
smallest  1.30385160446167e-008
biggest  0.9999999979045242
Standard deviation  0.2886927923066043
Chi distance 16669211.47615323
time taken  for everything 9.080913471989334
   

That is the raw chi squared distance for .randse
srvaldez
Posts: 2063
Joined: Sep 25, 2005 21:54

Re: FreeBASIC's PRNG #2

Postby srvaldez » Sep 07, 2018 20:47

just in case someone else is wondering where "PCG32II.bas" can be found https://freebasic.net/forum/viewtopic.p ... 28#p233633

@deltarho[1859], do you have a version of PCG32II.bas that does not use asm?
FB on macOS does not support intel syntax, it only supports att syntax.
deltarho[1859]
Posts: 1921
Joined: Jan 02, 2017 0:34
Location: UK

Re: FreeBASIC's PRNG #2

Postby deltarho[1859] » Sep 07, 2018 21:05

@dodicat

Did you download the zip I linked to on page 3 post #5?

@srvaldez

Thanks. There is a lot of reading to do following the link. I have started a Help file. I would rather spend time on a graphics programming course or go to the dentist but it seems that I am rather good at it according to feedback from earlier chm publications.
dodicat
Posts: 5913
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Postby dodicat » Sep 07, 2018 21:38

message
Oopsie,problem with file, try 'er again -- in a little powerbasic messagebox.
I did try 'er again, no joy.

I just
∑(0-e)^2/e
O = observed
e=mean
To get chi squared.
srvaldez
Posts: 2063
Joined: Sep 25, 2005 21:54

Re: FreeBASIC's PRNG #2

Postby srvaldez » Sep 07, 2018 21:44

@deltarho[1859], could you test the following for correctness?
in the test that you posted at the link above, the function that has the asm code is not invoked and I don't know what result to expect.
I translated the inline asm to gcc inline asm

compile with -asm att
PCG32II.bas ''att asm version

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 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
  Private:
  state As Ulongint
  seq   As Ulongint
   seed  As Ulongint
End Type

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

Function pcg32.rand() 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
  Return (xorshifted Shr rot) Or (xorshifted Shl ((-rot) And 31))
End Function

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\2
    Else
      this.state = Get64Bit
      this.seq = seq
    End If
  Else
    If seq = 0 Then
      this.state = seed
      this.seq = Get64Bit\2
    Else
      this.state = seed
      this.seq = seq
    End If
  End If
  This.Seed = This.state ' Save initial state
  ' Warm up generator - essential
  ' We probably don't need 200 samples but it only takes one micro-second on my machine
  For i As Integer = 1 To 200
    this.rand()
  Next i
End Sub

Function pcg32.randse() As Double
  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))
    Return TempVar/4294967296.0
End Function

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 by John Gleason @ PowerBASIC forums

   Asm
      "movl %[TEMPVAR$1],%%edx \n" _
      "movl %[ONE$1],%%ecx \n" _
      "movl %[TWO$1],%%eax \n" _
      "cmp %%eax,%%ecx \n" _
      "jl 0f \n" _
      "xchg %%ecx,%%eax \n" _
      "0: \n" _
      "Subl %%ecx,%%eax \n" _
      "inc %%eax \n" _
      "jz 1f \n" _
      "mul %%edx \n" _
      "Addl %%ecx,%%edx \n" _
      "1: \n" _
      "movl %%edx,%[TEMPVAR$1] \n" _
      ::[TEMPVAR]"m"(TEMPVAR), [ONE]"m"(ONE), [TWO]"m"(TWO)  _
      :"rax", "rdx", "rcx"
   End Asm
   function = Tempvar
End Function

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

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


PCG32II-test.bas

Code: Select all

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

Re: FreeBASIC's PRNG #2

Postby deltarho[1859] » Sep 07, 2018 22:09

@dodicat

Oh, dear. In my post, I wrote: "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>"

ENT and chiSquare2g use ∑(0-e)^2/e but then use the percentage points of the chi-squared distribution to give a % likelihood of getting what we got - otherwise the chi-squared value is meaningless.

Having said that I have found out why FB #4 is getting lousy test results. If we execute Rnd*32 only 24 bits are filled so if we dump Ulongs to a file every fourth byte is empty. No wonder anything which examines 32 bits is going to balk. We need to dump only the 24 bits.
deltarho[1859]
Posts: 1921
Joined: Jan 02, 2017 0:34
Location: UK

Re: FreeBASIC's PRNG #2

Postby deltarho[1859] » Sep 07, 2018 22:23

@srvaldez

I ran PCG32II-test.bas including your version of PCG32II.bas and got

Code: Select all

False: 310 mill/sec
True:  268 mill/sec 270 mill/sec
srvaldez
Posts: 2063
Joined: Sep 25, 2005 21:54

Re: FreeBASIC's PRNG #2

Postby srvaldez » Sep 07, 2018 23:18

thanks deltarho[1859], it completely escaped my mind that the asm code I posted is 64-bit specific, that is, the register mangling declaration declares rax, rdx and rcx as clobbered registers, should anyone wish to compile to 32-bit they would need to change those registers to their 32-bit counterpart.
deltarho[1859]
Posts: 1921
Joined: Jan 02, 2017 0:34
Location: UK

Re: FreeBASIC's PRNG #2

Postby deltarho[1859] » Sep 07, 2018 23:53

srvaldez wrote:is 64-bit specific

I didn't spot that. It compiled in 32-bit. The results above are with 32-bit.

Anyway, in 64-bit mode.

Code: Select all

False: 972 mill/sec
True:  973 mill/sec 972 mill/sec

That proves without a shadow of a doubt that we have thread safety - we would not go through 10^8 random numbers without quite a few collisions otherwise; resulting in the 'bandwidth' of both generators being severely impaired.

I tested the '-asm att' and the range is not working.

Code: Select all

#include "PCG32IIsravaldez.bas"
Dim pcg As pcg32
Dim As Ulong i
For i = 1 To 10
  Print pcg.range(0,255)
Next

should churn out 10 numbers from (0,255) but a signed integer is being outputted.

BTW, the PCG32II.bas that you are using is not the same as the one dodicat is using. You need to add 'Declare Constructor' in 'Type pcg32' and add

Code: Select all

Constructor pcg32
  This.MyRandomize
End Constructor

otherwise MyRandomize will not get automatically executed.

The thread safety code worked because I use MyRandomize - it was written before the Constructor was added.

Incidentally, why are you testing with '-asm att'?

Added: I just spotted your addition above.
@deltarho[1859], do you have a version of PCG32II.bas that does not use asm?

No.
srvaldez
Posts: 2063
Joined: Sep 25, 2005 21:54

Re: FreeBASIC's PRNG #2

Postby srvaldez » Sep 08, 2018 0:31

deltarho[1859] wrote:I tested the '-asm att' and the range is not working.

Code: Select all

#include "PCG32IIsravaldez.bas"
Dim pcg As pcg32
Dim As Ulong i
For i = 1 To 10
  Print pcg.range(0,255)
Next

should churn out 10 numbers from (0,255) but a signed integer is being outputted.

apparently gcc is smart enough to cast the 64-bit registers to their 32-bit counterparts
it works with optimization levels 0, s or fast but not with O levels 1 thru 3
srvaldez
Posts: 2063
Joined: Sep 25, 2005 21:54

Re: FreeBASIC's PRNG #2

Postby srvaldez » Sep 08, 2018 0:40

deltarho[1859] wrote:Incidentally, why are you testing with '-asm att'?

so I can test on macOS
dafhi
Posts: 1245
Joined: Jun 04, 2005 9:51

Re: FreeBASIC's PRNG #2

Postby dafhi » Sep 08, 2018 2:37

well i tested the period of a low modulus datatype as a noob way of extrapolating performance of the real thing. it gets me thinking: could it be possible to tweak just one datatype to get a higher period? actually i think the answer is yes, and it'd have to be the state.

my original idea was, allowing an early state repeat early might not be a bad thing in the long run. might not be possible but if my curiosity holds i might pursue
deltarho[1859]
Posts: 1921
Joined: Jan 02, 2017 0:34
Location: UK

Re: FreeBASIC's PRNG #2

Postby deltarho[1859] » Sep 08, 2018 6:29

@dodicat

On page 3, post #3 you write
Here is the qb (#4) source used in fb.

That code is fundamentally flawed.

Firstly, 'randomize timer' should not be in a loop. Secondly, and more to the point, it is seeding the default generator and is then a waste of time with '-lang fb'.

Worst of all the seed is hard-wired forcing the sequence outputted to be the same for each invocation; which I spotted in an earlier post.

The only reason for including crt.bi is so that we can use uint32_t.

Who wrote that?
dodicat
Posts: 5913
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: FreeBASIC's PRNG #2

Postby dodicat » Sep 08, 2018 8:12

Hi deltarho
The randomize in a loop is not required , I should have removed it.
I got a bit excited finding the flaw (More in the documentation).

I just left the crt.bi in, the code is from the fbc source code and is a translation from the .c file in the run time library (rtl)
But #4 ,#5 are unaffected by a seed given via randomize in any case.
v1ctor the compiler writer must have wrote it way back.
We rarely hear from v1ctor now although FreeBASIC is his baby.
Here is the math_rnd.bas from the fb source code.

Code: Select all

 /' rnd# function '/

#include "fb.bi"
'namespace rtl

#if defined (HOST_WIN32)
   #include "windows.bi"
   #include "win\wincrypt.bi"
#elseif defined (HOST_LINUX)
   #include "crt\fcntl.bi"
   #include "crt\unistd.bi"
#endif
namespace rtl
#define RND_AUTO      0
#define RND_CRT         1
#define RND_FAST      2
#define RND_MTWIST      3
#define RND_QB         4
#define RND_REAL      5

#define INITIAL_SEED   327680

#define MAX_STATE      624
#define PERIOD         397

extern "C"
declare function hRnd_Startup cdecl ( n as single ) as double
declare function hRnd_CRT cdecl ( n as single ) as double
declare function hRnd_QB cdecl ( n as single ) as double
Dim shared rnd_func as function cdecl ( as single ) as double = @hRnd_Startup
dim shared as uint32_t iseed = INITIAL_SEED
dim shared as uint32_t state(0 to MAX_STATE-1)
dim shared as uint32_t ptr p = NULL
dim shared as double last_num = 0.0

function hRnd_Startup cdecl ( n as single ) as double   
   select case __fb_ctx.lang
      case FB_LANG_QB:
         rnd_func = @hRnd_QB
         iseed = INITIAL_SEED
      case FB_LANG_FB_FBLITE or FB_LANG_FB_DEPRECATED:
         rnd_func = @hRnd_CRT
      case else
         fb_Randomize( 0.0, 0 )
   end select
   return fb_Rnd( n )
end function

function hRnd_CRT cdecl ( n as single ) as double
   if ( n = 0.0 ) then
      return last_num
   end if
   /' return between 0 and 1 (but never 1) '/
   return cast(double,rand( )) * ( 1.0 / ( cast(double, RAND_MAX) + 1.0 ) )
end function

function hRnd_FAST cdecl ( n as single ) as double
   /' return between 0 and 1 (but never 1) '/
   /' Constants from 'Numerical recipes in C' chapter 7.1 '/
   if ( n <> 0.0 ) then
      iseed = ( ( 1664525 * iseed ) + 1013904223 )
   end if

   return cast(double, iseed) / cast(double,4294967296ULL)
end function

function hRnd_MTWIST cdecl ( n as single ) as double
   if ( n = 0.0 ) then
      return last_num
   end if
   
   dim as uint32_t i, v, xor_mask(0 to 1) = { 0, &h9908B0DF }

   if ( p = NULL ) then
      /' initialize state starting with an initial seed '/
      fb_Randomize( INITIAL_SEED, RND_MTWIST )
   end if

   if ( p >= state(0) + MAX_STATE ) then
      /' generate another array of 624 numbers '/
      for i = 0 to MAX_STATE - PERIOD
         v = ( state(i) and &h80000000 ) or ( state(i + 1) and &h7FFFFFFF )
         state(i) = state(i + PERIOD) xor ( v shr 1 ) xor xor_mask(v and &h1)
      next
      
      for j as long =i to MAX_STATE - 1
         v = ( state(i) and &h80000000 ) or ( state(i + 1) and &h7FFFFFFF )
         state(i) = state(i + PERIOD - MAX_STATE) xor ( v shr 1 ) xor xor_mask(v and &h1)
      next
      v = ( state(MAX_STATE - 1) and &h80000000 ) or ( state(0) and &h7FFFFFFF  )
      state(MAX_STATE - 1) = state(PERIOD - 1) xor ( v shr 1 ) xor xor_mask(v and &h1)
      p = @state(0)
   end if

   v = *p + 1
   v xor= ( v shr 11 )
   v xor= ( ( v shl 7 ) and &h9D2C5680 )
   v xor= ( ( v shl 15 ) and &hEFC60000 )
   v xor= ( v shr 18 )

   return cast(double, v) / cast(double, 4294967296ULL)
end function

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

#if defined (HOST_WIN32) or defined (HOST_LINUX)
function hGetRealRndNumber cdecl ( ) as uinteger
   union _number
      as uinteger i
      as ubyte b(sizeof(uinteger))
   end union
   
   dim number as _number
   number.i = 0

   #if defined (HOST_WIN32)
   dim as HCRYPTPROV provider = 0
   if ( CryptAcquireContext( @provider, NULL, 0, PROV_RSA_FULL, 0 ) = TRUE ) then
      if ( CryptGenRandom( provider, sizeof(number), @number.b(0) ) = FALSE ) then
         number.i = 0
      end if
      CryptReleaseContext( provider, 0 )
   end if

   #elseif defined (HOST_LINUX)
   dim as long urandom = open_( "/dev/urandom", O_RDONLY )
   if ( urandom <> -1 ) then
      if ( read_( urandom, @number.b(0), sizeof(number) ) <> sizeof(number) ) then
         number.i = 0
      end if
      close_( urandom )
   end if
   #endif

   return number.i
end function

function hRnd_REAL cdecl ( n as single ) as double
   static as uinteger count = 0
   static as uinteger v
   dim as double mtwist

   mtwist = hRnd_MTWIST(n)
   if ( (count mod 256) = 0 ) then
      count = 1

      /' get new random number '/
      v = hGetRealRndNumber( )
   else
      count += 1
   end if

   if ( v = 0 ) then
      return mtwist
   end if
   v *= mtwist

   v xor= (v shr 11)
   v xor= ((v shl 7) and &h9D2C5680)
   v xor= ((v shl 15) and &hEFC60000)
   v xor= (v shr 18)

   return cast(double, v) / cast(double, 4294967296ULL)
end function
#endif

function fb_Rnd FBCALL ( n as single ) as double
   last_num = rnd_func( n )
   return last_num
end function

sub fb_Randomize FBCALL ( seed as double, algorithm as long )
   dim as long i

   union _dtoi
      as double d
      as uint32_t i(0 to 1)
   end union
   
   dim dtoi as _dtoi

   if ( algorithm = RND_AUTO ) then
      select case __fb_ctx.lang
         case FB_LANG_QB:
            algorithm = RND_QB
         case FB_LANG_FB_FBLITE or FB_LANG_FB_DEPRECATED:
            algorithm = RND_CRT
         case FB_LANG_FB:
            algorithm = RND_MTWIST
      end select
   end if

   if ( seed = -1.0 ) then
      /' Take value of Timer to ensure a non-constant seed.  The seeding
      algorithms (with the exception of the QB one) take the long value
      of the seed, so make a value that will change more than once a second '/

      dtoi.d = fb_Timer( )
      seed = cast(double, (dtoi.i(0) xor dtoi.i(1)))
   end if

   select case algorithm
      case RND_CRT:
         rnd_func = @hRnd_CRT
         srand( cast(uinteger, seed) )
         rand( )
         
      case RND_FAST:
         rnd_func = @hRnd_FAST
         iseed = cast(uint32_t, seed)

      case RND_QB:
         rnd_func = @hRnd_QB
         dtoi.d = seed
         dim as uint32_t s = dtoi.i(1)
         s xor= ( s shr 16 )
         s = ( ( s and &hFFFF ) shl 8 ) or ( iseed and &hFF )
         iseed = s

      #if defined (HOST_WIN32) or defined (HOST_LINUX)
      case RND_REAL:
         rnd_func = @hRnd_REAL
         state(0) = cast(uinteger,seed)
         for i = 1 to MAX_STATE
            state(i) = ( state(i - 1) * 1664525 ) + 1013904223
         next
         p = @state(0) + MAX_STATE
      #endif
      case RND_MTWIST:
         rnd_func = @hRnd_MTWIST
         state(0) = cast(uinteger, seed)
         for i = 1 to MAX_STATE
            state(i) = ( state(i - 1) * 1664525 ) + 1013904223
         next
         p = @state(0) + MAX_STATE
   end select
end sub
end extern

end namespace
   

Which includes #4.
If you download the fb source code then much of what goes on will be revealed.
What's wrong with your chisquare2g .exe from powerbasic?
I very rarely run linked .exe files with win 10.
But my confidence was returning and I decided to run yours.
But the failure has set me back months.

Return to “General”

Who is online

Users browsing this forum: No registered users and 30 guests