Bernard Widynski's Middle Square Weyl Sequence RNG (MsWs)

[MrSwiss]

It is also much faster than dipping into double and multiplying.
The uniformity seems OK. Perhaps the randomness needs further testing.

Agreed, but only (IMO), if you can use the generated value unmodified.
Putting it into a range, creates the problems ...

[deltarho[1859]]

Here is a very easy way to test the validity of an integer range function instead of dividing the range into equal units and counting how many 'hits' the units get looking for a reasonable balance. If we have a uniform distribution then the average will be equal to (first + last)/2 for an infinite test but fairly close with a 'decent' size test.

Code: Select all

Function rnd_range (first As Longint, last As Longint) As LongInt
  Function = Rnd*(last - first + 1) - .5 + first
End Function
 
Randomize
 
Dim As Ulong i
Dim As Longint tot
 
For i = 1 To 10^8
  tot += rnd_range(-123,456)
Next
Print tot/10^8
 
Sleep

Examples:

Code: Select all

(-5,5)     =>   0    Progamatically I got 1.118e-005
(-10,-1)   =>  -5.5   -do- I got -5.49982073
(0,9)      =>   4.5   -do- I got 4.49988167
(1,10)     =>   5.5   -do- I got 5.50016871
(347,692)  => 519.5   -do- I got 519.50917841
(-123,456) => 166.5   -do- I got 166.48569356

That very quickly proved that '-(first>=0)' was actually wrong.

@dodicat

Still not entirely sure what you mean but will this do? It is simply the above adapted to take a random 64-bit value as a further parameter.

Code: Select all

Function rnd_range (RandomNumber as UlongInt, first As Longint, last As Longint) As LongInt
  Function = RandomNumber/2^64*(last - first + 1) - .5 + first
End Function

[dodicat]

similar results with no float operators, only uongint/longint operations.

Code: Select all

 
'==========================
type Rand
   a as ulongint
   b as ulongint
   c as ulongint
   d as ulongint
end type

function ranulong(x as rand) as ulongint
   dim e as ulongint = x.a - ((x.b shl 7) or (x.b shr (64 - 7)))
   x.a = x.b xor ((x.c shl 13) or (x.c shr (64 - 13)))
   x.b = x.c + ((x.d shl 37) or (x.d shr (64 - 37)))
   x.c = x.d + e
   x.d = e + x.a
   return x.d
end function

'function randouble(x as rand) as double ''NOT NEEDED
    'return ranulong(x)/18446744073709551615ull
   ' end function

 sub init(x as rand, byval seed as ulongint=100)
   dim i as ulongint
    x=type(4058668781,seed,seed,seed)
    for i as ulongint=1 to 20
        ranulong(x)
        next
end sub

'=========
dim shared as rand z
init(z,200)
'========


function range(f as longint,l as longint) as longint 'for any ranges, ulongint for + only ranges
    return (ranulong(z) mod ((l-f)+1)) + f
end function


Dim As Ulong i
Dim As Longint tot
 dim as double t=timer
For i = 1 To 10^8
  tot += range(-5,5)
Next
Print tot/10^8
tot=0
For i = 1 To 10^8
  tot += range(-10,-1)
Next
Print tot/10^8
tot=0
For i = 1 To 10^8
  tot += range(0,9)
Next
Print tot/10^8
tot=0
For i = 1 To 10^8
  tot += range(1,10)
Next
Print tot/10^8
tot=0
For i = 1 To 10^8
  tot += range(347,692)
Next
Print tot/10^8
tot=0
For i = 1 To 10^8
  tot += range(-123,456)
Next
Print tot/10^8
print "Time taken  ";timer-t
sleep
   

 

Code: Select all

 -0.00051191
-5.50043617
 4.49987626
 5.50025758
 519.49842485
 166.4962325
Time taken   2.242408916354179
 

[sean_vn]

Nice RNG algorithm. I must remember it.
There is a 128 bit multiply in the AMD 64 instruction set that you can use for fixed point calculations. A random 64 bit number by (n+1) will give a result between 0 and n in the register that contains the high 64 bits of the 128 bit result. Not perfect but generally okay.

[deltarho[1859]]

@dodicat

You are using a few rotations there and they remind me of Sebastiano Vigna's xorshiro algorithms where he used a C macro

Code: Select all

static inline uint64_t rotl(const uint64_t x, int k) {
  return (x << k) | (x >> (64 - k));
}

I still haven't got my head around

Code: Select all

Return (ranulong(z) Mod ((l-f)+1)) + f

yet but I will get there. <smile>

It does not look like it will compete against Knuth64 on speed but I bet it has better randomness than Knuth64 which leaves a lot to be desired.

[dodicat]

deltarho
The difference is ( last-first +1)
then first + random number mod difference puts the final number in the output range.
The randomising algo can be half decent ulongint generator.
The range function simply avoids doubles.
The mod operator works in the ulongint range, although the help file says integer.
You can overload for longint and ulongint.
pass parameters with ll or ull subscripts.
The line
return (ranulong(z) mod ((l-f)+1)) + f
is generally
return (someulongintgenerator(z) mod ((l-f)+1)) + f
You own generator would be very fast I think.

[deltarho[1859]]

You own generator would be very fast I think.

Blindingly fast but it failed the fast validity test. [1]

I used this in the opening post:

Code: Select all

Function MsWsParams.Range( First As Longint, Last As Longint ) As Longint
  This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
  This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
  Function = This.Seed/2^64 Mod ( Last - First + 1)  + First
End Function

Added: [1]. For the life of me, I cannot see how Mod can possibly work in this context.

[dodicat]

For your generator I use

Code: Select all


#define Get64Bit Cast( Ulongint, Rnd*2^64 )

Type MsWsParams
    Declare Constructor
    Declare Sub MyRandomize( Byval seed As Ulongint = 0, Byval sequence As Ulongint = 0 )
    Declare Function Rand() As Ulongint
    Declare Function RandSE() As Double
    Declare Function Range Overload ( First As Double, Last As Double ) As Double
    Declare Function Range Overload ( First As Longint, Last As Longint ) As Longint
    Declare Function Range Overload ( First As Ulongint, Last As Ulongint ) As Ulongint
    Seed As Ulongint
    Sequence As Ulongint
End Type

Sub MsWsParams.MyRandomize( Byval seed As Ulongint = 0, Byval sequence As Ulongint = 0 )
    Randomize , 5
    If seed = 0 Then
        This.Seed = Get64Bit
        This.Sequence = Get64Bit
    Else
        If sequence = 0 Then
            This.Seed = seed
            This.Sequence = Get64Bit
        Else ' For a sequence repeat
            This.Seed = seed
            This.Sequence = sequence
        End If
    End If
    ' Warm up generator - essential for any PRNG
    For i As Ulong = 1 To 1000
        This.rand()
    Next
End Sub

' 18446744073709551557 is the largest prime number < 2^64

Function MsWsParams.Rand() As Ulongint
    This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
    This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
    Return( This.Seed )
End Function

Function MsWsParams.RandSE() As Double
    This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
    This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
    Return This.Seed/2^64
End Function

Function MsWsParams.Range( First As Double, Last As Double ) As Double
    This.Seed *= This.Seed : This.Sequence += 18446744073709551557 : This.Seed += This.Sequence
    This.Seed = ( This.Seed Shr 32 ) Or ( This.Seed Shl 32 )
    Function = This.Seed/2^64*( Last - First ) + First
End Function

'alterations
Function MsWsParams.Range( First As Longint, Last As Longint ) As Longint
    Return (this.Rand() Mod ((last-first)+1)) + first
End Function

Function MsWsParams.Range( First As Ulongint, Last As Ulongint ) As Ulongint
    Return (this.Rand() Mod ((last-first)+1)) + first
End Function

Constructor MsWsParams
This.MyRandomize()
End Constructor

Dim  As MsWsParams z

Dim As Ulong i
Dim As Longint tot
Dim As Double t=Timer
For i = 1 To 10^8
    tot += z.range(-5ll,5)
Next
Print tot/10^8
tot=0
For i = 1 To 10^8
    tot += z.range(-10ll,-1)
Next
Print tot/10^8
tot=0
For i = 1 To 10^8
    tot += z.range(0ull,9)
Next
Print tot/10^8
tot=0
For i = 1 To 10^8
    tot += z.range(1ull,10)
Next
Print tot/10^8
tot=0
For i = 1 To 10^8
    tot += z.range(347ull,692)
Next
Print tot/10^8
tot=0
For i = 1 To 10^8
    tot += z.range(-123ll,456)
Next
Print tot/10^8

Print "pure ulongint example"

For i = 1 To 10
    Print z.range(18446744073709551610ull,18446744073709551613)
Next


Print "Time taken  ";Timer-t
Sleep




 

result

Code: Select all

-6.353e-005
-5.50021221
 4.50014817
 5.49993317
 519.49835457
 166.50732974
pure ulongint example
18446744073709551611
18446744073709551610
18446744073709551613
18446744073709551613
18446744073709551610
18446744073709551612
18446744073709551612
18446744073709551611
18446744073709551612
18446744073709551613
Time taken   2.048540160991252
 

[deltarho[1859]]

@dodicat

I was using This.Seed/2^64 when I should have used just This.Seed. Dear, oh dear. Your method is producing the same results as dafhi's and my Int method before that.

The reason why I used This.Seed, as opposed to This.Rand, is because the engine only takes up a few lines so I use the engine for every function to avoid calling another function. As it turns out the difference between the engine plus This.Seed and This.Rand is negligible. FB keeps surprising me.

Anyway, my original Int method, using '-gen gcc -Wc -O3' in 32-bit mode, was beaten by dafhi's method from 34MHz to 63MHz. Your method storms in at 135MHz.

The most extraordinary thing is using '-gen gcc -Wc -O3' in 64-bit mode saw my Int method and dafhi's method coming in at the same speed and now your method is the same. Goodness knows what the binary looks like and I find it hard to believe that all three methods have been reduced to the same native code. However, they are coming in at the same speed nonetheless.

Since you MODification has improved the 32-bit results I have revised the LongInt in the opening post yet again for folks who use 32-bit.

I must now sit down and fathom out how the Mod actually works. How did you come up with that? Is it from your archives or did you come up with it during the life of this thread? This is important because if you got it from your archives then you qualify for a fortnight in Clacton-On-sea. On the other hand, if you came up with it during the life of this thread you qualify for a week in Clacton-On-Sea. The burning question is are you going to tell 'porkies' to avoid going to Clacton for a fortnight.

I wonder what MrSwiss thinks about your Mod approach. We cannot get much more elegant than that and it is very close to being four times faster than my original Int method. That should definitely go into the manual but it will, I am sure, leave a lot of folks scratching their heads.

Actually, I am being unfair about Clacton because I have stayed in a few Bed & Breakfasts there and I have no complaints on that score. It is probably very good for families as there is plenty of stuff going on to keep the kids fully occupied. It is over 25 years since I was last there but there was a superb fish restaurant that I used to frequent owned by a divorced lady who I took a shine too. She knocked out a haddock to die for. <smile>

[deltarho[1859]]

Remembering that this thread is about MsWs it is worth noting that in 64-bit mode MsWs' range function is about 11% faster than PCG32II. However, for RandSE PCG32II is about 12% faster than MsWs. In 'real life' scenarios it may not be easy to separate them. There is a difference in that PCG32II is a 32-bit generator and MsWs is a 64-bit generator with PCG32II's RandSE having 32-bit granularity and MsWs' RandSE having 53-bit granularity. It could be argued then that MsWs has elbowed PCG32II off the podium.

[dodicat]

My ideal fortnight would be absolute silence and solitude where I can practice my medication techniques, thus I would choose Morecambe.

[MrSwiss]

dodicat, just make very shure to NOT run out of medicine or have a established supply-chain there.<lol>

[dafhi]

range using modulus may fit unevenly into 2^64. the remainder goes to the least significant digits

[deltarho[1859]]

Code: Select all

Dim As Longint lo = 20, hi = 147, R, ans
Dim As Double x

R = 7547257594037927936
ans =  R Mod (hi - lo + 1)  + 20 ' dodicat
Print ans
x = R/2^64*(hi - lo + 1) - 0.5 + 20
ans = R/2^64*(hi - lo + 1) - 0.5 + 20  ' dafhi
Print ans, x

Sleep

gives

Code: Select all

 20
 72            71.86961971048731

dafhi's ans is correct

Perhaps we need to write our own Mod function if it is fast enough.

[dafhi]

also, you could get rid of the hi (+ 1) if you change params to double

Perhaps we need to write our own Mod function if it is fast enough.

:-)

i'm having a blast with this RNG stuff