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: The CMWC4827 RNG, an improvement on MWC4691.
Posted:
Oct 26, 2010 6:43 PM
|
|
On Oct 25, 7:21 pm, "christian.bau" <christian....@cbau.wanadoo.co.uk> wrote: > On Oct 22, 12:18 pm, geo <gmarsag...@gmail.com> wrote: > > > CMWC4827 RNG snipped > > George Marsaglia > > I tried to figure out the maths behind this, and I ended up with a > just slightly different algorithm: > > The important number is p = 4095 * (2^32)^4827 + 1, which is a prime. > Let M = (2^32)^4827, then p = 4095 * M + 1. We can start with any > number 1 <= x < p, then replace x with (4095 x) modulo p, and again > and again. Each value x produces 4827 random integers of 32 bits. The > result x will never be 0. > > A slight change: Instead of storing x, we store x-1. So instead of > replacing x with (4095 x) modulo p, we add 1 to x, calculate (4095 x) > modulo p, and subtract 1. So we calculate ((4095 (x+1)) modulo p) - 1 > = (4095 x + 4094) modulo p. Really the same thing, but we now have 0 > <= x < 4095 * M. > > Now assume x = X + c * M, where 0 <= X < M, and 0 <= c < 4095. Then > > (4095 x + 4094) modulo p = > (4095 * (X + c * M) + 4094) modulo p = > (4095 * X + 4094 + c * (4095 M)) modulo p = > (4095 * X + 4094 + c * (p - 1)) modulo p = > (4095 * X + 4094 - c) modulo p = > 4095 * X + (4094 - c). > > We store X in an array of 4827 32-bit integers and keep c separate. > Initialiase X to random integers and c to a random value from 0 to > 4094. The random number calculation then is > > unsigned long CMWC4827(void) > { > unsigned long t, > x; > static int j = 4827; > if (j < 4826) ++j; else { j = 0; c = 4094 - c; } > x = Q[j]; > t = (x << 12) + carry; > carry = (x >> 20) - (t < x); > return (Q[j] = (t - x)); > > } > > So the change is that we replace c with 4094-c when we restart at the > beginning of the array, and don't do the "not" operation in the > calculation of the next component of x. You could change this to use > 64 bit arithmetic and write > > uint32_t CMWC4827 (void) > { > static long j = 4827; > static uint64_t carry = 1271; > if (++j >= 4827) { j = 0; carry = 4094 - carry; } > uint64_t x = Q [j] * (uint64_t) 4095 + carry; > carry = x >> 32; > return (Q [j] = x); > > > > }- Hide quoted text - > > - Show quoted text -
With b=2^32 and p=4095*b^4827+1, your suggestion will generate a sequence of 4827 "digits" forming, in reverse order, the base-b expansion of j1/(p-2) for some j1, then jump to another set of 4827 in the expansion of j2/(p-2) for some other j2, and so on. Such "digits" may well serve as satisfactory 32-bit random integers, but we cannot be sure of the period, or that we will avoid partial overlapping with string segments already generated.
If you change your suggested method so as to keep the old carry c---rather than reset c before refilling the 4827 array---and go ahead with the MWC rather than CMWC method, you will generate, in reverse order, the full base-b MWC expansion of j/m, with m=4095*b^4827-1 and period the order of b for the composite m.
The CMWC method produces, in reverse order, the full base-b=2^32 digits in the expansion of some j/p, with p=4095*b^4827+1 and j determined by the 4827 seed values and initial c in 0<=c<4095.
The order of b for the MWC method applied to m probably differs from the CMWC order, 4095*2^154458, by at most a few dozen powers of 10, and, with 10^46500, we have 46,500 in hand.
But finding that IGO, (Industrial-Grade-Order), which amounts to finding all, or all the up-to-30-digit prime factors of m=4095*b^4827-1, say p1,p2,...,pn, then IGO = lcm(order(b,p1),order(b,p2)...,order(b,pn)) is difficult.
I have considered that approach, but using m=a*b^4096-1, or even m=a*b^8192-1 with a=2^i-1, so that one can more easily effect incrementing the index of the array element. With a couple of my PCs working on finding up-to-30-digit factors of such composites, I hope to report some interesting IGOs soon.
Can any of you, with time and factoring programs available, find all---or all the up-to-30-digit---prime factors of a*b^4096-1 or a*b^8192-1 with b=2^32 and a=2^i-1, for MWC?
Or the same for CMWC primes or composites, a*b^4096+1 or a*b^8192+1?
George Marsaglia
|
|
|
|