Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.


Math Forum
»
Discussions
»
sci.math.*
»
sci.math
Notice: We are no longer accepting new posts, but the forums will continue to be readable.
Topic:
KISS4691, a potentially topranked RNG.
Replies:
108
Last Post:
Dec 6, 2010 1:29 AM



geo
Posts:
38
Registered:
6/19/08


Re: KISS4691, a potentially topranked 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 a1=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^1334071 baseb digits of the expansion of some j/p, we jumpafter an average of b=2^32 stepsto 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^n1 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



