;; from: GNU Scientific Library
;;
;; Generator: gsl_rng_cmrg
;; This is a combined multiple recursive generator by L'Ecuyer. Its
;; sequence is,
;;
;;   z_n = (x_n - y_n) mod m_1
;;
;; where the two underlying generators x_n and y_n are,
;;
;;   x_n = (a_1 x_{n-1} + a_2 x_{n-2} + a_3 x_{n-3}) mod m_1
;;   y_n = (b_1 y_{n-1} + b_2 y_{n-2} + b_3 y_{n-3}) mod m_2
;;
;; with coefficients a_1 = 0, a_2 = 63308, a_3 = -183326, b_1 = 86098,
;; b_2 = 0, b_3 = -539608, and moduli m_1 = 2^31 - 1 = 2147483647 and
;; m_2 = 2145483479.
;;
;; The period of this generator is 2^205 (about 10^61). It uses 6
;; words of state per generator. For more information see,
;;
;; * P. L'Ecuyer, "Combined Multiple Recursive Random Number
;; Generators," Operations Research, 44, 5 (1996), 816--822.
;;

(define (make-combined-multiple-recursive-generator)
  (let ((x-1 151) (x-2 2711) (x-3 6719)
        (y-1 7333) (y-2 2447) (y-3 433)
        (a1 0) (a2 63308) (a3 -183326)
        (b1 86098) (b2 0) (b3 -539608)
        (m1 2147483647) (m2 2145483479))
    (lambda ()
      (let ((x (modulo (+ (* a1 x-1) (* a2 x-2) (* a3 x-3)) m1))
            (y (modulo (+ (* b1 y-1) (* b2 y-2) (* b3 y-3)) m2)))
        (set! x-3 x-2)
        (set! x-2 x-1)
        (set! x-1 x)
        (set! y-3 y-2)
        (set! y-2 y-1)
        (set! y-1 y)
        (modulo (- x y) m1)))))

;; -----------------------------------------------------------------
;; tests
(use gauche.time)

(define cmr-max (- (expt 2 31) 1))
(define (rng/float)
  (/ (rng) cmr-max))

(define rng (make-combined-multiple-recursive-generator))

(define (test-rng rng)
  (dotimes (n 3)
    (let ((timer (make <real-time-counter>))
          (t 100000))
      (with-time-counter timer (dotimes (i t) (rng)))
      (format #t "t=~10d, ~d rng/s\n"
              t (/ t (time-counter-value timer))))))

(test-rng rng)