//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author John Miller
* @version 1.2
* @date Sat Mar 22 14:39:30 EDT 2014
* @see LICENSE (MIT style license file).
*
* CMRG (Combined Multiple Recursive Generator) using 32-bit Int's
*/
package scalation.random
import scalation.util.time
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `Random2` class generates random real numbers in the range (0, 1).
* It implements, using 32-bit integers (Int's), the 'MRG31k3p' generator
* developed by L'Ecuyer and Touzin, described in "FAST COMBINED MULTIPLE
* RECURSIVE GENERATORS WITH MULTIPLIERS OF THE FORM a = 2^q +/- 2^r".
* MRG31k3p is a Combined Multiple Recursive Generator (CMRG) shown to have good
* performance and statistical properties for simulations. It has a period of
* about 2^185 and is considered to be a faster alternative to the popular
* 'MRG32k3' generator. MRG31k3p combines MRG1 and MRG2.
*
* MRG1: x_i = (0 + a_12 x_i-2 + a_13 x_i-3) % M1
* MRG2: x_i = (a_21 x_i-1 + 0 + a_23 x_i-3) % M2
*
* where a_12 = 2^22, a_13 = 2^7+1, a_21 = 2^15 and a_23 = 2^15+1.
* @see http://www.informs-sim.org/wsc00papers/090.PDF
* @see http://www.iro.umontreal.ca/~simardr/ssj/doc/pdf/guiderng.pdf
* @param stream the random number stream index
*/
case class Random2 (stream: Int = 0)
extends RNG (stream)
{
private val M1 = 2147483647 // modulus for MRG1 (2^31 - 1)
private val M2 = 2147462579 // modulus for MRG2 (2^31 - 21069)
private val MASK12 = 511 // mask to extract 9 lsb's (2^9 - 1)
private val MASK13 = 16777215 // mask to extract 24 lsb's (s^24 - 1)
private val MASK21 = 65535 // mask to extract 16 lsb's (s^16 - 1)
private val MULT2 = 21069 // multiplier
private val NORM = 4.656612873077393e-10 // 1.0 / (M1 + 1.0) normalization to (0, 1)
private val x = RandomSeeds.seeds (stream) // 6-dimensional vector of seed values
private var x11 = x(0) // x_i-1 for MRG1
private var x12 = x(1) // x_i-2 for MRG1
private var x13 = x(2) // x_i-3 for MRG1
private var x21 = x(3) // x_i-1 for MRG2
private var x22 = x(4) // x_i-2 for MRG2
private var x23 = x(5) // x_i-3 for MRG2
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** Return the next random number as a real `Double` in the interval (0, 1).
* This calculation uses 32-bit integers `Int`.
*/
def gen: Double =
{
// Calculate MRG1 (first component)
var y1 = ((x12 & MASK12) << 22) + (x12 >>> 9) + ((x13 & MASK13) << 7) + (x13 >>> 24)
if (y1 < 0 || y1 >= M1) y1 -= M1 //must also check overflow
y1 += x13
if (y1 < 0 || y1 >= M1) y1 -= M1
x13 = x12; x12 = x11; x11 = y1
// Calculate MRG2 (second component)
y1 = ((x21 & MASK21) << 15) + (MULT2 * (x21 >>> 16))
if (y1 < 0 || y1 >= M2) y1 -= M2
var y2 = ((x23 & MASK21) << 15) + (MULT2 * (x23 >>> 16))
if (y2 < 0 || y2 >= M2) y2 -= M2
y2 += x23
if (y2 < 0 || y2 >= M2) y2 -= M2
y2 += y1
if (y2 < 0 || y2 >= M2) y2 -= M2
x23 = x22; x22 = x21; x21 = y2
// Combine MRG1 and MRG2 (must never return either 0 or M1+1)
if (x11 <= x21) (x11 - x21 + M1) * NORM else (x11 - x21) * NORM
} // gen
//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** Return the next stream value as an integer.
* This calculation uses 32-bit integers `Int`.
*/
def igen: Int =
{
// Calculate MRG1 (first component)
var y1 = ((x12 & MASK12) << 22) + (x12 >>> 9) + ((x13 & MASK13) << 7) + (x13 >>> 24)
if (y1 < 0 || y1 >= M1) y1 -= M1 //must also check overflow
y1 += x13
if (y1 < 0 || y1 >= M1) y1 -= M1
x13 = x12; x12 = x11; x11 = y1
// Calculate MRG2 (second component)
y1 = ((x21 & MASK21) << 15) + (MULT2 * (x21 >>> 16))
if (y1 < 0 || y1 >= M2) y1 -= M2
var y2 = ((x23 & MASK21) << 15) + (MULT2 * (x23 >>> 16))
if (y2 < 0 || y2 >= M2) y2 -= M2
y2 += x23
if (y2 < 0 || y2 >= M2) y2 -= M2
y2 += y1
if (y2 < 0 || y2 >= M2) y2 -= M2
x23 = x22; x22 = x21; x21 = y2
// Combine MRG1 and MRG2 (must never return either 0 or M1+1)
if (x11 <= x21) x11 - x21 + M1 else x11 - x21
} // igen
} // Random2 class