Ensuring the long period of the KISS4691 RNG.

G

geo

My posting for KISS4691 has become so mixed with
diverse comments over various newsgroups that I
thought it better to use a separate posting to
bring suggested changes to interested potential
users via sci.math, comp.lang.c, comp.lang.fortran.

At least two readers have pointed out special cases that
might arise from the code I suggested for the KISS4691.

christain.bau pointed out that (x<<13)+c
can exceed 2^32 when the rightmost 19 bits of x are 1's
and c is its maximum possible, a-1=8192=2^13.

Changing
c=(t<x)+(x>>19); to c=(t<=x)+(x>>19); fixes that.

But io_x points out that this overlooks the case x=c=0,
for which c=(t<=x)+(x>>19) would incorrectly make c=1
rather than 0.
This is perhaps best handled by merely, when x=0,
setting Q[j]=c; c=0; return Q[j];

Below is a revised listing of the C version;
changes should be transparent for versions
kindly provided by interested viewers---in Fortran,
PL/1 and others.

While neither correction is likely to have an
appreciable effect on use of KISS4691 as an RNG,
both corrections are necessary to ensure that
the period of MWC() will be the value (> 10^45191)
called for by the underlying mathematics of MWC,
(assuming the two period-1 seed sets,
c=0 and each Q=0
c=8192 and each Q=2^32-1
are not used).
And the results after 10^9 calls to MWC() then
10^9 call to KISS as MWC()+CNG+XS,
confirmed in various versions, are not changed.

George Marsaglia
----------------------------------------------------------
static unsigned long xs=521288629,xcng=362436069,Q[4691];

unsigned long MWC(void) /*238 million/second*/
{unsigned long t,x; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j]; if(x==0) {Q[j]=c; c=0; return 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=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);
}
---------------------------------------------------------
 
R

robin

| My posting for KISS4691 has become so mixed with
| diverse comments over various newsgroups that I
| thought it better to use a separate posting to
| bring suggested changes to interested potential
| users via sci.math, comp.lang.c, comp.lang.fortran.
|
| At least two readers have pointed out special cases that
| might arise from the code I suggested for the KISS4691.
|
| christain.bau pointed out that (x<<13)+c
| can exceed 2^32 when the rightmost 19 bits of x are 1's
| and c is its maximum possible, a-1=8192=2^13.
|
| Changing
| c=(t<x)+(x>>19); to c=(t<=x)+(x>>19); fixes that.
|
| But io_x points out that this overlooks the case x=c=0,
| for which c=(t<=x)+(x>>19) would incorrectly make c=1
| rather than 0.
| This is perhaps best handled by merely, when x=0,
| setting Q[j]=c; c=0; return Q[j];
|
| Below is a revised listing of the C version;
| changes should be transparent for versions
| kindly provided by interested viewers---in Fortran,
| PL/1 and others.

Below is my PL/I version, using unsigned arithmetic,
which I confirm delivers the same test results as before.

(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); if x = 0 then do; Q(j) = c; c = 0; return(Q(j)); end;
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 );
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;
 
R

robin

| Below is a revised listing of the C version;
| changes should be transparent for versions
| kindly provided by interested viewers---in Fortran,
| PL/1 and others.

Below is my PL/I version, using signed arithmetic,
which I confirm delivers the same test results as previously,
and which are the same results as the signed version.

(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 fixed binary (31);

if j < hbound(Q,1) then j = j + 1; else j = 0;
x = Q(j); if x = 0 then do; Q(j) = c; c = 0; return(Q(j)); end;
t = isll(x,13)+c+x;
/* Tests x >= t by computing x - t in two parts. */
t1 = isrl(x, 2) - isrl(t, 2);
if t1 = 0 then t1 = iand(x, 3) - iand(t, 3);
if t1 >= 0 then t1 = 1; else t1 = 0;
c = t1 + 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;
 
R

robin

|
|| Below is a revised listing of the C version;
|| changes should be transparent for versions
|| kindly provided by interested viewers---in Fortran,
|| PL/1 and others.
|
| Below is my PL/I version, using signed arithmetic,
| which I confirm delivers the same test results as previously,
| and which are the same results as the signed version.

um...unsigned version.
 
D

Dann Corbit

My posting for KISS4691 has become so mixed with
diverse comments over various newsgroups that I
thought it better to use a separate posting to
bring suggested changes to interested potential
users via sci.math, comp.lang.c, comp.lang.fortran.

At least two readers have pointed out special cases that
might arise from the code I suggested for the KISS4691.

christain.bau pointed out that (x<<13)+c
can exceed 2^32 when the rightmost 19 bits of x are 1's
and c is its maximum possible, a-1=8192=2^13.

Changing
c=(t<x)+(x>>19); to c=(t<=x)+(x>>19); fixes that.

But io_x points out that this overlooks the case x=c=0,
for which c=(t<=x)+(x>>19) would incorrectly make c=1
rather than 0.
This is perhaps best handled by merely, when x=0,
setting Q[j]=c; c=0; return Q[j];

Below is a revised listing of the C version;
changes should be transparent for versions
kindly provided by interested viewers---in Fortran,
PL/1 and others.

While neither correction is likely to have an
appreciable effect on use of KISS4691 as an RNG,
both corrections are necessary to ensure that
the period of MWC() will be the value (> 10^45191)
called for by the underlying mathematics of MWC,
(assuming the two period-1 seed sets,
c=0 and each Q=0
c=8192 and each Q=2^32-1
are not used).
And the results after 10^9 calls to MWC() then
10^9 call to KISS as MWC()+CNG+XS,
confirmed in various versions, are not changed.

George Marsaglia
----------------------------------------------------------
static unsigned long xs=521288629,xcng=362436069,Q[4691];

unsigned long MWC(void) /*238 million/second*/
{unsigned long t,x; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j]; if(x==0) {Q[j]=c; c=0; return Q[j];}
t=(x<<13)+c+x; c=(t<=x)+(x>>19);


Here is a slightly modified version of the code (for instance, default
types were used for c and j, which is no longer allowed in ANSI/ISO C):

-----------------------------------------------------------------
#include <stdio.h>

static unsigned long xs = 521288629, xcng = 362436069, Q[4691];

unsigned long MWC (void)
{
unsigned long t, x;
static unsigned long c = 0, j = 4691;
j = (j < 4690) ? j + 1 : 0;
x = Q[j];
if (x == 0)
{
Q[j] = c;
c = 0;
return Q[j];
}
t = (x << 13) + c + x;
c = (t <= x) + (x >> 19);
return (Q[j] = t);
}

int main (void)
{
unsigned long i, x=0;
for (i = 0; i < 4691; i++)
Q = (xcng = 69069 * xcng + 123) + (xs ^= (xs << 13), xs ^=
(xs >> 17), xs ^= (xs << 5));
for (i = 0; i < 1000000000; i++)
x = MWC ();
printf (" MWC result=3740121002 ?\n%22u\n", x);
for (i = 0; i < 1000000000; i++)
x = (MWC () + (xcng = 69069 * xcng + 123) + (xs ^= (xs << 13),
xs ^= (xs >> 17), xs ^= (xs << 5)));
printf ("KISS result=2224631993 ?\n%22u\n", x);
}
/*
Possible output:
c:\tmp>kiss5a
MWC result=3740121002 ?
3740121002
KISS result=2224631993 ?
2224631993
*/

-----------------------------------------------------------------
I think that there are two important difficulties with the code.

First, much of the calculation is performed inline. That means that we
are avoiding function call overhead, and so comparisons with fully
encapsulated PRNG functions are skewed. In order to make fully accurate
comparisons, we should have all code related to the generator fully
encapsulated in a function and called in that way.

Second, the use of static variables makes the code non-reentrant for SMP
use. It will be easy to overcome this in C++ with constructors and
class member variables, but a bit more tedious in C without the use of
some sort of gate operations.
 
T

Tim Little

Do you really have to ask that?

1. On a 64 bit processor, 64 bit types are faster than 32 bit types.

Not normally. On most 64-bit processors, 32-bit operations are
usually the same speed (or sometimes faster) than 64-bit operations.
They also produce less cache and memory bandwidth pressure.

I do agree with your overall point though: putting conceptually
separate variables at constant offset in an array is almost always a
stupid coding practice.

2. Your compiler knows that an access to an array element and an
access to a scalar struct member cannot access the same object.

An array access with suitable index can easily end up pointing at a
struct field, especially if the array is a field in the struct. Of
course, the compiler is free to ignore that as such behaviour is
undefined.

This enables many optimisations, for example a struct element can be
kept in a register throughout a lengthy computation. By using an
array, you force those registers to be written back to memory every
time another array element is accessed.

That's not actually true in practice. Modern compilers track data
access more finely than that, especially since in this case the
elements are being accessed at constant offset from a constant base.
(This doesn't stop it from being a bad idea for other reasons)


- Tim
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Members online

No members online now.

Forum statistics

Threads
473,818
Messages
2,569,732
Members
45,691
Latest member
Dick331194

Latest Threads

Top