G

#### geo

(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

*=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

/* 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