G

#### geo

period I cannot determine---nor is it likely that

anyone else can---but that unknown period is almost

certainly greater than 10^40000000, i.e, 10^(40million).

Versions for either 32- and 64-bit integers are given,

as well as suggestions for similar versions with

periods as great or even greater magnitude.

MWC RNGs are ordinarily based on primes of the form

p=ab^r-1, with b=2^32 or 2^64 and a<b. They generate,

32- or 64-bits at a time, and in reverse order, the

binary expansions of rationals j/p, 0<j<p and the

period of that binary expansion is the order of 2 mod p.

The key part of the MWC process is: given 32-bit a,x,c,

extract the top and bottom 32 bits of a 64 bit a*x+c,

or, for given 64-bit a,x,c, extract the top and bottom

64 bits of a 128-bit a*x+c.

I have put various computers to work establishing the

order of 2 for MWC primes of the form ab^r-1, b=2^32,

or CMWC (Complimentary-Multiply-With-Carry), p=ab^r+1.

For example, it took 24 days on an early laptop to

establish that the order of 2 for p=54767*2^1337287+1

is (p-1)/2^3, and just last week I put a Samsung 64-bit

laptop to finding the order of 2 for the 11th largest

known prime, p=27653*2^9167433+1; in CMWC form,

p=14158336*b^286482+1. I expect it may take months.

Such searches are motivated by more than a Mt. Everest

syndrome, as extremely-long-period RNGs not only seem to

perform very well on tests of randomness, but have uses

for cryptography. For example, suppose, from a random

set of seeds, you have observed 100000 successive 32-bit

numbers from a CMWC RNG based on p=14158336*b^286482+1.

You will then be somewhere in an immense loop that has

over 10^(2million) appearances of exactly that same

sequence of 100000, leaving you virtually no chance of

finding your particular location and thus able to

predict forthcoming values.

Of course, if you can observe 286482 successive values,

then, after a little modular arithmetic to determine

the carry, you are OK. It is ensuring a large exponent

for b that makes the RNG desirable---the larger the better.

But in seeking large exponents for b and record periods,

why bother with an Everest when there are ranges with

peaks far beyond surveyor's skills but certain

to be several times as tall?

Such ranges come from considering, for very large r,

composites of the form m=ab^r-1, rather than primes of

that form. The example for this posting: the composite

m=(2^28-1)*2^(2^27)-1=(2^28-1)*b^(2^22)-1.

m has no prime divisors less than 21000000, and order(2,m)

is the lcm of the orders mod the prime divisors of m.

We cannot---and may never be able to--factor such an

m>10^40403570, but for primes q, the average value of

order(2,q)/q is around .572. Thus, even if m has, say,

40 prime factors, since (.572)^40>10^(-10), we can only

expect to "lose" around ten 10's, and perhaps ten more

through common factors for the lcm, providing an estimate

around 10^40403550 for the order of 2 mod m.

That the order of 2 mod m is less than 10^(40million)

seems extremely unlikely; we would have to "lose" over

400thousand 10's, when something in the range

of three to sixty 10's seems more likely.

So users should feel confident that the periods of the

following MWC RNGs far exceed 10^(40million). Based on

m=(2^28-1)b^(2^22)-1=(2^28-1)b^4194304-1 with b=2^32,

m=(2^28-1)B^(2^21)-1=(2^28-1)B^2091752-1 with B=2^64,

they are simple and very fast.

You would have to observe more than 134 million

successive bits produced be these MWC RNGs

in order to get close to cracking them.

Note that choosing a=2^i-1 in m=ab^r-1 makes finding

the top and bottom 32 bits of a purported 64 bit

a*x+c feasible using only 32-bit arithmetic,

or, within 64-bit arithmetic, finding the top

and bottom 64 bits of a purported 128 bit a*x+c.

Here are C versions using, for 32-bit integers,

an unsigned long array Q[41943104], and for 64-bits,

an unsigned long long array Q[2091752].

Try them and see if you get the same results I do.

--------------------32-bit MWC version-------------------

#include <stdio.h>

static unsigned long Q[4194304],carry=0;

unsigned long b32MWC(void)

{unsigned long t,x; static int j=4194303;

j=(j+1)&4194303;

x=Q[j]; t=(x<<28)+carry;

carry=(x>>4)-(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 KISS ( b32MWC()+CNG+XS )

int main(void)

{unsigned long int i,x,cng=123456789,xs=362436069;

/* First seed Q[] with CNG+XS: */

for(i=0;i<4194304;i++) Q

*=CNG+XS;*

/* Then generate 10^9 b32MWC()s */

for(i=0;i<1000000000;i++) x=b32MWC();

printf("Does x=2769813733 ?\n x=%lu\n",x);

/* followed by 10^9 KISSes: */

for(i=0;i<1000000000;i++) x=KISS;

printf("Does x=3545999299 ?\n x=%lu\n",x);

return 0;

}

---------------- 64-bit MWC version ---------------------

#include <stdio.h>

static unsigned long long Q[2097152],carry=0;

unsigned long long B64MWC(void)

{unsigned long long t,x; static int j=2097151;

j=(j+1)&2097151;

x=Q[j]; t=(x<<28)+carry;

carry=(x>>36)-(t<x);

return (Q[j]=t-x);

}

#define CNG ( cng=6906969069LL*cng+13579 )

#define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<43) )

#define KISS ( B64MWC()+CNG+XS )

int main(void)

{unsigned long long i,x,cng=123456789987654321LL,

xs=362436069362436069LL;

/* First seed Q[] with CNG+XS: */

for(i=0;i<2097152;i++) Q

/* Then generate 10^9 b32MWC()s */

for(i=0;i<1000000000;i++) x=b32MWC();

printf("Does x=2769813733 ?\n x=%lu\n",x);

/* followed by 10^9 KISSes: */

for(i=0;i<1000000000;i++) x=KISS;

printf("Does x=3545999299 ?\n x=%lu\n",x);

return 0;

}

---------------- 64-bit MWC version ---------------------

#include <stdio.h>

static unsigned long long Q[2097152],carry=0;

unsigned long long B64MWC(void)

{unsigned long long t,x; static int j=2097151;

j=(j+1)&2097151;

x=Q[j]; t=(x<<28)+carry;

carry=(x>>36)-(t<x);

return (Q[j]=t-x);

}

#define CNG ( cng=6906969069LL*cng+13579 )

#define XS ( xs^=(xs<<13), xs^=(xs>>17), xs^=(xs<<43) )

#define KISS ( B64MWC()+CNG+XS )

int main(void)

{unsigned long long i,x,cng=123456789987654321LL,

xs=362436069362436069LL;

/* First seed Q[] with CNG+XS: */

for(i=0;i<2097152;i++) Q

*=CNG+XS;*

/* Then generate 10^9 B64MWC()s */

for(i=0;i<1000000000;i++) x=B64MWC();

printf("Does x=13596816608992115578 ?\n x=%llu\n",x);

/* followed by 10^9 KISSes: */

for(i=0;i<1000000000;i++) x=KISS;

printf("Does x=5033346742750153761 ?\n x=%llu\n",x);

return 0;

}

---------------------------------------------------------

Note that I advocate the KISS, (Keep-It-Simple-Stupid),

principle. Even though the MWC RNGs perform very well on

tests of randomness, combining with Congruential (CNG)

and Xorshift (XS) probably provides extra insurance

at the cost of a few simple instructions, (and making

the resulting KISS even harder to crack).

Also note that for unsigned integers such as in Fortran,

the C instruction -(t<x) means: subtract 1 if t<s,

and can be implemented for signed integers as t'<s',

where the ' means: flip the sign bit.

Full random seeding of the carry and the Q[] array

calls for more than 17 megabytes of random bits.

For ordinary randomness testing, filling, as above,

the Q[] array with Congruential+Xorshift may serve, but

for stronger crypto uses, many more than the random 64-

or 128-bits that CNG and XS require are needed.

Random seeding of the carry and the Q[] array amounts to

producing, either 32- or 64-bits at a time, the reverse-

order binary expansion of j/m for some 0<=j<=m. Indeed,

for m=ab^r-1, given the carry and the Q[] array, that j is

j=carry+a*sum(Q/* Then generate 10^9 B64MWC()s */

for(i=0;i<1000000000;i++) x=B64MWC();

printf("Does x=13596816608992115578 ?\n x=%llu\n",x);

/* followed by 10^9 KISSes: */

for(i=0;i<1000000000;i++) x=KISS;

printf("Does x=5033346742750153761 ?\n x=%llu\n",x);

return 0;

}

---------------------------------------------------------

Note that I advocate the KISS, (Keep-It-Simple-Stupid),

principle. Even though the MWC RNGs perform very well on

tests of randomness, combining with Congruential (CNG)

and Xorshift (XS) probably provides extra insurance

at the cost of a few simple instructions, (and making

the resulting KISS even harder to crack).

Also note that for unsigned integers such as in Fortran,

the C instruction -(t<x) means: subtract 1 if t<s,

and can be implemented for signed integers as t'<s',

where the ' means: flip the sign bit.

Full random seeding of the carry and the Q[] array

calls for more than 17 megabytes of random bits.

For ordinary randomness testing, filling, as above,

the Q[] array with Congruential+Xorshift may serve, but

for stronger crypto uses, many more than the random 64-

or 128-bits that CNG and XS require are needed.

Random seeding of the carry and the Q[] array amounts to

producing, either 32- or 64-bits at a time, the reverse-

order binary expansion of j/m for some 0<=j<=m. Indeed,

for m=ab^r-1, given the carry and the Q[] array, that j is

j=carry+a*sum(Q

**b^i),i=0..r-1.*

There are phi(m) such j's with period order(2,m).

Since m has no prime factors less than 21000000, with

probability close to 1 a random seeding of Q[] should

produce a j that is relatively prime to m and thus

produce the full period---unknown but almost certainly

far in excess of 10^(40million). Even if random

seeding of carry and Q[] provides a j with gcd(j,m)>1,

the period is still likely to exceed 10^(40million).

There are two exceptions:

All QThere are phi(m) such j's with period order(2,m).

Since m has no prime factors less than 21000000, with

probability close to 1 a random seeding of Q[] should

produce a j that is relatively prime to m and thus

produce the full period---unknown but almost certainly

far in excess of 10^(40million). Even if random

seeding of carry and Q[] provides a j with gcd(j,m)>1,

the period is still likely to exceed 10^(40million).

There are two exceptions:

All Q

*=0, carry=0 yields j=0 and the binary expansion*

0/m=.00000000000000000...

All Q0/m=.00000000000000000...

All Q

*=b-1=2^32-1 and carry=a-1=(2^28-2) yields j=m*

and the binary expansion

m/m=.11111111111111111...

Almost any choice of m=(2^i-1)*2^(2^27)-1 is likely

to provide an immense order for 2 mod m.

Choices i=16,18,22,28 resulted in m's that have no

factors<10^7. Memory seems to be the key restriction

for the resulting Q[] array. With enough memory,

one might seek i's for which m=(2^i-1)*2^(2^28)-1

has no small factors, or even m=(2^i-1)*2^(2^29)-1,

boosting the unknown and unknowable periods of

MWC RNGs to ranges beyond 10^(40million).

George Marsagliaand the binary expansion

m/m=.11111111111111111...

Almost any choice of m=(2^i-1)*2^(2^27)-1 is likely

to provide an immense order for 2 mod m.

Choices i=16,18,22,28 resulted in m's that have no

factors<10^7. Memory seems to be the key restriction

for the resulting Q[] array. With enough memory,

one might seek i's for which m=(2^i-1)*2^(2^28)-1

has no small factors, or even m=(2^i-1)*2^(2^29)-1,

boosting the unknown and unknowable periods of

MWC RNGs to ranges beyond 10^(40million).

George Marsaglia