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 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.
Please forgive me.
George Marsaglia
|
|
|
|
