Re: RNGs with periods exceeding 10^(40million).

["robin" <robin51@dodo.mapson.com.au>]

"geo" <gmarsaglia@gmail.com> wrote in message news:603ebe15-a32f-4fbb-ba44-6c73f7919a33@t35g2000yqj.googlegroups.com...
|
| I offer here a MWC (Multiply-With-Carry) RNG whose
| period I cannot determine---nor is it likely that
| anyone else can---but that unknown period is almost
| certainly greater than 10^40000000, i.e, 10^(40million).

Here is the 64-bit version using PL/I :--

/* 64-bit version in PL/I using unsigned arithmetic  18/1/2011 */
(nosize):
rng64: procedure options (main);
   declare (Q(0:2097151), carry initial (0)) unsigned fixed binary (64)
                                             static;

B64MWC: procedure () returns (unsigned fixed binary (64));
   declare (t,x) unsigned fixed binary (64);
   declare j fixed binary (31) static initial (2097151);
   j=iand(j+1, 2097151);
   x=Q(j); t=isll(x, 28)+carry;
   carry=isrl(x, 36)-(t<x);
   Q(j)=t-x;
   return (Q(j));
end B64MWC;

KISS: procedure() returns (unsigned fixed binary (64));
   declare t unsigned fixed binary (64);
   t = B64MWC();
   cng=unsigned(6906969069, 64)*cng+unsigned(13579);
   xs = ieor(xs, isll(xs, 13));
   xs = ieor(xs, isrl(xs, 17));
   xs = ieor(xs, isll(xs, 43));
   return (t + cng + xs);
end KISS;

   declare (i,x,cng static initial (123456789987654321),
                 xs static initial (362436069362436069))
                    unsigned fixed binary (64);
   declare (start_time, finish_time) float (18);
   start_time = secs();

/* First seed Q with CNG+XS:    */
   do i=0 to 2097151;
      xs = ieor(xs, isll(xs, 13));
      xs = ieor(xs, isrl(xs, 17));
      xs = ieor(xs, isll(xs, 43));
      cng=unsigned(6906969069, 64)*cng+unsigned(13579);
      Q(i)=CNG+XS;
   end;

/* Then generate 10^9 B64MWC()s */
   do i=1 to 1000000000; x=B64MWC(); end;
   put skip list (x);
   if decimal(x) ^= 13596816608992115578 then
      put list ('is not the expected value 13596816608992115578');

/* followed by 10^9 KISSes:   */
   do i=1 to 1000000000; x=KISS(); end;
   put skip list (x);
   if decimal(x) ^= 5033346742750153761 then
      put list ('is not the expected value 5033346742750153761');

   Finish_time = secs();
   put skip edit (Finish_time - Start_time, ' secs') (F(10,2), A);

end rng64;