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

General FreeBASIC programming questions.
deltarho
Posts: 1906
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
Posts: 1906
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: 3216
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

deltarho 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 maximumPrint nllmax Mod 19Print pllmax Mod 19Print : Print "done ... ";Sleep`
result wrote:-1
17
Added: the neg. max. referres to binary = set bit's (value -1)
dodicat
Posts: 5913
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 zDim answer As String   'THE ANSWER STRING  '_______ 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) Thenpart1=Mid\$(s,1,position-1)part2=Mid\$(s,position)s=part1+char+part2End if#endmacroinsert(answer,".",decpos)  answer=thepoint+zeros+answerIf dpflag="raw" Then    answer=Mid(answer,1,decimal_places)    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 DIVISIONDim As integer lennint,lenddec,lend,lenn,difflensplit(number,".",var1,var2)lennint=Len(var1)split(divisor,".",var1,var2)lenddec=Len(var2)If Instr(number,".") Then     Removepoint(number)    number=var1+var2    End ifIf Instr(divisor,".") Then     Removepoint(divisor)    divisor=var1+var2    End ifDim As integer numzerosnumzeros=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-lennDim decpos As integer=lenddec+lennint-lend+2-numzeros 'THE POSITION INDICATORDim _sgn As Byte=-Sgn(decpos)If _sgn=0 Then _sgn=1Dim As String thepoint=String(_sgn,".") 'DECIMAL AT START (IF)Dim As String zeros=String(-decpos+1,"0")'ZEROS AT START (IF) e.g. .0009if dpflag<>"mod" thenIf Len(zeros) =0 Then dpflag="s"end ifDim As integer runlengthIf Len(zeros) Then      runlength=decimal_places     answer=String(Len(zeros)+runlength+10,"0")    If dpflag="raw" Then         runlength=1        answer=String(Len(zeros)+runlength+10,"0")        If decimal_places>Len(zeros) Then            runlength=runlength+(decimal_places-Len(zeros))            answer=String(Len(zeros)+runlength+10,"0")            End If            End IfElsedecimal_places=decimal_places+decposrunlength=decimal_placesanswer=String(Len(zeros)+runlength+10,"0")End if'___________DECIMAL POSITION DETERMINED  _____________'SET UP THE VARIABLES AND START UP CONDITIONSnumber=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 'notIf Lens>lenf Thentemp= "done"Exit DoEnd ifIf divisor>topstring Then temp= "done"Exit DoEnd ifEnd 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       answer[runcount-1]=count+47       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       answer=Rtrim(answer,"0")       answer=Rtrim(answer,".")       answer=Ltrim(answer,"0")       If answer="" Then Return "0"   Return sign+answerEnd Function#define mod_(a,b) _divide(a,b,,"mod")#define div_(a,b) _divide(a,b,-2)dim as ulongint R=2^32-1   + 1000000dim as ulongint n=234print R Mod n' returns zero if R >= 2^32-1.dim as ulongint u=-1for 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))    printnextsleep `

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

Code: Select all

` Dim As Longint lo = 20, hi = 147, R, ansDim As Double x   R = 7547257594037927936ans =  R Mod (hi - lo + 1)  + 20 ' dodicatPrint ansx = R/2^64*(hi - lo + 1) - 0.5 + 20ans = R/2^64*(hi - lo + 1) - 0.5 + 20  ' dafhiPrint ans, xprintprintprint "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   nextSleep `

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: 3216
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
Posts: 1906
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: 3216
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
Posts: 1906
Joined: Jan 02, 2017 0:34
Location: UK

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

No further comments until admin get back to me.
MrSwiss
Posts: 3216
Joined: Jun 02, 2013 9:27
Location: Switzerland

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

deltarho wrote:No further comments until admin get back to me.

Fine with me. No problem.
dodicat
Posts: 5913
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 functiondim as ulongint lo=0,hi=10,y,y2,y3dim as ulongint R=7547257594037927936y=(R/2^64)*(hi-lo+1)-.5 +loprint y;" by dafhi"y= R Mod (hi - lo + 1)  + loprint y; " by mod"y= map(0,2^64-1,R,lo,hi)print y;" by mapping"printprint "press a key"sleepprintr=0print "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
Posts: 1906
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,y3Dim As Ulongint R=7547257594037927936y=(R/2^64)*(hi-lo+1)-.5 +loPrint y;" by dafhi"y= R Mod (hi - lo + 1)  + loPrint y; " by Mod"PrintPrint "press a key"SleepPrint 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=hiSleep`

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 dafhi58 by Mod press a key R                            dafhi               Mod0                            50                  50368934881474191032           51                  65737869762948382064           52                  801106804644422573096          53                  951475739525896764128          54                  591844674407370955160          55                  742213609288845146192          56                  892582544170319337224          57                  532951479051793528256          58                  683320413933267719288          59                  833689348814741910320          60                  984058283696216101352          61                  624427218577690292384          62                  774796153459164483416          63                  925165088340638674448          64                  565534023222112865480          65                  715902958103587056512          66                  866271892985061247544          67                  506640827866535438576          68                  657009762748009629608          69                  807378697629483820640          70                  957747632510958011672          71                  598116567392432202704          72                  748485502273906393736          73                  898854437155380584768          74                  539223372036854775800          75                  689592306918328966832          76                  839961241799803157864          77                  9810330176681277348896         78                  6210699111562751539928         79                  7711068046444225730960         80                  9211436981325699921992         81                  5611805916207174113024         82                  7112174851088648304056         83                  8612543785970122495088         84                  5012912720851596686120         85                  6513281655733070877152         86                  8013650590614545068184         87                  9514019525496019259216         88                  5914388460377493450248         89                  7414757395258967641280         90                  8915126330140441832312         91                  5315495265021916023344         92                  6815864199903390214376         93                  8316233134784864405408         94                  9816602069666338596440         95                  6216971004547812787472         96                  7717339939429286978504         97                  9217708874310761169536         98                  5618077809192235360568         99                  7118446744073709551600         100                 86`
dodicat
Posts: 5913
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
Posts: 1906
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: 1243
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 8End Functionfunction range(hi as single=1, lo as single=0)as ubyte  return lcg mod (hi+1-lo) + loEnd Functiondim as long bins(255)for i as long = 1 to 256 * 100  bins( lcg ) += 1Nextvar vlo = bins(0), vhi = vlofor i as long = 1 to 255  if vlo > bins(i) then vlo = bins(i)  if vhi < bins(i) then vhi = bins(i)  bins(i) = 0Next: 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: 5913
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,_miniend typefunction 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 mend functionfunction 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 accend functionfunction rangedbl(lo as long=0, hi as long=1) as long  return (rand()/rand_max)*(hi+1-lo) -.5 + loEnd Functionfunction rangemod(lo as long=0, hi as long=1) as long  return (rand() mod (hi-lo +1))+ loEnd Functionsrand(timer)'================   TEST ============dim as long lim=5000000dim as long lo=0dim as long hi=1redim 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))nextfor n as long=lbound(a) to 20    print n, a(n),b(n)nextprint "..."printdim as range Rdim as double mvar c=chisquared(a(),R,m)print "Mod"print "Mean ";m;"  Standard dev. ";R.sd;" Min ";R._min;" Max ";R._maxprint "chi squared distance ";cprintprintprint "double"c=chisquared(b(),R,m)print "Mean ";m;"  Standard dev. ";R.sd;" Min ";R._min;" Max ";R._maxprint "chi squared distance ";cprintsleep  `