Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Marsaglia RNG in Mathematica

• To: mathgroup at smc.vnet.net
• Subject: [mg38973] Re: [mg38958] Marsaglia RNG in Mathematica
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Thu, 23 Jan 2003 08:04:14 -0500 (EST)
• References: <200301221114.GAA05344@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

flip wrote:
> 
> Hello All,
> 
> I was wondering if someone here would be nice enough to convert this C to
> Mathematica code.  A module would be best I suppose, so users can just call
> it just like Random[].
> 
> I'd like to experiment with this RNG in Mathematica.  In fact, IIRC Mathematica uses some variant of SWC also.
> 
> Thank you to anyone that can help.
> 
> Flip    (remove "_alpha" to email me directly)
> 
> By the way, the following words are from another thread in another NG
> regarding a RNG.
> 
> ***  Words below from Dr. Marsaglia ***
> 
> Here is a MWC example with period 3056868392^33216-1, a mere
> 10^4005 times as long as that of the Mersenne twister, yet faster,
> far simpler, and seemingly at least as well-performing in tests
> of randomness:
> 
> static unsigned long Q[1038],c=362436;
> 
> unsigned long MWC1038(void){
> static unsigned long i=1037;
> unsigned long long t, a=611373678LL;
> t=a*Q[i]+c; c=(t>>32);
> if(--i) return(Q[i]=t);
> i=1037; return(Q[0]=t);
>                            }
> 
> You need to assign random 32-bits seeds to the static array Q[1038]
> and to the initial carry c, with 0<=c<61137367.
> 
> George Marsaglia

I think the code below does what is intended. This assumes a long long
is 64 bits and a long is 32 bits. These are both quite platform
dependent for C; maybe Marsaglia specified a particular platform? (a URL
to his note would be helpful). For seeding the q and c variables I use
built-in Random. I also changed the variable names to avoid single
characters and capitalized names.

SeedRandom[1111];
qarray = Table[Random[Integer,{0,2^32-1}], {1038}];
cc = Random[Integer, {0,61137367}];
index = 1037;
aconst = 611373678;

mwc1038[] := Module[
  {},
  tt = aconst*qarray[[index]] + cc;
  cc = Developer`BitShiftRight[tt,32]; (*Could instead use Quotient *)
  tt = Mod[tt,2^32];
  index--;
  If [index!=0, qarray[[index]] = tt; Return[tt]];
  index = 1037;
  qarray[[0]] = tt
  ]

Example:

Table[mwc1038[], {10}]

Out[48]= {1392539486, 2580953375, 977837705, 1393978643, 2243451427, 
  374027900, 2689334407, 2953422566, 41924440, 2262165557}


Daniel Lichtblau
Wolfram Research

• References:
  • Marsaglia RNG in Mathematica
    • From: "flip" <flip_alpha@safebunch.com>