Math Forum Discussions

geo

Posts: 38
Registered: 6/19/08

The CMWC4827 RNG, an improvement on MWC4691.
Posted: Oct 22, 2010 7:18 AM

Plain Text

In August 2010 I posted descriptions of a MWC
(Multiply-With-Carry) RNG based on, with b=2^32,
the prime p=8193b^4691-1. With an apparently
very long period and multiplier a=(2^13+1), the
basic MWC operation could be carried out using only
32-bit operations: given 32-bit c and x, imagine
t=a*x+c in 64 bits, determine the new c and new x
as the top and bottom 32-bits of t.

There were two drawbacks: rare nuisance cases that
had to be accommodated, and no proof of the order of b
for the prime p, the period of the RNG. We can be
almost---but not quite---certain that the order of
b mod p is k=(p-1)/2, since b^k mod p = 1 and for each
attainable prime divisor q of k---that is, with q one
of 431,18413, 15799501,1505986643549,
2883504568596254032007909, b^(k/q) mod p is not 1.
(My thanks to Darío Alejandro Alpern of Buenos Aires for
providing the larger, and estimates on future, factors.)

But there are at least two more, presently unattainable,
prime divisors of k, each with at least 30 digits.
With z a primitive root of p and b=z^j, the chances of j
being a multiple of one of those unknown prime divisors
of k seem less than 1/10^30. So we can be reasonably---
indeed, very---confident that the Industrial Grade Order.
(IGO) of b is (p-1)/2 for the prime p=8193b^4691-1.

So, with nuisance cases and uncertainties in mind, I put
three computers to work searching for a prime of the
form p=a*b^r+1, with a=2^i-1, to avoid the nuisance
cases arising from a=2^13+1, and p=(2^i-1)*b^r+1 so
that factoring p-1 would be feasible when r>4000.

After about 5 days, one of them found p=4095*b^4827+1.
And it turned out that the order of b mod p was the
maximum possible, k=(p-1)/2^6, (since we "lose" one 2
because 2 cannot be primitive, and must lose at least
five more because b=2^5 and k is 4095*2^154458).
Thus b^k mod p = 1, and for each prime divisor q of k,
q=2,3,5,7,13: b^(k/q) mod p is not 1.

Unlike p=ab^r-1, a prime of the form p=ab^r+1 leads
to CMWC (Complimentary-Multiply-With-Carry) RNGs, in
which we again imagine forming t=a*x+c in 64 bits and
again seek c as the top 32 bits, but rather than
x=(t mod b) for MWC, the new x is x=(b-1)-(t mod b),
that is x=~(t mod b), using C's ~ op.

Here is a C version of the resulting CMWC method for
p=4095*b^4827+1, using only 32-bit arithmetic, with
verified period 4095*2^154458, faster and simpler than
the previously posted MWC4691() and readily adapted
for either signed or unsigned integers and for Fortran
or other programming languages:
_________________________________________________________

#include <stdio.h>

static unsigned long Q[4827],carry=1271;

unsigned long CMWC4827(void)
{unsigned long t,x; static int j=4827;
j=(j<4826)? j+1:0;
x=Q[j]; t=(x<<12)+carry;
carry=(x>>20)-(t<x);
return (Q[j]=~(t-x));
}

#define CNG ( cng=69069*cng+13579 )
#define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<5) )
#define KISS4827 ( CMWC4827()+CNG+XS )

int main(void)
{unsigned long int i,x,cng=123456789,xs=362436069;
/* First seed Q[] with CNG+XS: */
for(i=0;i<4827;i++) Q[i]=CNG+XS;
/* Then generate 10^9 CMWC4827()s */
for(i=0;i<1000000000;i++) x=CMWC4827();
printf("Does x=1346668762?\n x=%lu\n",x);
/* followed by 10^9 KISS4827s: */
for(i=0;i<1000000000;i++) x=KISS4827;
printf("Does x=4041198809?\n x=%lu\n",x);

return 0;
}
______________________________________________________

Getting 10^9 CMWC4827s should take less than four seconds,
while 10^9 KISS4827s should take less than seven seconds.

For signed integers, replace that single (t<x) by (t'<x'),
where t' means flip the sign bit: t^(1<<31) for C
or ieor(t,ishft(1,-31)) for Fortran,
(or use hex 10000000 to specify the bit to be flipped).

CMWC4827() used alone will pass all tests in
the Diehard Battery of Tests of Randomness
www.cs.hku.hk/~diehard/
But as any RNG based on a single mathematical structure
is likely to have potential flaws---such as, but
perhaps not as striking as congruential RNGs "falling
mainly in the planes"---I advocate the KISS
(Keep-It-Simple-Stupid) approach: combine with
Congruential and Xorshift sequences invoked inline.

Why not splurge at a cost of 3 nanoseconds?

George Marsaglia