Prev: solutions book
Next: real kind declaration
From: geo on 25 Jul 2010 09:53
On Jul 24, 11:34 pm, Gib Bogle <g.bo...(a)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
From: Gib Bogle on 25 Jul 2010 20:00
geo wrote:
> 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.

That produces different c values from those generated by the method based on the
value of (t<x), therefore it deviates from the C code. This is what I used:

m = shiftl(x,13) + c
if (sign(1,m) == sign(1,x)) then
c = sign(1,x) + shiftr(x,19)
else
c = 1 - sign(1,m) + shiftr(x,19)
endif
From: robin on 25 Jul 2010 22:38
"geo" <gmarsaglia(a)gmail.com> wrote in message news:a82cebe3-cdb9-48af-8080-bca935eeb9b1(a)l14g2000yql.googlegroups.com...
|I have been asked to recommend an RNG
| (Random Number Generator) that ranks
| at or near the top in all of the categories:
| performance on tests of randomness,
| length of period, simplicity and speed.
| The most important measure, of course, is
| performance on extensive tests of randomness, and for
| those that perform well, selection may well depend
| on those other measures.
|
| The following KISS version, perhaps call it KISS4691,
| seems to rank at the top in all of those categories.
| It is my latest, and perhaps my last, as, at age 86,
| I am slowing down.
|
| Compiling and running the following commented
| C listing should produce, within about four seconds,
| 10^9 calls to the principal component MWC(), then
| 10^9 calls to the KISS combination in another ~7 seconds.
|
| Try it; you may like it.
|
| George Marsaglia

Here's the PL/I version:

(NOSIZE, NOFOFL):
RNG: PROCEDURE OPTIONS (MAIN, REORDER);

declare (xs initial (521288629), xcng initial (362436069),
Q(0:4690) ) static fixed binary (32) unsigned;

MWC: procedure () returns (fixed binary (32) unsigned);
/*takes about 4.2 nanosecs or 238 million/second*/
declare (t,x,i) fixed binary (32) unsigned;
declare (c initial (0), j initial (4691) ) fixed binary (32) unsigned static;

if j < 4690 then j = j + 1; else j = 0;
x = Q(j);
t = isll(x,13)+c+x; c = (t<x)+isrl(x,19);
Q(j)=t;
return (t);
end MWC;

CNG: procedure returns (fixed binary (32) unsigned);
xcng=bin(69069)*xcng+bin(123);
return (xcng);
end CNG;

XXS: procedure returns (fixed binary (32) unsigned);
xs = ieor (xs, isll(xs, 13) );
xs = ieor (xs, isrl(xs, 17) );
xs = ieor (xs, isll(xs, 5) );
return (xs);
end XXS;

KISS: procedure returns (fixed binary (32) unsigned);
return ( MWC()+CNG+XXS ); /*138 million/sec*/
end KISS;

declare (i,x) fixed binary (32) unsigned;
/* Initialize: */
do i = 0 to 4691-1; Q(i) = CNG+XXS; end;
put skip list (q(0), q(4690));
put skip list ('initialized'); put skip;
do i = 0 to 1000000000-1; x=MWC(); end;
put skip edit (" MWC result=3740121002 ",x) (a, f(23));
do i = 0 to 1000000000-1; x=KISS; end;
put skip edit ("KISS result=2224631993 ",x) (a, f(23));

end RNG;


From: robin on 26 Jul 2010 09:32
"Gib Bogle" <g.bogle(a)auckland.no.spam.ac.nz> wrote in message news:i2ij74$kd6$1(a)speranza.aioe.org...
| geo wrote:
| > 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(x<<13+c)+(x>>19)
| >
| > Sorry for that error.
|
| That produces different c values from those generated by the method based on the
| value of (t<x), therefore it deviates from the C code. This is what I used:
|
| m = shiftl(x,13) + c
| if (sign(1,m) == sign(1,x)) then
| c = sign(1,x) + shiftr(x,19)
| else
| c = 1 - sign(1,m) + shiftr(x,19)
| endif

Maybe I missed something,
but isn't this exactly equivalent to what George wrote?
Just substitute x<<13+c for m in your last two assignments ...



From: Geoff on 26 Jul 2010 12:25
On Mon, 26 Jul 2010 23:32:27 +1000, "robin" <robin51(a)dodo.com.au>
wrote:

>"Gib Bogle" <g.bogle(a)auckland.no.spam.ac.nz> wrote in message news:i2ij74$kd6$1(a)speranza.aioe.org...
>| geo wrote:
>| > 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(x<<13+c)+(x>>19)
>| >
>| > Sorry for that error.
>|
>| That produces different c values from those generated by the method based on the
>| value of (t<x), therefore it deviates from the C code. This is what I used:
>|
>| m = shiftl(x,13) + c
>| if (sign(1,m) == sign(1,x)) then
>| c = sign(1,x) + shiftr(x,19)
>| else
>| c = 1 - sign(1,m) + shiftr(x,19)
>| endif
>
>Maybe I missed something,
>but isn't this exactly equivalent to what George wrote?
>Just substitute x<<13+c for m in your last two assignments ...
>
>

I am no FORTRAN hacker but I think there's a difference between
sign(x) and sign(1,x).
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4 5 6 7 8 9 10 11 12
Prev: solutions book
Next: real kind declaration