Prev: books on numerical programming in Ada
Next: Making measurements (Was: compiler settings in AdaGIDE)
From: Gene on 24 Jul 2010 20:49
On Jul 24, 9:02 am, geo <gmarsag...(a)gmail.com> wrote:
> 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)
> */

Here's an implementation in Ada, verified to produce the same answers
as the C code.

The Ada version of gcc compiles this to nearly the same code as
produced
for the C. Function call overhead adds about 40% to the run time.

What's a good way to take an arbitrary, single 32-bit seed value to
a complete state initialization? Accept that lots of people will
pick very small numbers.

-- kiss_random_numbers.ads
-- Specification for George Marsaglia's KISS random number generator.
package KISS_Random_Numbers is

type Result_Type is mod 2 ** 32;

type State_Type is private;

function Next_Random_Value(State : access State_Type) return
Result_Type;

private

type State_Index_Type is mod 4691;
type State_Vector_Type is array (State_Index_Type) of Result_Type;
Initial_Q : State_Vector_Type; -- set in package body

type Substate_Type is
record
Xs : Result_Type := 521288629;
Xcng : Result_Type := 362436069;
end record;

type State_Type is
record
Sub : aliased Substate_Type;
C : Result_Type := 0;
Q : State_Vector_Type := Initial_Q;
J : State_Index_Type := State_Index_Type'Last;
end record;

end KISS_Random_Numbers;

-- kiss_random_numbers.adb
-- Implementation of George Marsaglia's KISS random number generator.
package body KISS_Random_Numbers is

function Mwc (State : access State_Type) return Result_Type is
T, X : Result_Type;
begin
State.J := State.J + 1;
X := State.Q(State.J);
T := X * 2**13 + State.C + X;
if T < X then
State.C := X / 2**19 + 1;
else
State.C := X / 2**19;
end if;
State.Q(State.J) := T;
return T;
end Mwc;
pragma Inline(Mwc);

function Cng (State : access Substate_Type) return Result_Type is
begin
State.Xcng := 69069 * State.Xcng + 123;
return State.Xcng;
end Cng;
pragma Inline(Cng);

function Xs (State : access Substate_Type) return Result_Type is
Xs : Result_Type;
begin
Xs := State.Xs;
Xs := Xs xor (Xs * 2**13);
Xs := Xs xor (Xs / 2**17);
Xs := Xs xor (Xs * 2**5);
State.Xs := Xs;
return Xs;
end Xs;
pragma Inline(Xs);

function Next_Random_Value(State : access State_Type) return
Result_Type is
begin
return Mwc(State) + Cng(State.Sub'Access) +
Xs(State.Sub'Access);
end Next_Random_Value;

begin
declare
S : aliased Substate_Type;
begin
for I in Initial_Q'Range loop
Initial_Q(I) := Cng(S'access) + Xs(S'Access);
end loop;
end;
end KISS_Random_Numbers;

-- test_kiss.adb
-- Simple test of George Marsaglia's KISS random number generator.
with Ada.Text_IO; use Ada.Text_IO;
with KISS_Random_Numbers;
use KISS_Random_Numbers;

procedure Test_Kiss is
X : Result_Type;
S : aliased State_Type;
begin
for I in 1 .. 1000000000 loop
X := Next_Random_Value(S'Access);
end loop;
Put(Result_Type'Image(X));
New_Line;
end;
From: robin on 25 Jul 2010 22:50
"Gene" <gene.ressler(a)gmail.com> wrote in message
news:9ea9185a-23e7-4265-acde-097c5c14ca14(a)y11g2000yqm.googlegroups.com...
"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 27 Jul 2010 01:46
"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.

I have already posted a PL/I version using unsigned arithmetic.

Here is another version, this time using signed arithmetic :--

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

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

MWC: procedure () returns (fixed binary (31));
declare (t,x,i) fixed binary (31);
declare (c initial (0), j initial (4691) ) fixed binary (31) static;
declare (t1, t2, t3) fixed binary (31);

if j < hbound(Q,1) then j = j + 1; else j = 0;
x = Q(j);
t = isll(x,13)+c+x;
t1 = iand(x, 3) - iand(t, 3);
t2 = isrl(x, 2) - isrl(t, 2);
if t2 = 0 then t2 = t1;
if t2 > 0 then t3 = 1; else t3 = 0;
c = t3 + isrl(x, 19);
Q(j)=t;
return (t);
end MWC;

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

XXS: procedure returns (fixed binary (31));
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 (31));
return ( MWC()+CNG+XXS );
end KISS;

declare (i,x) fixed binary (31);
declare y fixed decimal (11);

Q = CNG+XXS; /* Initialize. */
do i = 1 to 1000000000; x=MWC(); end;
put skip edit (" Expected MWC result = 3740121002", 'computed =', x)
(a, skip, x(12), a, f(11));
y = iand(x, 2147483647);
if x < 0 then y = y + 2147483648;
put skip edit (y) (x(11), f(22)); put skip;
do i = 1 to 1000000000; x=KISS; end;
put skip edit ("Expected KISS result = 2224631993", 'computed =', x)
(a, skip, x(12), a, f(11));
y = iand(x, 2147483647);
if x < 0 then y = y + 2147483648;
put skip edit (y) (x(11), f(22));

end RNG;


From: robin on 27 Jul 2010 06:19
"jacob navia" <jacob(a)spamsink.net> wrote in message news:i2fir2$op4$1(a)speranza.aioe.org...

| 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

For a 64-bit machine (using 64-bit integer arithmetic),
you'd need to truncate each result to 32 bits. That not
only applies to the multiplication, it also applies to addition, etc.
On a 32-bit machine, these extra bits are discarded,
but in 64-bit arithmetic, they are retained,
and unless they are similarly discarded,
you won't get the same results.
I suggest using IAND(k, 2*2147483647+1)
for the truncation.

With such modifications in the program,
it should then produce the same results on both 32-bit and
64-bit machines.

P.S. the product 2*2... is best obtained using ISHFT.

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



From: Uno on 30 Jul 2010 04:33
robin wrote:
> "jacob navia" <jacob(a)spamsink.net> wrote in message news:i2fir2$op4$1(a)speranza.aioe.org...
>
> | 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
>
> For a 64-bit machine (using 64-bit integer arithmetic),
> you'd need to truncate each result to 32 bits. That not
> only applies to the multiplication, it also applies to addition, etc.
> On a 32-bit machine, these extra bits are discarded,
> but in 64-bit arithmetic, they are retained,
> and unless they are similarly discarded,
> you won't get the same results.
> I suggest using IAND(k, 2*2147483647+1)
> for the truncation.
>
> With such modifications in the program,
> it should then produce the same results on both 32-bit and
> 64-bit machines.
>
> P.S. the product 2*2... is best obtained using ISHFT.
>
> | Compiling with 32 bit machine yields:
> | MWC result=3740121002 ?
> | 3740121002
> | KISS result=2224631993 ?
> | 2224631993
>
>
>

First of all, I think we're talking to the actual George Marsaglia here.
Thank you so much for posting. You may have displaced Terence as our
senior member.

Are you creating bigger numbers just to accomodate your age?:-)

$ gcc -Wall -Wextra geo1.c -o out
geo1.c: In function �MWC�:
geo1.c:5: warning: type defaults to �int� in declaration of �c�
geo1.c:5: warning: type defaults to �int� in declaration of �j�
geo1.c:5: warning: unused variable �i�
geo1.c: In function �main�:
geo1.c:21: warning: format �%22u� expects type �unsigned int�, but
argument 2 has type �long unsigned int�
geo1.c:23: warning: format �%22u� expects type �unsigned int�, but
argument 2 has type �long unsigned int�
geo1.c:24: warning: control reaches end of non-void function

$ ./out
MWC result=3740121002 ?
3740121002
KISS result=2224631993 ?
2224631993
$ cat geo1.c
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);
}

// gcc -Wall -Wextra geo1.c -o out
$

So, what is all this? In particular, is there something special about
the value of 3.7 billion?
--
Uno
 |  Next  |  Last
Pages: 1 2 3 4
Prev: books on numerical programming in Ada
Next: Making measurements (Was: compiler settings in AdaGIDE)