geo
Posts:
38
Registered:
6/19/08


Implementing the KISS4691 RNG
Posted:
Sep 5, 2010 5:04 PM


For those interested in the KISS4691 RNG, I would like to recommend a final form for the MWC (MultiplyWithCarry) 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=ab1, for which the order of b mod p is (p1)/2 = 5*2^311 = 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^r1, and I as able to find the prime (2^13+1)*b^46911 for the KISS4691 version.
But the simpler version with p=(2^2+1)b1 permits going through a full cycle (using only 32bit 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^321,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 (p1)/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 overflowindeed, 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 integersif 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

