Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
Drexel University or The Math Forum.
|
|
geo
Posts:
38
Registered:
6/19/08
|
|
Re: KISS4691, a potentially top-ranked RNG.
Posted:
Aug 19, 2010 9:11 AM
|
|
On Aug 18, 6:08 pm, "christian.bau" <christian....@cbau.wanadoo.co.uk>
wrote:
> Hi George,
>
> the problem is not that the carry could match or exceed the multiplier
> 8193.
> The problem is that the carry can be as large as 8192.
>
> In your code you write
>
> t=(x<<13)+c+x; c=(t<x)+(x>>19);
>
> To make it a bit more obvious what happens, I can change this to the
> 100% equivalent
>
> t = (x << 13) + c;
> c = x >> 19;
> t += x; if (t < x) ++c;
>
> The last line is the usual trick to add numbers and check for carry: t
> will be less than x exactly if the addition t += x wrapped around and
> produced a carry.
>
> The problem is that (x << 13) + c can also produce a carry, and that
> goes unnoticed. If c = 8192, and the lowest 19 bits in x are all set,
> then (x << 13) has the highest 19 bits set and 13 bits zero. Adding c
> "overflows" and gives a result of 0, and adding x leaves x. The
> comparison (t < x) is false and produces zero, even though (x<<13)+c+x
> should have produced a carry.
>
> This is very, very rare. It doesn't actually happen in your original
> source code with a total of 2 billion random numbers. If you produce
> four billion random numbers, then it happens at i = 3596309491 and
> again at i = 3931531774.
>
> I am more interested in 64 bit performance, so I just made c, t, and x
> 64 bit numbers and wrote
>
> uint64_t x = Q [i];
> uint64_t t = (x << 13) + x + c;
> c = (t >> 32);
>
> That might help in Fortran as well, assuming 64 bit integers are
> available, signed or unsigned doesn't matter.
Thanks for pointing that out.
The correction (t<x) should be (t<=x),
to take care of that rare, but inevitable, event
when the last 19 bits of x are 1's
and the carry c is a-1=2^13, making (x<<13)+c = 0,
and thus t=(x<<13)+c+x = x,
(mod 2^32, of course, the assumed underlying
integer arithmetic of the processor).
The entire MWC() routine should have that
(t<x) to (t<=x) correction:
unsigned long MWC(void)
{ unsigned long t,x,i; static c=0,j=4691;
j=(j<4690)? j+1:0;
x=Q[j];
t=(x<<13)+c+x; c=(t<=x)+(x>>19);
return (Q[j]=t);
}
and for interested readers who may have pasted
the original, please change that (t<x) to (t<=x).
The test values from 10^9 calls will still be as indicated,
but strict adherence to the underlying mathematics requires
the change from (t<x) to (t<=x) to ensure the full period.
Keeping that (t<x) may be viewed as, rather than making a
circular transition through the 8193*2^133407-1 base-b digits
of the expansion of some j/p, we jump---after an average of
b=2^32 steps---to another point in the immense circle.
It could be that the period is increased rather than
decreased, but it remains a curiosity. Perhaps light could
be shed with primes p=ab^n-1 for which the actual distribution
of such jumps, averaging every b steps, could be examined.
Or perhaps it could remain a mystery
and be further fodder for those who tend to
equate lack of understanding to "true" randomness.
George Marsaglia
|
|
|
|
