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


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


In August 2010 I posted descriptions of a MWC (MultiplyWithCarry) RNG based on, with b=2^32, the prime p=8193b^46911. With an apparently very long period and multiplier a=(2^13+1), the basic MWC operation could be carried out using only 32bit operations: given 32bit c and x, imagine t=a*x+c in 64 bits, determine the new c and new x as the top and bottom 32bits 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 almostbut not quitecertain that the order of b mod p is k=(p1)/2, since b^k mod p = 1 and for each attainable prime divisor q of kthat 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, veryconfident that the Industrial Grade Order. (IGO) of b is (p1)/2 for the prime p=8193b^46911.
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^i1, to avoid the nuisance cases arising from a=2^13+1, and p=(2^i1)*b^r+1 so that factoring p1 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=(p1)/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^r1, a prime of the form p=ab^r+1 leads to CMWC (ComplimentaryMultiplyWithCarry) 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=(b1)(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 32bit 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]=~(tx)); }
#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 flawssuch as, but perhaps not as striking as congruential RNGs "falling mainly in the planes"I advocate the KISS (KeepItSimpleStupid) approach: combine with Congruential and Xorshift sequences invoked inline.
Why not splurge at a cost of 3 nanoseconds?
George Marsaglia



