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

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

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

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

and the -0.5

Double has been taken care of.

In any case, we should focus on one thing at a time. Mod is the issue.
deltarho[1859]
Posts: 2186
Joined: Jan 02, 2017 0:34
Location: UK

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

dodicat wrote:The mod operator works in the ulongint range, although the help file says integer.

'R Mod n' returns zero if R >= 2^32-1.

If we write our own Mod function it looks like we need extended precision, ie 80-bit, but we don't have it. We should have since the FPU uses 80-bit. PowerBASIC has extended precision.
MrSwiss
Posts: 3348
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

deltarho[1859] wrote:'R Mod n' returns zero if R >= 2^32-1.
This sort of statement is the reason, that I don't trust code of yours.

Proof that above statement is nonsense (tested FBC 32/64, 1.05.0):

Code: Select all

Const As LongInt    pllmax = &h7FFFFFFFFFFFFFFFll, _    ' positive 64 bit maximum
nllmax = &hFFFFFFFFFFFFFFFFll       ' negative 64 bit maximum

Print nllmax Mod 19
Print pllmax Mod 19

Print : Print "done ... ";

Sleep
result wrote:-1
17
Added: the neg. max. referres to binary = set bit's (value -1)
dodicat
Posts: 6155
Joined: Jan 10, 2006 20:30
Location: Scotland

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

mod seems OK for > 2^32-1
Why do you reckon it is not?

Here are the string mods and fb mods for ulongints high in the range:

Code: Select all

Function _divide(n1 As String,n2 As String,decimal_places As integer=10,dpflag As String="s") As String
Dim As String number=n1,divisor=n2
dpflag=lcase(dpflag)
'For MOD
dim as integer modstop
if dpflag="mod" then
if len(n1)<len(n2) then return n1
if len(n1)=len(n2) then
if n1<n2 then return n1
end if
modstop=len(n1)-len(n2)+1
end if
if dpflag<>"mod" then
If dpflag<>"s"  Then dpflag="raw"
end if
Dim runcount As integer
'_______  LOOK UP TABLES ______________
Dim Qmod(0 To 19) As Ubyte
Dim bool(0 To 19) As Ubyte
For z As Integer=0 To 19
Qmod(z)=(z Mod 10+48)
bool(z)=(-(10>z))
Next z

'_______ SET THE DECIMAL WHERE IT SHOULD BE AT _______
Dim As String part1,part2
#macro set(decimal)
#macro insert(s,char,position)
If position > 0 And position <=Len(s) Then
part1=Mid\$(s,1,position-1)
part2=Mid\$(s,position)
s=part1+char+part2
End if
#endmacro
If dpflag="raw" Then
End if
#endmacro
'______________________________________________
'__________ SPLIT A STRING ABOUT A CHARACTRR __________
Dim As String var1,var2
Dim pst As integer
#macro split(stri,char,var1,var2)
pst=Instr(stri,char)
var1="":var2=""
If pst<>0 Then
var1=Rtrim(Mid(stri,1,pst),".")
var2=Ltrim(Mid(stri,pst),".")
Else
var1=stri
End if
#endmacro

#macro Removepoint(s)
split(s,".",var1,var2)
#endmacro
'__________ GET THE SIGN AND CLEAR THE -ve __________________
Dim sign As String
If Left(number,1)="-" Xor Left (divisor,1)="-" Then sign="-"
If Left(number,1)="-" Then  number=Ltrim(number,"-")
If Left (divisor,1)="-" Then divisor=Ltrim(divisor,"-")

'DETERMINE THE DECIMAL POSITION BEFORE THE DIVISION
Dim As integer lennint,lenddec,lend,lenn,difflen
split(number,".",var1,var2)
lennint=Len(var1)
split(divisor,".",var1,var2)
lenddec=Len(var2)

If Instr(number,".") Then
Removepoint(number)
number=var1+var2
End if
If Instr(divisor,".") Then
Removepoint(divisor)
divisor=var1+var2
End if
Dim As integer numzeros
numzeros=Len(number)
number=Ltrim(number,"0"):divisor=Ltrim (divisor,"0")
numzeros=numzeros-Len(number)
lend=Len(divisor):lenn=Len(number)
If lend>lenn Then difflen=lend-lenn
Dim decpos As integer=lenddec+lennint-lend+2-numzeros 'THE POSITION INDICATOR
Dim _sgn As Byte=-Sgn(decpos)
If _sgn=0 Then _sgn=1
Dim As String thepoint=String(_sgn,".") 'DECIMAL AT START (IF)
Dim As String zeros=String(-decpos+1,"0")'ZEROS AT START (IF) e.g. .0009
if dpflag<>"mod" then
If Len(zeros) =0 Then dpflag="s"
end if
Dim As integer runlength
If Len(zeros) Then
runlength=decimal_places
If dpflag="raw" Then
runlength=1
If decimal_places>Len(zeros) Then
runlength=runlength+(decimal_places-Len(zeros))
End If
End If

Else
decimal_places=decimal_places+decpos
runlength=decimal_places
End if
'___________DECIMAL POSITION DETERMINED  _____________

'SET UP THE VARIABLES AND START UP CONDITIONS
number=number+String(difflen+decimal_places,"0")
Dim count As integer
Dim temp As String
Dim copytemp As String
Dim topstring As String
Dim copytopstring As String
Dim As integer lenf,lens
Dim As Ubyte takeaway,subtractcarry
Dim As integer n3,diff
If Ltrim(divisor,"0")="" Then Return "Error :division by zero"
lens=Len(divisor)
topstring=Left(number,lend)
copytopstring=topstring
Do
count=0
Do
count=count+1
copytemp=temp

Do
'___________________ QUICK SUBTRACTION loop _________________

lenf=Len(topstring)
If  lens<lenf=0 Then 'not
If Lens>lenf Then
temp= "done"
Exit Do
End if
If divisor>topstring Then
temp= "done"
Exit Do
End if
End if

diff=lenf-lens
temp=topstring
subtractcarry=0

For n3=lenf-1 To diff Step -1
takeaway= topstring[n3]-divisor[n3-diff]+10-subtractcarry
temp[n3]=Qmod(takeaway)
subtractcarry=bool(takeaway)
Next n3
If subtractcarry=0 Then Exit Do
If n3=-1 Then Exit Do
For n3=n3 To 0 Step -1
takeaway= topstring[n3]-38-subtractcarry
temp[n3]=Qmod(takeaway)
subtractcarry=bool(takeaway)
if subtractcarry=0 then exit do
Next n3
Exit Do

Loop 'single run
temp=Ltrim(temp,"0")
If temp="" Then temp= "0"
topstring=temp
Loop Until temp="done"
' INDIVIDUAL CHARACTERS CARVED OFF ________________
runcount=runcount+1
If count=1 Then
topstring=copytopstring+Mid(number,lend+runcount,1)
Else
topstring=copytemp+Mid(number,lend+runcount,1)
End If
copytopstring=topstring
topstring=Ltrim(topstring,"0")
if dpflag="mod" then
if runcount=modstop then
if topstring="" then return "0"
return mid(topstring,1,len(topstring)-1)
end if
end if
If topstring="" And runcount>Len(n1)+1 Then
Exit Do
End if
Loop Until runcount=runlength+1

' END OF RUN TO REQUIRED DECIMAL PLACES
set(decimal) 'PUT IN THE DECIMAL POINT
'THERE IS ALWAYS A DECIMAL POINT SOMEWHERE IN THE ANSWER
'NOW GET RID OF IT IF IT IS REDUNDANT
End Function

#define mod_(a,b) _divide(a,b,,"mod")
#define div_(a,b) _divide(a,b,-2)

dim as ulongint R=2^32-1   + 1000000
dim as ulongint n=234
print R Mod n' returns zero if R >= 2^32-1.

dim as ulongint u=-1

for n as ulongint=u-50 to u-1
dim as ulongint z=rnd*u
print n; " mod ";z;"  =  "
print "fb",,"string"
print n mod z,mod_(str(n),str(z))
print
next
sleep

Regarding the previuos code, the mod method preserves some granularity.

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

print
print
print "Mod","Double"
for z as longint=R-5 to R + 130
ans =  z Mod (hi - lo + 1)  + lo :print ans ,"   ";
ans = z/2^64*(hi - lo + 1) - 0.5 + lo :print ans
print

next

Sleep

The question I ask is not the integrity of mod but:
Both methods put the output in the range uniformly as previously tested.
Which method is more random.
(Tested 64 and 32 bit official fb)
MrSwiss
Posts: 3348
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

dodicat wrote:Both methods put the output in the range uniformly as previously tested. Which method is more random.
Not shure to understand properly.
The randomness should come from the gererator (used), not from the ranging
methods process ... ? (as long, as the generated value itself, is used 'as-is')
deltarho[1859]
Posts: 2186
Joined: Jan 02, 2017 0:34
Location: UK

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

@MrSwiss

I wrote a code snippet and tested it with various values. What surprised me was R>=2^32-1. I would have thought R>2^32-1.

Obviously, I did something wrong but I used 'Quick Run' on the snippet and am not able to see what I did wrong.

Now, I go to a lot of effort when posting on this forum to make sure that what I say or code is accurate. The following statement is totally out of order.
This sort of statement is the reason, that I don't trust code of yours.

You remind me of some people who say "You are always doing ....." or "You never do ...." when referring to an event which may never have happened before.
dodicat wrote:mod seems OK for > 2^32-1
Why do you reckon it is not?

That is how you should have responded.

I should tell you that I have just reported your statement to admin. I will say no more until I get a response.
MrSwiss
Posts: 3348
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

@deltarho, this must be a joke of yours.

Because, I stayed with it (even if 'sworn in'), by any court of law ...
I simply can't help it, if you can't take the truth.
That's what it is, no more, no less.
deltarho[1859]
Posts: 2186
Joined: Jan 02, 2017 0:34
Location: UK

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

MrSwiss
Posts: 3348
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

Fine with me. No problem.
dodicat
Posts: 6155
Joined: Jan 10, 2006 20:30
Location: Scotland

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

Anyway
My question was regarding the randomness.
for simplicity the out values are 0 to 10
Using Dafhi's scenario, it is equivalent to a mapping of all of 0 to 18446744073709551615 into 0 to 10.
So as a number ranges from 0 to 18446744073709551615 another number ranges from 0 to 10.
For many values of a number starting at 0, dafhi's output remains the same until the number passes the first 1/10 th. part of the ulongint range, then the 2/10 th. part e.t.c.

With the mod method every number has a different output as it ranges upwards.
0 to 10 then 0 to 10 ... ....
But because the outputs are 0 to 10 by each method they seem similar, but they couldn't be more different.
Which method produces the best randomness?

Roughly:

Code: Select all

function map(a as double,b as double,x as double,c as double,d as double) as double
return ((d)-(c))*((x)-(a))/((b)-(a))+(c)
end function

dim as ulongint lo=0,hi=10,y,y2,y3
dim as ulongint R=7547257594037927936
y=(R/2^64)*(hi-lo+1)-.5 +lo
print y;" by dafhi"
y= R Mod (hi - lo + 1)  + lo
print y; " by mod"
y= map(0,2^64-1,R,lo,hi)
print y;" by mapping"
print
print "press a key"
sleep
print

r=0

print "R";tab(30);"dafhi";tab(40);"map";tab(50);"mod"
do
y=(R/2^64)*(hi-lo+1)-.5 +lo
y2=map(0,2^64-1,R,lo,hi)
y3=R Mod (hi - lo + 1)  + lo

print r;tab(30);y;tab(40);y2;tab(50);y3
r+=  18446744073709551615\10
loop until r> 18446744073709551610 or y=10
sleep

The mod method happens to produce the same outputs, I wasn't really expecting that.
deltarho[1859]
Posts: 2186
Joined: Jan 02, 2017 0:34
Location: UK

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

@dodicat

Formula wise the only difference between dafhi and mapping is the divisor being '2^64' and '2^64-1' which will not be seen computationally.

Expanding 'function map' the ordering of computation is different between dafhi and mapping and the difference of 'R' compared with 'R-a' will be negligible for most cases. So, I have dropped mapping and concentrated on dafhi and Mod.

I have changed the divisions from hard-code to '\(hi-lo) and the loop limit from hard-code to 'y = hi'.

Code: Select all

Dim As Ulongint lo=50,hi=100,y,y2,y3
Dim As Ulongint R=7547257594037927936
y=(R/2^64)*(hi-lo+1)-.5 +lo
Print y;" by dafhi"
y= R Mod (hi - lo + 1)  + lo
Print y; " by Mod"
Print
Print "press a key"
Sleep
Print

r=0

Print "R";Tab(30);"dafhi";Tab(50);"Mod"
Do
y=(R/2^64)*(hi-lo+1)-.5 +lo
y3=R Mod (hi - lo + 1)  + lo
Print r;Tab(30);y;Tab(50);y3
r+=  18446744073709551615\(hi-lo)            ' \(hi-lo)
Loop Until r> 18446744073709551610 Or y=hi     ' y=hi
Sleep

With 'lo = 0' and 'hi = 10' we get as above.

However, with 'lo = 50' and 'hi = 100' dafhi has a clear mapping correspondence but Mod is all over the place. Small ranges hide the difference. A 'lo = 20' and 'hi = 50' sees Mod locking in at 20 for all r which makes me suspicious of my alterations.

Code: Select all

70 by dafhi
58 by Mod

press a key

R                            dafhi               Mod
0                            50                  50
368934881474191032           51                  65
737869762948382064           52                  80
1106804644422573096          53                  95
1475739525896764128          54                  59
1844674407370955160          55                  74
2213609288845146192          56                  89
2582544170319337224          57                  53
2951479051793528256          58                  68
3320413933267719288          59                  83
3689348814741910320          60                  98
4058283696216101352          61                  62
4427218577690292384          62                  77
4796153459164483416          63                  92
5165088340638674448          64                  56
5534023222112865480          65                  71
5902958103587056512          66                  86
6271892985061247544          67                  50
6640827866535438576          68                  65
7009762748009629608          69                  80
7378697629483820640          70                  95
7747632510958011672          71                  59
8116567392432202704          72                  74
8485502273906393736          73                  89
8854437155380584768          74                  53
9223372036854775800          75                  68
9592306918328966832          76                  83
9961241799803157864          77                  98
10330176681277348896         78                  62
10699111562751539928         79                  77
11068046444225730960         80                  92
11436981325699921992         81                  56
11805916207174113024         82                  71
12174851088648304056         83                  86
12543785970122495088         84                  50
12912720851596686120         85                  65
13281655733070877152         86                  80
13650590614545068184         87                  95
14019525496019259216         88                  59
14388460377493450248         89                  74
14757395258967641280         90                  89
15126330140441832312         91                  53
15495265021916023344         92                  68
15864199903390214376         93                  83
16233134784864405408         94                  98
16602069666338596440         95                  62
16971004547812787472         96                  77
17339939429286978504         97                  92
17708874310761169536         98                  56
18077809192235360568         99                  71
18446744073709551600         100                 86
dodicat
Posts: 6155
Joined: Jan 10, 2006 20:30
Location: Scotland

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

Thanks deltarho.
You would expect mod to be all over the place in almost any range and ulongint split accordingly, I was surprised to see it otherwise in 0 to 10.
Also range 20,30 and 40 and some others see mod a constant regurgitation.
It is all a bit like dafhi's method being akin to a large shell and the mod method a gatling gun.
Both causing equal carnage in different ways.
Anyway, I'll leave it at that for now.
deltarho[1859]
Posts: 2186
Joined: Jan 02, 2017 0:34
Location: UK

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

@dodicat

It is amazing how one piece of code tells us that we are on a winner and another tells us that all hell has broken loose. I think that we can agree that dafhi's method is more robust. I had already changed the opening post back to dafhi's method. I should tell you that there are issues with linear mapping in general and it is often the case of employing damage control.

I will also leave it a that.
dafhi
Posts: 1324
Joined: Jun 04, 2005 9:51

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

i'm testing the mod method on a ubyte LCG. just need "a few hours" debugging my test.

i kind of already know the result but .. gotta code it

...

Boom.

Code: Select all

function LCG as ubyte:  const mul = 23, add = 3
static as ushort state:  state = mul * state + add
return state shr 8
End Function

function range(hi as single=1, lo as single=0)as ubyte
return lcg mod (hi+1-lo) + lo
End Function

dim as long bins(255)

for i as long = 1 to 256 * 100
bins( lcg ) += 1
Next

var vlo = bins(0), vhi = vlo
for i as long = 1 to 255
if vlo > bins(i) then vlo = bins(i)
if vhi < bins(i) then vhi = bins(i)
bins(i) = 0
Next: bins(0) = 0

? "low count", "high count";  "  ...  normal range"
? vlo, vhi
?
?
? "halves sum, from range midpoint"
?
? "low", "high", "midpoint"
for range_hi as long = 10 to 30 step 2

for i as long = 1 to 256 * 19
bins( range( range_hi ) ) += 1
Next

dim as ulong m = (range_hi) \ 2
for i as long = 1 to m
bins(0) += bins(i)
bins(range_hi) += bins(range_hi-i)
Next

? bins(0), bins(range_hi), m

for i as long = 0 to range_hi
bins(i)=0
next

next:  sleep
dodicat
Posts: 6155
Joined: Jan 10, 2006 20:30
Location: Scotland

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

Hi dafhi.
Yes, it makes sense to use ubyte.
18446744073709551615 is really beyond comprehension for testing.
But I see you have range hi to lo, was that a deliberate trap, it caught me out at first anyhow (but that is not difficult)
Here is the crt rand(), the top half of short on this box.
There seems to be little difference in using the double method and the mod method, only speed.
We can only assume that rand() has been tried and tested for randomness.

Code: Select all

#include "crt.bi"

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

function chisquared(a() as double,R as range,byref m as double) as double
m=mean(a(),R)
dim as double acc
for n as long=lbound(a) to ubound(a)
acc+=(a(n)-m)^2/m
next
return acc
end function

function rangedbl(lo as long=0, hi as long=1) as long
return (rand()/rand_max)*(hi+1-lo) -.5 + lo
End Function

function rangemod(lo as long=0, hi as long=1) as long
return (rand() mod (hi-lo +1))+ lo
End Function

srand(timer)
'================   TEST ============
dim as long lim=5000000
dim as long lo=0
dim as long hi=1

redim as double a(1 to lim),b(1 to lim)
for n as long=1 to lim
a(n)=(rangemod(lo,hi))
b(n)=(rangedbl(lo,hi))
next

for n as long=lbound(a) to 20
print n, a(n),b(n)
next
print "..."

print

dim as range R
dim as double m
var c=chisquared(a(),R,m)
print "Mod"
print "Mean ";m;"  Standard dev. ";R.sd;" Min ";R._min;" Max ";R._max
print "chi squared distance ";c
print
print
print "double"
c=chisquared(b(),R,m)
print "Mean ";m;"  Standard dev. ";R.sd;" Min ";R._min;" Max ";R._max
print "chi squared distance ";c
print
sleep