Velocity Reviews > KISS4691, a potentially top-ranked RNG.

# KISS4691, a potentially top-ranked RNG.

geo
Guest
Posts: n/a

 07-24-2010
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

/*
The KISS4691 RNG, a Keep-It-Simple-Stupid combination of
MWC (Multiply-With-Carry), Congruential and Xorshift sequences.
Expressed as a simple MWC() function and C macros
#define CNG ( xcng=69069*xcng+123 )
#define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
#define KISS ( MWC()+CNG+XS )
but easily expressed in other languages, with a slight
complication for signed integers.

With its immense period, >10^45000, and speed: MWC()s at
around 238 million/sec or full KISSes at around 138 million,
there are few RNGs that do as well as this one
on tests of randomness and are comparable in even one
of the categories: period, speed, simplicity---not to
mention comparable in all of them.

The MWC4691 sequence x[n]=8193*x[n-4691]+carry mod b=2^32
is based on p=8193*b^4691-1, period ~ 2^150124 or 10^45192
For the MWC (Multiply-With-Carry) process, given a current
x and c, the new x is (8193*x+c) mod b,
the new c is the integer part of (8193*x+c)/b.

The routine MWC() produces, in reverse order, the base-b=2^32
expansion of some j/p with 0<j<p=8193*2^150112-1 with j
determined by the 64 bits of seeds xcng and xs, or more
generally, by 150112 random bits in the Q[] array.
*/

static unsigned long xs=521288629,xcng=362436069,Q[4691];

unsigned long MWC(void) /*takes about 4.2 nanosecs or 238 million/
second*/
{unsigned long t,x,i; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j];
t=(x<<13)+c+x; c=(t<x)+(x>>19);
return (Q[j]=t);
}

#define CNG ( xcng=69069*xcng+123 )
#define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
#define KISS ( MWC()+CNG+XS ) /*138 million/sec*/

#include <stdio.h>
int main()
{unsigned long i,x;
for(i=0;i<4691;i++) Q[i]=CNG+XS;
for(i=0;i<1000000000;i++) x=MWC();
printf(" MWC result=3740121002 ?\n%22u\n",x);
for(i=0;i<1000000000;i++) x=KISS;
printf("KISS result=2224631993 ?\n%22u\n",x);
}

/*
This RNG uses two seeds to fill the Q[4691] array by
means of CNG+XS mod 2^32. Users requiring more than two
seeds will need to randomly seed Q[] in main().
By itself, the MWC() routine passes all tests in the
Diehard Battery of Tests, but it is probably a good
idea to include it in the KISS combination.

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)
*/

jacob navia
Guest
Posts: n/a

 07-24-2010
geo a écrit :
> 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
>
>
> /*
> The KISS4691 RNG, a Keep-It-Simple-Stupid combination of
> MWC (Multiply-With-Carry), Congruential and Xorshift sequences.
> Expressed as a simple MWC() function and C macros
> #define CNG ( xcng=69069*xcng+123 )
> #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
> #define KISS ( MWC()+CNG+XS )
> but easily expressed in other languages, with a slight
> complication for signed integers.
>
> With its immense period, >10^45000, and speed: MWC()s at
> around 238 million/sec or full KISSes at around 138 million,
> there are few RNGs that do as well as this one
> on tests of randomness and are comparable in even one
> of the categories: period, speed, simplicity---not to
> mention comparable in all of them.
>
> The MWC4691 sequence x[n]=8193*x[n-4691]+carry mod b=2^32
> is based on p=8193*b^4691-1, period ~ 2^150124 or 10^45192
> For the MWC (Multiply-With-Carry) process, given a current
> x and c, the new x is (8193*x+c) mod b,
> the new c is the integer part of (8193*x+c)/b.
>
> The routine MWC() produces, in reverse order, the base-b=2^32
> expansion of some j/p with 0<j<p=8193*2^150112-1 with j
> determined by the 64 bits of seeds xcng and xs, or more
> generally, by 150112 random bits in the Q[] array.
> */
>
> static unsigned long xs=521288629,xcng=362436069,Q[4691];
>
> unsigned long MWC(void) /*takes about 4.2 nanosecs or 238 million/
> second*/
> {unsigned long t,x,i; static c=0,j=4691;
> j=(j<4690)? j+1:0;
> x=Q[j];
> t=(x<<13)+c+x; c=(t<x)+(x>>19);
> return (Q[j]=t);
> }
>
> #define CNG ( xcng=69069*xcng+123 )
> #define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
> #define KISS ( MWC()+CNG+XS ) /*138 million/sec*/
>
> #include <stdio.h>
> int main()
> {unsigned long i,x;
> for(i=0;i<4691;i++) Q[i]=CNG+XS;
> for(i=0;i<1000000000;i++) x=MWC();
> printf(" MWC result=3740121002 ?\n%22u\n",x);
> for(i=0;i<1000000000;i++) x=KISS;
> printf("KISS result=2224631993 ?\n%22u\n",x);
> }
>
> /*
> This RNG uses two seeds to fill the Q[4691] array by
> means of CNG+XS mod 2^32. Users requiring more than two
> seeds will need to randomly seed Q[] in main().
> By itself, the MWC() routine passes all tests in the
> Diehard Battery of Tests, but it is probably a good
> idea to include it in the KISS combination.
>
> 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)
> */

This doesn't work with systems that have unsigned long as a 64 bit quantity.

I obtain:

MWC result=3740121002 ?
4169348530
KISS result=2224631993 ?
1421918629

Compiling with 32 bit machine yields:
MWC result=3740121002 ?
3740121002
KISS result=2224631993 ?
2224631993

Gib Bogle
Guest
Posts: n/a

 07-25-2010
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

geo
Guest
Posts: n/a

 07-25-2010
On Jul 24, 11:34*pm, Gib Bogle <(E-Mail Removed)> 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

Dick Hendrickson
Guest
Posts: n/a

 07-25-2010
On 7/25/10 8:53 AM, geo wrote:
[snip]
>
> I still like inline functions in Fortan,
> so would tend to define
> isign(x)=ishft(x,-31)

I don't know anything about random numbers, but "sign" is an intrinsic
function in Fortran and compilers should generate near perfect code
inline for things like sign(1,x). The generic sign function was added
to Fortran 90. In FORTRAN 77, you'd need to use the specific ISIGN(1,X)
function.

Dick Hendrickson

[snip]
>
> gm

Gib Bogle
Guest
Posts: n/a

 07-26-2010
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

robin
Guest
Posts: n/a

 07-26-2010
"geo" <(E-Mail Removed)> wrote in message news:(E-Mail Removed)...
|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;

robin
Guest
Posts: n/a

 07-26-2010
"Gib Bogle" <(E-Mail Removed)> wrote in message news:i2ij74\$kd6\$(E-Mail Removed)...
| 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 ...

Harald Anlauf
Guest
Posts: n/a

 07-26-2010
On Jul 24, 3:02*pm, geo <(E-Mail Removed)> wrote:
> This RNG uses two seeds to fill the Q[4691] array by
> means of CNG+XS mod 2^32. Users requiring more than two
> seeds will need to randomly seed Q[] in main().
> By itself, the MWC() routine passes all tests in the
> Diehard Battery of Tests, but it is probably a good
> idea to include it in the KISS combination.

Does this mean that using different seeds will lead to
streams that are always statistically independent
(as long as one does not exhaust the RNG's period)?
Or are there restrictions on the possible combinations
of seeds?

I am currently using D.E.Knuth's generator from TAOCP,
which IIRC allows for 2^30-2 independent streams, and
which are asserted to be independent, but being able to
extend the "limit" would be nice.

Harald

Gib Bogle
Guest
Posts: n/a

 07-26-2010
robin wrote:
> "Gib Bogle" <(E-Mail Removed)> wrote in message news:i2ij74\$kd6\$(E-Mail Removed)...
> | 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 hope so. Maybe I didn't express myself clearly enough. I'll try again.
Using my implementation of George's corrected code, I get results from the
Fortran code that differ from those generated by his C code.