G
geo
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=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
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=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