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 18, 2010 8:25 AM

On Aug 17, 6:07 pm, "christian.bau" <christian....@cbau.wanadoo.co.uk>
wrote:
> I think there is a problem in the C code.
>
> MWC () is supposed to implement multiplication of a huge number by
> 8193, with the huge number split into 32 bit items.
> We could implement this by calculating (Q[i] * 8193 + carry) with 64
> bit arithmetic, storing 32 bits, moving the upper 32 bits into carry.
> In 32 bit arithmetic, we calculate the 32 bits to be stored as (Q[i] +
> (Q[i]<<13) + carry), that part is fine.
> The new carry should be (Q[i] >> 19), plus the carry in adding (Q[i] +
> (Q[i]<<13) + carry) in 32 bit arithmetic.
>
> The carry is found by adding all three items, and checking whether the
> sum is less than Q [i].
> That is incorrect when (Q[i]<<13) + carry produces a carry in 32 bit
> arithmetic.
> And that happens when we start with carry = 0x800 and Q [i] =
> 0xffffffff.
>
> I have no idea how this would influence the randomness, but it changes
> the random number generator from something that can be examined
> mathematically to something that is quite random.

I commend you on your analysis rather than concern over
programming style, but in the MWC method, with modulus b,(=2^32)
multiplier a<b, prime p=ab^n-1, (n=4196)
and 0 <= x < b and 0 <= c < a,

new x = a*x+c mod b, new c = integerpartof((a*x+c)/b)

and the new c can never reach, or exceed, a.

Formally, we may consider the set Z of a*b^4169 "vectors" z:

z=[c;x_1,x_2,...,x_{4169}], 0 <= c < a and 0 <= x_i < 2^32

and a function f(z): with t=a*x+c,

f([c;x_1,x_2,...,x_{4169}])=
[trunc(t/b);x_2,x_3,...,x_{4169},t mod b],

If k is the order of b mod p, the directed graph:z -> f(z)
has (p-1)/k loops of size k and two "loops" of size 1:

f(z)=z for z=[0;0,0,...,0] and z=[b-1,a-1,a-1,...,a-1].

The concatenated sequence made by taking the rightmost x
from each z forms, in reverse order,
the base-b expansion of j/p for some 0 <= j <= p.

When z=[0;0,0,...,0]
we get 0/p=.0000000...
When z=[a-1;d,d,d,...,d] with d=b-1,
we get p/p=.dddd...

In our case, we have four possible loops,
two of size (p-1)/2, two of size 1,
for a total of 2(p-1)/2 + 2 = p+1 = ab^4169,
related by the mappings
c <-> a-1-c, x <-> b-1-x, j/p <-> (p-j)/p.

The function f has an inverse: with s=c*b+x_{4169},

f^{-1}(z)=[s mod a;trunc(s/a),x_1,x_2,...,x_{4168}],

so that with g(z)=f^{-1}(z),
the trailing element in each of the vectors
g(z),g(g(z)),g(g(g(z))),...
forms, in normal order, the base-b expansion of some j/p.
But arithmetic in g(z) is not well-suited for computer
implementation, while that of f(z) is eminently so, the
reason I created it---the method, not the arithmetic.

By finding the prime p=ab^n-1 with
b=2^32, a=8193=2^13+1 and n=4169,
I as able to ensure that the necessary arithmetic
could be simply---and quickly--carried out with
the base b=2^32 arithmetic available to most of us.

In my original post I neglected to point out that the
initial choice of c should satisfy 0 <= c < a=8193
and that, choosing
c=0 and all the Q's zero, or
c=a-1 and all the Q's b-1
would makes the MWC generator have period 1
and the resulting KISS period a mere 2^64-2^32
rather than one exceeding 10^40000.

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