Re: Implementing the KISS4691 RNG

Re: Implementing the KISS4691 RNG

[robin]

"geo" <gmarsaglia@gmail.com> wrote in message news:93b6ca11-f4ae-4598-98f8-9a10be69ca3a@j18g2000yqd.googlegroups.com...
|
| For those interested in the KISS4691 RNG, I would
| like to recommend a final form for the MWC
| (Multiply-With-Carry) component, one that provides
| a pattern for applications with signed integers,
| takes care of the nuisance cases where (x<<13)+c
| overflows, and, with different choice of multiplier
| and prime but with the same mathematics, permits
| going through a full cycle to verify the period.
|
| First, a simple version for confirming the period:
| Here we still have b=2^32, but multiplier a=5 and
| prime p=ab-1, for which the order of b mod p is
|    (p-1)/2 = 5*2^31-1 = 10737418239
|
| Note that, as with the full MWC() function below, we
| deliberately seek multipliers 'a' of the form 2^i+1
| to ensure the ability to carry out the MWC operations
| in mod 2^32 arithmetic.  This required a search for primes
| of the form (2^i+1)*b^r-1, and I as able to find
| the prime (2^13+1)*b^4691-1 for the KISS4691 version.
|
| But the simpler version with p=(2^2+1)b-1 permits going
| through a full cycle (using only 32-bit operations)
| in less that a minute.   Here we have  a function
|
| f([x,c])=[(5x+c) mod 2^32,trunc((5x+c)/2^32)]
|
| with inverse
|
| f^{-1}([x,c])=[trunc((c*2^32+x)/5),(c*2^32+x) mod 5]
|
| Initial pairs  [x,c]=[0,0] or [2^32-1,4] will make the
| sequence f([x,c]),f(f([x,c])),f(f(f([x,c]))),... have
| period 1. All others, with 0<=x<2^32, 0<=c<5, will have
| period (p-1)/2=10737418239 and the (hexed) x's will form,
| in reverse order, the hex expansion of j/p for some 0<j<p.
|
| Running this C program:
| -----------------------------------------------
| #include <stdio.h>
| int main()
| { unsigned long x=123456789,c=3,t;
|  unsigned long long k;
|  for(k=0;k<10737418239LL;k++)
|  {t=(x<<2)+c;
|   if(t<c) { c=(x>>30)+1; x=t+x;}
|   else {c=(x>>30); x=t+x; c+=(x<t);}
|  }
| printf("%u,%u\n",x,c);
| }
| -----------------------------------------------
| will go through a full cycle and, after some
| 35 seconds, reproduce the initial x,c pair:
| 123456789,3
| Try it.
|
| For the KISS4691 version, the MWC multiplier is a=2^13+1,
| and I had a mental image of c taking 13 bits with the
| rightmost 13 bits of (x<<13) all 0's, so that (x<<13)+c
| could not overflow---indeed, that was the reason I sought
| multipliers of the form a=2^i+1.  But I overlooked the
| rare case where c could occupy 14 bits: c=8193, and cause
| an overflow when the rightmost 19 bits of x were all 1's.
| Astute readers pointed that out.  The above version, in
| which we first test for overflow after t=(x<<2)+c, then,
| if necessary, test with t=t+x, covers all cases and has the
| advantage of permitting similar structures for programs
| using signed integers---if we use a clever device advocated
| by Glen Herrmannsfeldt when providing a Fortran version of
| KISS4691 (comp.lang.fortran Aug 18, 2010).
|
| With t=x+y for integers, the C version of the overflow for
| signed integers, (t<x), may be replaced by (t'<x') for signed
| integers, where the prime means: flip the sign bit.  Thus,
| in C, with m=(1<<31), the overflow is ((t^m)<(x^m)), while
| for Fortran, with m=ishft(1,31), the logic is
|  ieor(t,m) .lt. ieor(x,m).
|
|
| Thus I suggest using this listing for the MWC function
| in KISS4691:
| ------------------------------------------------
|
| unsigned long MWC(void)
| {unsigned long t,x; static c=0,j=4691;
|   j=(j<4690)? j+1:0;
|   x=Q[j]; t=(x<<13)+c;
|     if(t<c) {c=(x>>19)+1; t=t+x;}
|     else {t=t+x; c=(x>>19)+(t<x);}
|   return (Q[j]=t);
| }
|
| ------------------------------------------------
|
| and recommend that it be the pattern for implementations in
| PL/I, asm, Fortran or others.  For signed integers, flip the
| sign bits so that t<x becomes (t^m)<(x^m) or equivalents.
| (C versions using signed integers should also avoid the
| problem of right shifts, replacing (x>>19) by
| ((unsigned)x>>19) or ((unsigned)xs>>17) in the XS component.)
|
|
| In case you do not have, or don't want to fetch, the original,
| here is the entire listing, with the recommended form for MWC():
|
| ------------------------------------------------------------
| static unsigned long xs=521288629,xcng=362436069,Q[4691];
|
| unsigned long MWC(void)
| {unsigned long t,x; static c=0,j=4691;
|  j=(j<4690)? j+1:0;
|  x=Q[j]; t=(x<<13)+c;
|    if(t<c) {c=(x>>19)+1; t=t+x;}
|    else {t=t+x; c=(x>>19)+(t<x);}
|  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 )
|
| #include <stdio.h>
| int main()
| {unsigned long i,x;
| for(i=0;i<4691;i++)  Q[i]=CNG+XS;
| printf("Does  MWC result=3740121002 ?\n");
| for(i=0;i<1000000000;i++) x=MWC();
|   printf("%27u\n",x);
|   printf("Does KISS result=2224631993 ?\n");
| for(i=0;i<1000000000;i++) x=KISS;
|   printf("%27u\n",x);
| }
| ------------------------------------------------------------
|
| Unfortunately, even at the rate of 238 million per second
| for MWC(), it would take over 10^40407 years to verify the
| period by running through a complete loop.
|
| George Marsaglia


Here is the PL/I version, using 32-bit unsigned arithmetic :-

(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);
   declare (t,x,i) fixed binary (32) unsigned;
   declare (c initial (0), j initial (4691) ) fixed binary (32) unsigned static;

   if j < hbound(Q,1) then j = j + 1; else j = 0;
   x = Q(j); t = isll(x,13) + c;
   if t<c then
      do; c = isrl(x, 19) + 1; t = t + x; end;
   else
      do; t = t + x; c = isrl(x, 19) + (t<x); end;
   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 );
end KISS;

   declare (i,x) fixed binary (32) unsigned;

   Q = CNG + XXS; /* Initialize: */
   do i = 1 to 1000000000; x=MWC(); end;
   put skip edit (" MWC result=3740121002",x) (a, skip, f(22));
   do i = 1 to 1000000000; x=KISS; end;
   put skip edit ("KISS result=2224631993",x) (a, skip, f(22));

end RNG; 

Re: Implementing the KISS4691 RNG

[robin]

"geo" <gmarsaglia@gmail.com> wrote in message news:93b6ca11-f4ae-4598-98f8-9a10be69ca3a@j18g2000yqd.googlegroups.com...
|
| For signed integers, flip the
| sign bits so that t<x becomes (t^m)<(x^m) or equivalents.
| (C versions using signed integers should also avoid the
| problem of right shifts, replacing (x>>19) by
| ((unsigned)x>>19) or ((unsigned)xs>>17) in the XS component.)

In the case of signed integers, both  t<c  AND  t<x
must have their sign bits flipped, not just  t<x.

Thus, the signed version in PL/I is:

(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 m fixed binary (31) static initial ((isll(1,31)));

   if j < hbound(Q,1) then j = j + 1; else j = 0;
   x = Q(j); t = isll(x,13) + c;
   if ieor(t,m) < ieor(c,m) then
      do; c = isrl(x, 19) + 1; t = t + x; end;
   else
      do; t = t + x; c = isrl(x, 19) + (ieor(t,m) < ieor(x,m)); end;
   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 (10);

   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;