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 topranked 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^n1, (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 (p1)/k loops of size k and two "loops" of size 1:
f(z)=z for z=[0;0,0,...,0] and z=[b1,a1,a1,...,a1].
The concatenated sequence made by taking the rightmost x from each z forms, in reverse order, the baseb expansion of j/p for some 0 <= j <= p.
When z=[0;0,0,...,0] we get 0/p=.0000000... When z=[a1;d,d,d,...,d] with d=b1, we get p/p=.dddd...
In our case, we have four possible loops, two of size (p1)/2, two of size 1, for a total of 2(p1)/2 + 2 = p+1 = ab^4169, related by the mappings c <> a1c, x <> b1x, j/p <> (pj)/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 baseb expansion of some j/p. But arithmetic in g(z) is not wellsuited for computer implementation, while that of f(z) is eminently so, the reason I created itthe method, not the arithmetic.
By finding the prime p=ab^n1 with b=2^32, a=8193=2^13+1 and n=4169, I as able to ensure that the necessary arithmetic could be simplyand quicklycarried 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=a1 and all the Q's b1 would makes the MWC generator have period 1 and the resulting KISS period a mere 2^642^32 rather than one exceeding 10^40000.
Please forgive me.
George Marsaglia



