make-MRG32k3a.rkt

#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.
;;   Ported to Scheme by Cametan
;;   ---------------------------------------------------------

(define (make-MRG32k3a)
  (let ((norm 2.328306549295728e-10)
        (m1   4294967087.0)
        (m2   4294944443.0)
        (a12     1403580.0)
        (a13n     810728.0)
        (a21      527612.0)
        (a23n    1370589.0)
        
        
        (SEED 12345))
    
    (letrec ((stream1
              (lambda (s10 s11 s12)
                (stream-cons s10 (stream1 s11 s12 (Component1 s10 s11)))))
             
             (Component1
              (lambda (s10 s11)
                (let ((p1 (modulo (- (* a12 s11) (* a13n s10)) m1)))
                  (if (negative? p1)
                      (+ p1 m1)
                      p1))))
             
             (stream2
              (lambda (s20 s21 s22)
                (stream-cons s20 (stream2 s21 s22 (Component2 s20 s22)))))
             
             (Component2
              (lambda (s20 s22)
                (let ((p2 (modulo (- (* a21 s22) (* a23n s20)) m2)))
                  (if (negative? p2)
                      (+ p2 m2)
                      p2))))
             
             (Combination
              (lambda (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)))))
      
      (Combination (stream1 SEED SEED SEED)
                   (stream2 SEED SEED SEED)))))