Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: KISS4691, a potentially top-ranked RNG.
Replies: 108   Last Post: Dec 6, 2010 1:29 AM

 Messages: [ Previous | Next ]
 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

Date Subject Author
7/24/10 geo
7/24/10 Geoff
7/24/10 jacob navia
7/27/10 robin
7/30/10 Uno
7/24/10 Gene
7/24/10 Gib Bogle
7/25/10 geo
7/25/10 Dick Hendrickson
7/25/10 Gib Bogle
7/26/10 robin
7/26/10 Geoff
7/26/10 Gib Bogle
7/26/10 robin
7/27/10 Geoff
7/27/10 nmm1@cam.ac.uk
7/26/10 Gib Bogle
7/29/10 Uno
7/29/10 Gib Bogle
7/29/10 Uno
7/30/10 Gib Bogle
7/30/10 glen herrmannsfeldt
7/30/10 Gib Bogle
7/30/10 robin
7/30/10 robin
7/30/10 Gib Bogle
7/25/10 robin
7/26/10 Harald Anlauf
7/26/10 nmm1@cam.ac.uk
7/26/10 glen herrmannsfeldt
7/29/10 Uno
7/29/10 glen herrmannsfeldt
7/29/10 Ron Shepard
7/29/10 glen herrmannsfeldt
7/30/10 Uno
7/30/10 nmm1@cam.ac.uk
7/30/10 Ron Shepard
7/31/10 Ilmari Karonen
7/31/10 Uno
7/31/10 sturlamolden
7/30/10 David Bernier
7/30/10 orz
7/30/10 glen herrmannsfeldt
8/1/10 Ron Shepard
8/1/10 glen herrmannsfeldt
7/30/10 orz
7/30/10 orz
7/30/10 Gib Bogle
7/30/10 orz
7/30/10 Gib Bogle
7/31/10 orz
7/31/10 Gib Bogle
8/4/10 orz
8/5/10 Uno
8/10/10 Gib Bogle
8/10/10 orz
8/10/10 Gib Bogle
8/10/10 Dann Corbit
8/11/10 Richard Maine
8/28/10 robin
8/12/10 robin
8/17/10 Gib Bogle
8/17/10 glen herrmannsfeldt
8/17/10 gnasher729
8/18/10 geo
8/18/10 gnasher729
8/19/10 gnasher729
8/19/10 geo
8/19/10 glen herrmannsfeldt
8/19/10 orz
8/19/10 orz
8/21/10 Guest
8/21/10 Guest
8/17/10 Gib Bogle
8/18/10 glen herrmannsfeldt
8/18/10 Gib Bogle
8/18/10 baf
8/19/10 Gib Bogle
8/19/10 baf
8/19/10 Gib Bogle
8/18/10 Gib Bogle
8/28/10 robin
8/28/10 robin
8/10/10 orz
8/11/10 glen herrmannsfeldt
12/6/10 Michael Press
7/26/10 Harald Anlauf
7/26/10 nmm1@cam.ac.uk
7/27/10 robin
7/30/10 Uno
8/3/10 robin
8/1/10 sturlamolden
8/2/10 Georg Beyerle
8/2/10 Uno
8/2/10 sturlamolden
8/3/10 Dann Corbit
8/3/10 sturlamolden
8/4/10 orz