* Re: Implementing the KISS4691 RNG
[not found] <93b6ca11-f4ae-4598-98f8-9a10be69ca3a@j18g2000yqd.googlegroups.com>
@ 2010-09-07 16:51 ` robin
2010-09-08 3:15 ` robin
1 sibling, 0 replies; 2+ messages in thread
From: robin @ 2010-09-07 16:51 UTC (permalink / raw)
"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;
^ permalink raw reply [flat|nested] 2+ messages in thread
* Re: Implementing the KISS4691 RNG
[not found] <93b6ca11-f4ae-4598-98f8-9a10be69ca3a@j18g2000yqd.googlegroups.com>
2010-09-07 16:51 ` Implementing the KISS4691 RNG robin
@ 2010-09-08 3:15 ` robin
1 sibling, 0 replies; 2+ messages in thread
From: robin @ 2010-09-08 3:15 UTC (permalink / raw)
"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;
^ permalink raw reply [flat|nested] 2+ messages in thread
end of thread, other threads:[~2010-09-08 3:15 UTC | newest]
Thread overview: 2+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
[not found] <93b6ca11-f4ae-4598-98f8-9a10be69ca3a@j18g2000yqd.googlegroups.com>
2010-09-07 16:51 ` Implementing the KISS4691 RNG robin
2010-09-08 3:15 ` robin
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox