cametan001

#lang racket

(require srfi/41)

;;   32-bits Random number generator U(0,1): MRG32k3a
;;   Author: Pierre L'Ecuyer,
;;   Source: Good Parameter Sets for Combined Multiple Recursive Random
;;           Number Generators,
;;           Shorter version in Operations Research,
;;           47, 1 (1999), 159--164.

cametan001

> (stream-ref s 100)
0.07967177768510994
> (stream-ref s 1000)
0.9014158012565428
> (stream-ref s 10000)
0.1110427223837204

cametan001

> (define s (Combination (stream1 SEED SEED SEED)
                         (stream2 SEED SEED SEED)))
> s
#<stream>
> (stream->list 10 s)
'(0.12701112204657714
  0.3185275653967945
  0.3091860155832701
  0.8258468629271136
  0.2216299157820229

cametan001

(define (Combination stream1 stream2)
  (let ((stream1 (stream-cdr (stream-cdr (stream-cdr stream1))))
        (stream2 (stream-cdr (stream-cdr (stream-cdr stream2)))))
    (stream-map (lambda (p1 p2)
                  (* (- p1 (if (<= p1 p2)
                               (- p2 m1)
                               p2)) norm))
                stream1 stream2)))

cametan001

(define (stream2 s20 s21 s22)
  (stream-cons s20 (stream2 s21 s22 (Component2 s20 s22))))

(define (Component2 s20 s22)
  (let ((p2 (modulo (- (* a21 s22) (* a23n s20)) m2)))
    (if (negative? p2)
        (+ p2 m2)
        p2)))

cametan001

> (define s1 (stream1 SEED SEED SEED))
> s1
#<stream>
> (stream->list 10 s1)
'(12345
  12345
  12345
  3023790853.0
  3023790853.0
  3385359573.0

cametan001

(define (Component1 s10 s11)
  (let ((p1 (modulo (- (* a12 s11) (* a13n s10)) m1)))
    (if (negative? p1)
        (+ p1 m1)
        p1)))

cametan001

   p1 = a12 * s11 - a13n * s10;
   k = p1 / m1;
   p1 -= k * m1;
   if (p1 < 0.0)
      p1 += m1;

cametan001

(define (stream1 s10 s11 s12)
  (stream-cons s10 (stream1 s11 s12 (Component1 s10 s11))))

cametan001

   /* Component 1 */
   p1 = a12 * s11 - a13n * s10;
   k = p1 / m1;
   p1 -= k * m1;
   if (p1 < 0.0)
      p1 += m1;
   s10 = s11;
   s11 = s12;
   s12 = p1;