G

#### geo

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

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