Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: KISS4691, a potentially top-ranked RNG.
Replies: 108   Last Post: Dec 6, 2010 1:29 AM

 Messages: [ Previous | Next ]
 geo Posts: 38 Registered: 6/19/08
Re: KISS4691, a potentially top-ranked RNG.
Posted: Jul 25, 2010 9:53 AM

On Jul 24, 11:34 pm, Gib Bogle <g.bo...@auckland.no.spam.ac.nz> wrote:
> geo wrote:
>

> > Languages requiring signed integers should use equivalents
> > of the same operations, except that the C statement:
> >        c=(t<x)+(x>>19);
> > can be replaced by that language's version of
> >   if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
> >            else c=1-sign(t)+(x>>19)
> > */

>
> Hi George,
>
> I translated this into Fortran, and found that I get different results than with
> C.  I've tracked the difference into MWC().  The following Fortran code, with my
> laborious comparison of two signed integers treating them as unsigned, gives
> correct results.  If I comment out the line
> c = tLTx + SHIFTR(x,19)
> and uncomment the following lines that implement your suggestion above to
> compute c, I get different results.
>
> integer function MWC()
> integer :: t, x, i
> integer, save :: c = 0, j = 4691
> integer :: tLTx
>
> if (j < 4690) then
>      j = j + 1
> else
>      j = 0
> endif
> x = Q(j)
> t = SHIFTL(x,13) + c + x
> if ((t >= 0 .and. x >= 0) .or. (t < 0 .and. x < 0)) then
>      if (t < x) then
>          tLTx = 1
>      else
>          tLTx = 0
>      endif
> elseif (x < 0) then
>      tLTx = 1
> elseif (t < 0) then
>      tLTx = 0
> endif
>
> c = tLTx + SHIFTR(x,19)
>
> !if (sign(1,SHIFTL(x,13)+c) == sign(1,x)) then
> !    c = sign(1,x) + SHIFTR(x,19)
> !else
> !    c = 1 - sign(1,t) + SHIFTR(x,19)
> !endif
> Q(j) = t
> MWC = t
> end function
>
> Is it the case that although your suggested workaround gives different results
> from the C code in this case, it is still equivalent as a RNG?
>
> Cheers
> Gib

Thanks very much for the Fortran version.
I made a mistake in the comment on versions
for signed integers. This:

Languages requiring signed integers should use equivalents
of the same operations, except that the C statement:
c=(t<x)+(x>>19);
can be replaced by that language's version of
if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
else c=1-sign(t)+(x>>19)

should have been:

Languages requiring signed integers should use equivalents
of the same operations, except that the C statement:
c=(t<x)+(x>>19);
can be replaced by that language's version of
if sign(x<<13+c)=sign(x) then c=sign(x)+(x>>19)
else c=1-sign(x<<13+c)+(x>>19)

Sorry for that error.

I still like inline functions in Fortan,
so would tend to define
isign(x)=ishft(x,-31)
and
m=ishft(x,13)+c
if(isign(m).eq.isign(x)) then c=isign(x)+ishft(x,-19)
else c=1-isign(m)+ishft(x,-19)
and
Q[j]=m+x

If calculating the carry c of the MWC operation fails
to fix that extra increment properly, then rather than
a systematic expansion, in reverse order, 32 bits at a time,
of the binary representation of j/(1893*2^150112-1) for some
j determined by the random seeds, we will be jumping around
in that expansion, and we can't be sure that the period will
still be the order of b=2^32 for the prime p=1893*b^4196-1.

gm

Date Subject Author
7/24/10 geo
7/24/10 Geoff
7/24/10 jacob navia
7/27/10 robin
7/30/10 Uno
7/24/10 Gene
7/24/10 Gib Bogle
7/25/10 geo
7/25/10 Dick Hendrickson
7/25/10 Gib Bogle
7/26/10 robin
7/26/10 Geoff
7/26/10 Gib Bogle
7/26/10 robin
7/27/10 Geoff
7/27/10 nmm1@cam.ac.uk
7/26/10 Gib Bogle
7/29/10 Uno
7/29/10 Gib Bogle
7/29/10 Uno
7/30/10 Gib Bogle
7/30/10 glen herrmannsfeldt
7/30/10 Gib Bogle
7/30/10 robin
7/30/10 robin
7/30/10 Gib Bogle
7/25/10 robin
7/26/10 Harald Anlauf
7/26/10 nmm1@cam.ac.uk
7/26/10 glen herrmannsfeldt
7/29/10 Uno
7/29/10 glen herrmannsfeldt
7/29/10 Ron Shepard
7/29/10 glen herrmannsfeldt
7/30/10 Uno
7/30/10 nmm1@cam.ac.uk
7/30/10 Ron Shepard
7/31/10 Ilmari Karonen
7/31/10 Uno
7/31/10 sturlamolden
7/30/10 David Bernier
7/30/10 orz
7/30/10 glen herrmannsfeldt
8/1/10 Ron Shepard
8/1/10 glen herrmannsfeldt
7/30/10 orz
7/30/10 orz
7/30/10 Gib Bogle
7/30/10 orz
7/30/10 Gib Bogle
7/31/10 orz
7/31/10 Gib Bogle
8/4/10 orz
8/5/10 Uno
8/10/10 Gib Bogle
8/10/10 orz
8/10/10 Gib Bogle
8/10/10 Dann Corbit
8/11/10 Richard Maine
8/28/10 robin
8/12/10 robin
8/17/10 Gib Bogle
8/17/10 glen herrmannsfeldt
8/17/10 gnasher729
8/18/10 geo
8/18/10 gnasher729
8/19/10 gnasher729
8/19/10 geo
8/19/10 glen herrmannsfeldt
8/19/10 orz
8/19/10 orz
8/21/10 Guest
8/21/10 Guest
8/17/10 Gib Bogle
8/18/10 glen herrmannsfeldt
8/18/10 Gib Bogle
8/18/10 baf
8/19/10 Gib Bogle
8/19/10 baf
8/19/10 Gib Bogle
8/18/10 Gib Bogle
8/28/10 robin
8/28/10 robin
8/10/10 orz
8/11/10 glen herrmannsfeldt
12/6/10 Michael Press
7/26/10 Harald Anlauf
7/26/10 nmm1@cam.ac.uk
7/27/10 robin
7/30/10 Uno
8/3/10 robin
8/1/10 sturlamolden
8/2/10 Georg Beyerle
8/2/10 Uno
8/2/10 sturlamolden
8/3/10 Dann Corbit
8/3/10 sturlamolden
8/4/10 orz