MsWsII PRNG plus Help file

General FreeBASIC programming questions.
hhr
Posts: 206
Joined: Nov 29, 2019 10:41

Re: MsWsII PRNG plus Help file

Post by hhr »

In my computer getduo does not work with gas64:

test.asm: Assembler messages:
test.asm:783: Error: unsupported instruction `mov'

I located the error in line 140: u.Output64 = Temp Xor This.Seed2

Changing the line to

u.Output64 = Temp
u.Output64 Xor= This.Seed2

eliminates the error message.

gcc and gas32 are OK.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: MsWsII PRNG plus Help file

Post by deltarho[1859] »

@hhr

Wow. That is something for SARG to think about.

@srvaldez

Your machine is lethal. :)

Bernard Widynski wrote: “If 32-bit precision is adequate (which is most often the case) msws64 may be the fastest RNG for producing floating-point numbers.”

Well, it doesn't look like it to me.

I cannot see any way to make GetDuo faster. I will double-check Engine2, but any errors in that should have shown up with all the tests that I have done.

It never ceases to amaze me how 32-bit and 64-bit differ. With both srvaldez and me GetDuo is very much the worse in 32-bit but pretty much neck and neck in 64-bit.

I think I will have to send Bernard Widynski a very diplomatic email asking how he managed the first entry in his Table 1. Yes, I can be diplomatic when needed. :)
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: MsWsII PRNG plus Help file

Post by dodicat »

deltarho[1859] wrote: Apr 20, 2022 10:36 Thanks, dodicat.

What I was particularly interested in was a run as is with msws.getduo(x,y) and then to comment that with a run using x = Rnd and y = Rnd.

Your CPU is Q3'11 so will probably see GetDuo as much a disappointment as with my machine. I am hoping to see GetDuo do better than 2 x Rnd to justify leaving it in MsWsII. Widynski's architecture is only three years old.
about 17 seconds (I have edited my post) with rnd.
Also
Error: unsupported instruction `mov'
with -gen gas64 (like hhr)
But after the fix proposed by hhr, 31 seconds with -gen gas64.
SARG
Posts: 1756
Joined: May 27, 2005 7:15
Location: FRANCE

Re: MsWsII PRNG plus Help file

Post by SARG »

@hhr, @dodicat
Which version of fbc are you using (regular, beta) ?
I was trying to reproduce the error but unsuccefully.
Maybe I did something wrong when compiling the modules. Please tell me how I should do.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: MsWsII PRNG plus Help file

Post by deltarho[1859] »

@dodicat

So GetDuo is faster on your machine than 2 x Rnd (randse).

This is weird. :?

Anyway, I have sent Bernard Widynski an email.

Re the Monte Carlo. In 64-bit mode and 64-bit output, MsWsII is nearly twice as fast as Sebastiano Vigna's xoroshiro128**; which is a 64-bit output generator.

Just goes to show that MHz throughput testing does not tell us the real story – real-world applications do. :)
hhr
Posts: 206
Joined: Nov 29, 2019 10:41

Re: MsWsII PRNG plus Help file

Post by hhr »

I tried with FreeBASIC-1.09.0-winlibs-gcc-9.3.0 and FreeBASIC-1.09.0-win64 in Windows 7.

Code: Select all

#cmdline "-gen gas64"
'#cmdline "-gen gcc"

Print __fb_version__
Print __fb_signature__
Print __fb_backend__
#IFDEF __FB_64BIT__
Print "64Bit"
#ELSE
Print "32Bit"
#ENDIF

#Include "MsWsII.bas"

Dim As Double x,y
msws.getduo(x,y)
Print x,y

Sleep
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: MsWsII PRNG plus Help file

Post by dodicat »

Scrub all that deltarho.
I was careless again.
Here are the true results from here

Code: Select all



 6.566717800000333 seconds
Estimate of Pi: 3.141524532  with 64 bit fbc using  GetDuo

6.483655300005921 seconds  
Estimate of Pi: 3.141587388  with 64 bit fbc using rnd



 82.1120675000036 seconds
Estimate of Pi: 3.14162766  with 32 bit fbc using  GetDuo

 89.94279500000722 seconds
Estimate of Pi: 3.141613888 with 32 bit fbc using rnd


 
SARG
I am still using the official version 1.09.0, same as hhr.
srvaldez
Posts: 3373
Joined: Sep 25, 2005 21:54

Re: MsWsII PRNG plus Help file

Post by srvaldez »

hello SARG :)
here's an all-in-one for you, the problem is in line ~141

Code: Select all

' MsWsII.bas Version 1.0

#Include "Mysvalues.bi"

#define dodrange(f,l) Int(Rnd*(((l)+1)-(f))+(f)) ' By dodicat
#ifndef Uint32
  #define Uint32 Uinteger<32>
#endif
#ifndef Uint64
  #define Uint64 Uinteger<64>
#endif
#ifndef Int32
  #define Int32 Integer<32>
#endif

#macro Engine1
  This.Seed *= This.Seed : This.Sequence += svalue(0) : This.Seed += This.Sequence
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
#endmacro

#macro Engine2
  This.Seed1 *= This.Seed1
  This.Sequence1 += This.svalue(1)
  This.Seed1 += This.Sequence1
  Temp = This.Seed1
  This.Seed1 = ( This.Seed1 Shr 32 ) Or ( This.Seed1 Shl 32 )
  This.Seed2 *= This.Seed2
  This.Sequence2 += This.svalue(2)
  This.Seed2 += This.Sequence2
  This.Seed2 = ( This.Seed2 Shr 32 ) Or ( This.Seed2 Shl 32 )
#endmacro

Type MsWsParams
  Public:
  Declare Constructor
  Declare Sub Initialize()
  Declare Function rand() As Uint32
  Declare Function randse() As Double
  Declare Function randd() As Double
  Declare Sub GetDuo( Byref as Double, Byref as Double )
  Declare function range overload ( First as double, Last as double ) as double
  Declare function range overload ( First as Int32, Last as Int32 ) as Int32
  Declare Sub GetEngine1Snapshot()
  Declare Sub SetEngine1SnapShot()
  Declare Sub GetEngine2Snapshot()
  Declare Sub SetEngine2SnapShot()
  Declare Sub RestoreEngine1InitialState()
  Declare Sub RestoreEngine2InitialState()
  Declare Sub RestoreAllInitialState()
  Declare Function Gauss() as double
  'Declare Sub DumpState() ' Used during development
  Private:
  Declare Sub GetInitialState()
  As Uint64 svalue(0 to 2)
  As Uint64 SnapStateVector(0 to 5)
  As Uint64 InitialState(0 To 5)
  ' State vectors
  Seed As Uint64
  Seed1 As Uint64
  Seed2 As Uint64
  Sequence As Uint64
  Sequence1 As Uint64
  Sequence2 As Uint64
End Type

' ***********************************
' Seeding functions
Function HammingSeed64() As Uint64
Dim As Uint32 tscSeed0, tscSeed1, numBits
Dim As Uint64 Seed, CopySeed
  
  While Not ((numBits > 30) And (numBits < 34)) ' Gd quality entropy - we don't need exactly 32
    Asm
      rdtsc
      mov ecx, eax
      bswap eax ' fast moving to upper bits
      mov Dword Ptr [tscSeed0], eax
      add ecx, 1073741824 ' A quarter of a spin
      bswap ecx ' fast moving to upper bits
      mov Dword Ptr [tscSeed1], ecx
    End Asm
    Seed = (Cast( Uint64, tscSeed0 ) Shl 32) Or Cast( Uint64, tscSeed1 ) 
    CopySeed = Seed : numBits = 0
    While CopySeed <> 0
      CopySeed = CopySeed And CopySeed - 1
      numBits += 1
    Wend
  Wend
  Return Seed
End Function

Sub ShuffleUint64( Byref x As Uint64 )
' dodrange uses 'Randomize , 5' as set in the Constructor
  Union localUDT
    As Uint64 ul
    As Byte b(7)
  End Union
  Dim As localUDT l
  l.ul = x
  For i As Uint32 = 0 to 6
    Swap l.b(i), l.b(dodrange(i, 7))
  Next
  x = l.ul ' <- without this x will not change
End Sub
' ***********************************

' 32-bit random number
Function MsWsParams.rand() As Uint32
  Engine1
  Return This.Seed
End Function

' 32-bit granularity [0,1)
Private Function MsWsParams.randse() As Double
Dim As Uint32 Temp
  
  Engine1
  Temp = This.Seed
  Return Temp/2^32
End Function

' 53-bit granularity [0,1)
   Private Function MsWsParams.randd() As Double
   Dim As Uint64 Temp
      
      Engine2
      Return ((Temp xor This.Seed2) shr 11)/2^53
   end function

' Returns two 32-bit granularity [0,1)
Private Sub MsWsParams.GetDuo( ByRef One As Double, ByRef Two As Double )
Dim As Uint64 Temp
 
  Union GetUint32
    OutPut64 As Uint64
    Output32(0 To 1) As Uint32
  End Union
  Dim As GetUint32 u
  Engine2
  u.Output64 = Temp Xor This.Seed2
  'u.Output64 Xor= This.Seed2
  One = u.Output32(0)/2^32
  Two = u.OutPut32(1)/2^32
End Sub

' Floating point range
Private function MsWsParams.range( First as double, Last As Double ) As Double
Dim As Uint32 Temp
  
  Engine1
  Temp = This.Seed
  Function = Temp/2^32 * ( Last - First ) + First 
End Function

' Integral range - Int32 (Long)
Private Function MsWsParams.range(First As Int32, Last As Int32) As Int32
Dim As Uint32 Temp
  
  Engine1
  Temp = This.Seed
  Return clng( Temp Mod (Last-First+1)) + First
End Function

' Normal distribution
Private Function MsWsParams.Gauss As Double
Static As Int32 u2_cached
Static As Double u1, u2, x1, x2, w
  
  If u2_cached = True Then
    u2_cached = False
    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 = True
    Function = u1
  End If
End Function

Private Sub MsWsParams.GetEngine1Snapshot( )
  snapStateVector(0) = Seed
  snapStateVector(3) = Sequence
End Sub

Private Sub MsWsParams.SetEngine1Snapshot( )
  If snapStateVector(0)=0 Andalso snapStateVector(3)=0 then
    ' ie a GetSnapshot has not been executed so use InitialState()
    Seed = InitialState(0)
    Sequence = InitialState(3)
  Else
    Seed = SnapStateVector(0)
    Sequence = SnapStateVector(3)
  End If
End Sub

Private Sub MsWsParams.GetEngine2Snapshot( )
  snapStateVector(1) = Seed1
  snapStateVector(2) = Seed2
  snapStateVector(4) = Sequence1
  snapStateVector(5) = Sequence2
End Sub

Private Sub MsWsParams.SetEngine2Snapshot( )
  If snapStateVector(1)=0  Andalso snapStateVector(2)=0 Andalso snapStateVector(4)=0  Andalso snapStateVector(5)=0  then
    Seed1 = InitialState(1)
    Seed2 = InitialState(2)
    Sequence1 = InitialState(4)
    Sequence2 = InitialState(5)
  Else
    Seed1 = SnapStateVector(1)
    Seed2 = SnapStateVector(2)
    Sequence1 = SnapStateVector(4)
    Sequence2 = SnapStateVector(5)
  End If
End Sub

Sub MsWsParams.GetInitialState( )
  InitialState(0) = Seed
  InitialState(1) = Seed1
  InitialState(2) = Seed2
  InitialState(3) = Sequence
  InitialState(4) = Sequence1
  InitialState(5) = Sequence2
End Sub

Private Sub MsWsParams.RestoreEngine1InitialState( )
  Seed = InitialState(0)
  Sequence = InitialState(3)
End Sub

Private Sub MsWsParams.RestoreEngine2InitialState( )
  Seed1 = InitialState(1)
  Seed2 = InitialState(2)
  Sequence1 = InitialState(4)
  Sequence2 = InitialState(5)
End Sub

Private Sub MsWsParams.RestoreAllInitialState()
    RestoreEngine1InitialState
    RestoreEngine2InitialState
End sub

' Random seeding for all
Sub MsWsParams.Initialize( )
  ' For Engine1
  Seed = HammingSeed64 : ShuffleUint64( Seed )
  ' For Engine2
  Seed1  = HammingSeed64 : ShuffleUint64( Seed1 )
  Seed2 = HammingSeed64 : ShuffleUint64( Seed2 )
  ' For Engine1
  Sequence = HammingSeed64 : ShuffleUint64( Sequence )
  ' For Engine2
  Sequence1 = HammingSeed64 : ShuffleUint64( Sequence1 )
  Sequence2 = HammingSeed64 : ShuffleUint64( Sequence2 )
End Sub

Constructor MsWsParams
Dim as Uint32 ubnd, first, second, third
  
  Randomize , 5 ' Cryptographic and is a global setting
  ubnd = UBound( Mysvalues )
  first = dodrange(0,ubnd)
  do
    second = dodrange(0,ubnd)
  Loop until second <> first
  Do
    third = dodrange(0,ubnd)
  Loop until third <> first andalso third <> second
  This.svalue(0) = Mysvalues(first)
  This.svalue(1) = Mysvalues(second)
  This.svalue(2) = Mysvalues(third)
  Initialize
  GetInitialState
  Randomize , 3  ' This avoids 'Randomize , 5' being used in the following modules where you may want to revise it.
End Constructor

' Used during development 
'Private Sub MsWsParams.DumpState
'   Print Hex(This.Seed,16);" ";Bin(This.Seed,64)
'   Print Hex(This.Seed1,16);" ";Bin(This.Seed1,64)
'   Print Hex(This.Seed2,16);" ";Bin(This.Seed2,64)
'   Print Hex(This.Sequence,16);" ";Bin(This.Sequence,64)
'   Print Hex(This.Sequence1,16);" ";Bin(This.Sequence1,64)
'   Print Hex(This.Sequence2,16);" ";Bin(This.Sequence2,64)
'End Sub

Dim Shared As MsWsParams msws
#undef Rnd
#define Rnd msws.randse
#define RndD msws.randd
#define Range_ msws.range


'''#Include "MsWsII.bas"
 
Dim As Uint32 n = 10^9, i, lCount
Dim As Double t, x, y, z
 
lCount = 0
t = Timer
For i = 1 To n
  msws.getduo(x,y)
'  x = Rnd
'  y = Rnd
  If Sqr( x*x + y*y) <= 1 Then lCount += 1
Next
Print Timer - t;" seconds"
z = 4*lCount/n
Print "Estimate of Pi:";z
Sleep
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: MsWsII PRNG plus Help file

Post by deltarho[1859] »

Bernard Widynski got back to me – that was quick.

He pointed me to his timing code in his software download, which I hadn't spotted.

Code: Select all

#define two32 4294967296.
 
#include "msws64.h"
#include <stdint.h>
#include <stdio.h>
 
#define billion 1000000000
 
int main(void)
{
    int32_t i,n = billion, m = billion/2;
    double sum = 0;
 
    union {
       uint64_t i64;
       uint32_t i32[2];
    } u;
 
    for (i=0; i<m; i++) {
       u.i64 = msws64();
       sum += (u.i32[0] / two32);
       sum += (u.i32[1] / two32);
    }
 
    printf ("average = %lf\n",sum/n);
 
    return 0;
}
So, his timing loop is in a function.

In the Monte Carlo code, GetDuo is in a loop, so we have an accumulation of function overhead.

Thinking cap time. :D
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: MsWsII PRNG plus Help file

Post by dodicat »

...
Last edited by dodicat on Apr 21, 2022 8:19, edited 2 times in total.
SARG
Posts: 1756
Joined: May 27, 2005 7:15
Location: FRANCE

Re: MsWsII PRNG plus Help file

Post by SARG »

@hhr,dodicat and srvaldez
Thanks for the report and the help.

The bug is fixed. A case, due to the union structure (GetUint32), badly handled in the optimization function.
On my side the error didn't appear as the used -g option added an extra label (marking each basic line) which avoided the faulty optimization.

gas64 = 40.32s /gcc64 = 43.71s / gcc64 with O2 = 4.78s :shock:

@deltarho
Sorry for this part out of subject. :)
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: MsWsII PRNG plus Help file

Post by deltarho[1859] »

SARG wrote:@deltarho
Sorry for this part out of subject. :)
My code went down using gas64 and was reported by hhr. Before the day was out, you had fixed the bug. Great stuff SARG.

@dodicat

At Widynski's website, here, in section 8 '64-bit Output' a 64-bit engine is described. In MsWsII.bas this has been designated Engine2.

With Engine1, we obtain the lower half of the middle from the original middle square. This is done by using the lower 32 bits of the last 64-bit calculation. In randse This.Seed is 64 bit, and we equate that with 'Temp = This.Seed' where Temp is 32 bit.

At first glance, Engine2 looks like we are XORing two Engine1 lookalikes, but this would only give us 32 bits. There is a twist. That twist is in the statement 'Temp = This.Seed1', where in this case Temp is 64-bit, and it is used in 'Temp Xor This.Seed2'. None of that is discarded. This is first used in RndD (randd) where a 53-bit granularity [0,1) is returned.

Regarding your question: “Do You have a 64-bit integer random generator” the answer is yes; RomuTrio, xoroshiro128**, xoshiro256** and now we have another with Engine2 in MsWsII.

Added: I am going to ignore the code you posted as it is no longer relevant having answered your questions. Of course, it was an opportunity to promote your own 64-bit generator, which you should not be doing in this thread.

@SARG

When we have code which uses a fair amount of bit manipulation, as do many PRNGs, gcc does an extraordinary job of optimization. If the arithmetic is heavily biased toward 64-bit unsigned integers then, in 64-bit mode, the roof comes off from a performance perspective. Without optimization, my PRNG implementations are pedestrian. With -O2 we need to wear a seatbelt. With my Encrypternet program, from a performance perspective, there is no significant difference between a gas64 compilation and a gcc -O2 compilation. Sometimes gas64 binaries are a little faster, and sometimes gcc -O2 binaries are a little faster. The gas64 binaries are larger, but not enough to present us with grounds for concern.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: MsWsII PRNG plus Help file

Post by dodicat »

Actually deltarho, you posted the code in C.
I merely translated it to fb.
I couldn't find a 64 bit output to test it with among your methods.
engine2 in MsWsII.bas yes, but you have it in a macro.
I don't want to promote my 64 bit engine because it is pretty mediocre speed wise, but why don't you include a 64 bit output in MsWsII.bas.
I have removed my post using the translated c code thing, it was only a tester anyway, but it absolutely needs a 64 bit output to work properly.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: MsWsII PRNG plus Help file

Post by deltarho[1859] »

@dodicat

OK. I have added rand64 which outputs an Uint64.

The average of 10^9, with one run, came to 9.223382583316027e+018 which is accurate to nearly 6 significant figures.

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

Re: MsWsII PRNG plus Help file

Post by deltarho[1859] »

Back to the GetDuo saga.

I won't bore you with the details, but gcc was screwing up the timings. I may not be able to do this for all occasions, but I managed to stop gcc from doing that. The upshot is there is no difference between GetDuo and 2 x Rnds in 64-bit. In 32-bit GetDuo is slower.

We have then a 'hiding to nothing' scenario, so I have removed GetDuo from the package and revised the Help file.

Rather than simply remove GetDuo I suggest that you copy the whole MsWsII.bas source because I have done some tidying up to improve readability.

The chm zip in the opening post has been renamed to MsWsIIHelp.zip
Post Reply