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 KeepItSimpleStupid combination of > MWC (MultiplyWithCarry), 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, simplicitynot to > mention comparable in all of them. > > The MWC4691 sequence x[n]=8193*x[n4691]+carry mod b=2^32 > is based on p=8193*b^46911, period ~ 2^150124 or 10^45192 > For the MWC (MultiplyWithCarry) 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 baseb=2^32 > expansion of some j/p with 0<j<p=8193*2^1501121 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=1sign(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 32bit 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:9ea9185a23e74265acde097c5c14ca14(a)y11g2000yqm.googlegroups.com... "geo" <gmarsaglia(a)gmail.com> wrote in message news:a82cebe3cdb948af8080bca935eeb9b1(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 46911; Q(i) = CNG+XXS; end; put skip list (q(0), q(4690)); put skip list ('initialized'); put skip; do i = 0 to 10000000001; x=MWC(); end; put skip edit (" MWC result=3740121002 ",x) (a, f(23)); do i = 0 to 10000000001; 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:a82cebe3cdb948af8080bca935eeb9b1(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 64bit machine (using 64bit 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 32bit machine, these extra bits are discarded, but in 64bit 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 32bit and 64bit 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 64bit machine (using 64bit 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 32bit machine, these extra bits are discarded, > but in 64bit 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 32bit and > 64bit 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 nonvoid 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) 