; SRFI-27: Random bits
; 2004-01-11 / lth
;
; $Id$
;
; Larceny version, based on Sebastian Egner's fast all-Scheme
; reference implementation for Gambit-C 3.0, see copyrights and
; credits below.
;
; See also the temporary Larceny patch at the end.

;;; Copyright (C) Sebastian Egner (2002). All Rights Reserved.
;;;
;;; Permission is hereby granted, free of charge, to any person
;;; obtaining a copy of this software and associated documentation
;;; files (the Software), to deal in the Software without restriction,
;;; including without limitation the rights to use, copy, modify,
;;; merge, publish, distribute, sublicense, and/or sell copies of the
;;; Software, and to permit persons to whom the Software is furnished
;;; to do so, subject to the following conditions:
;;; 
;;; The above copyright notice and this permission notice shall be
;;; included in all copies or substantial portions of the Software.
;;; 
;;; THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND,
;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
;;; MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
;;; OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
;;; DEALINGS IN THE SOFTWARE.


(library (srfi :27 random-bits)

  (export random-integer random-real default-random-source
          make-random-source random-source?
          random-source-state-ref random-source-state-set!
          random-source-randomize! random-source-pseudo-randomize!
          random-source-make-integers random-source-make-reals)

  (import (rnrs base)
          (rnrs control)
         ;(only (rnrs arithmetic fixnums) greatest-fixnum)
          (rnrs arithmetic fixnums)
          (rnrs arithmetic flonums)
          (primitives
           bytevector-ref
           fxlogand fxlogior fxlsh fxrshl fxrsha)    ; FIXME: workaround
          (rnrs bytevectors)                         ; FIXME: workaround
          (rnrs r5rs)
          (srfi :9 records)
         ;(primitives r5rs:require current-utc-time) ; FIXME: see below
          (primitives current-seconds
                      memstats
                      memstats-elapsed-time
                      memstats-gc-total-elapsed-time))

; current-utc-time is defined by the time library,
; which uses Larceny's FFI, which doesn't seem to work anymore
; if Larceny was built on a Macintosh using a 64-bit MacOSX.
;
; Until that's fixed, it's safer to define a fake version of
; current-utc-time.
; 
;(define ignored (r5rs:require 'time))           ; CURRENT-UTC-TIME

(define (current-utc-time)
  (let ((stats (memstats)))
    (values (current-seconds)
            (remainder (+ (* 1000 (memstats-gc-total-elapsed-time stats))
                          (memstats-elapsed-time stats))
                       1000000))))

(define most-positive-fixnum greatest-fixnum)

; MODULE DEFINITION FOR SRFI-27
; =============================
; 
; Sebastian.Egner@philips.com, Mar-2002, in Gambit 3.0
;
; This file contains the top-level definition for the 54-bit flonum
; implementation of SRFI-27 for the Gambit 3.0 system.
;
; 1. The core generator is implemented in 'mrg32k3a-c.scm'.
; 2. The generic parts of the interface are in 'mrg32k3a.scm'.
; 3. The non-generic parts (record type, time, error) are here.
;
; loading the module:
;   (load "srfi-27-c.scm")
;
; history of this file:
;   SE, 25-Mar-2002: initial version
;   SE, 10-Apr-2002: incorporated Brad Lucier's optimizations

; (define-structure :random-source
;   state-ref
;   state-set!
;   randomize!
;   pseudo-randomize!
;   make-integers
;   make-reals)

(define-record-type :random-source 
  (:random-source-make a0 a1 a2 a3 a4 a5)
  :random-source?
  (a0 :random-source-state-ref)
  (a1 :random-source-state-set!)
  (a2 :random-source-randomize!)
  (a3 :random-source-pseudo-randomize!)
  (a4 :random-source-make-integers)
  (a5 :random-source-make-reals))

(define (:random-source-current-time)
  (let-values (((secs usecs) (current-utc-time)))
    (+ secs usecs)))


; 54-BIT FLONUM IMPLEMENTATION OF THE "MRG32K3A"-GENERATOR
; ========================================================
;
; Sebastian.Egner@philips.com, 2002.
;
; This code is a floating point implementation in Scheme of Pierre L'Ecuyer's 
; MRG32k3a generator. The implementation runs under the Gambit 3.0 Scheme
; system. It makes use of Gambit's f64vector. The code is based on Brad
; Lucier's code for this generator. Please refer to Brad's original posting
; in the SRFI-27 discussion archive and to his posting of an optimization.
;
; compliance:
;   Scheme R5RS with flonums covering at least {-2^53..2^53-1} exactly.
;   In addition: declare f64vector ##fixnum?
;
; loading this file into Gambit 3.0:
;   (load "mrg32k3a-2.scm")
;
; history of this file:
;   SE, 18-Mar-2002: initial version
;   SE, 25-Mar-2002: adapted to new interface
;   SE, 04-Apr-2002: debugged; simplified
;   SE, 10-Apr-2002: incorporated Brad Lucier's optimizations

; the actual generator

(define (mrg32k3a-random-m1 state) ; from Brad Lucier
  (let* ((x10 (- (* 1403580.0 (vector-ref state 1)) 
                 (* 810728.0 (vector-ref state 2))))
         (y10 (- x10 (* (floor (/ x10 4294967087.0)) 4294967087.0)))
         (x20 (- (* 527612.0 (vector-ref state 3)) 
                 (* 1370589.0 (vector-ref state 5))))
         (y20 (- x20 (* (floor (/ x20 4294944443.0)) 4294944443.0)))
         (dx (- y10 y20))
         (dy (- dx (* (floor (/ dx 4294967087.0)) 4294967087.0))))
    (vector-set! state 5 (vector-ref state 4))
    (vector-set! state 4 (vector-ref state 3))
    (vector-set! state 3 y20)
    (vector-set! state 2 (vector-ref state 1))
    (vector-set! state 1 (vector-ref state 0))
    (vector-set! state 0 y10)
    dy))

; interface to 'random.scm'

(define (mrg32k3a-pack-state unpacked-state)
  (list->vector (map exact->inexact (vector->list unpacked-state))))

(define (mrg32k3a-unpack-state state)
  (list->vector (map inexact->exact (vector->list state))))

(define (mrg32k3a-random-range)
  (most-positive-fixnum))

;;; Larceny: commented out in favor of patch at end.

(define (mrg32k3a-random-integer state range) ; from Brad Lucier
  (let* ((n (exact->inexact range))
         (q (floor (/ 4294967087.0 n)))
         (qn (* q n)))
    (do ((x (mrg32k3a-random-m1 state) (mrg32k3a-random-m1 state)))
        ((< x qn) 
         (inexact->exact (floor (/ x q)))))))

(define (mrg32k3a-random-real state)
  (* 0.0000000002328306549295728 (+ 1.0 (mrg32k3a-random-m1 state))))


; GENERIC PART OF MRG32k3a-GENERATOR FOR SRFI-27
; ==============================================
;
; Sebastian.Egner@philips.com, 2002.
;
; This is the generic R5RS-part of the implementation of the MRG32k3a
; generator to be used in SRFI-27. It is based on a separate implementation
; of the core generator (presumably in native code) and on code to
; provide essential functionality not available in R5RS (see below).
;
; compliance:
;   Scheme R5RS with integer covering at least {-2^53..2^53-1}.
;   In addition,
;     SRFI-23: error
;
; history of this file:
;   SE, 22-Mar-2002: refactored from earlier versions
;   SE, 25-Mar-2002: pack/unpack need not allocate
;   SE, 27-Mar-2002: changed interface to core generator
;   SE, 10-Apr-2002: updated spec of mrg32k3a-random-integer

; Generator
; =========
;
; Pierre L'Ecuyer's MRG32k3a generator is a Combined Multiple Recursive 
; Generator. It produces the sequence {(x[1,n] - x[2,n]) mod m1 : n} 
; defined by the two recursive generators
;
;   x[1,n] = (               a12 x[1,n-2] + a13 x[1,n-3]) mod m1,
;   x[2,n] = (a21 x[2,n-1] +                a23 x[2,n-3]) mod m2,
;
; where the constants are
;   m1       = 4294967087 = 2^32 - 209    modulus of 1st component
;   m2       = 4294944443 = 2^32 - 22853  modulus of 2nd component
;   a12      =  1403580                   recursion coefficients
;   a13      =  -810728
;   a21      =   527612
;   a23      = -1370589
;
; The generator passes all tests of G. Marsaglia's Diehard testsuite.
; Its period is (m1^3 - 1)(m2^3 - 1)/2 which is nearly 2^191.
; L'Ecuyer reports: "This generator is well-behaved in all dimensions
; up to at least 45: ..." [with respect to the spectral test, SE].
;
; The period is maximal for all values of the seed as long as the
; state of both recursive generators is not entirely zero.
;
; As the successor state is a linear combination of previous
; states, it is possible to advance the generator by more than one
; iteration by applying a linear transformation. The following
; publication provides detailed information on how to do that:
;
;    [1] P. L'Ecuyer, R. Simard, E. J. Chen, W. D. Kelton:
;        An Object-Oriented Random-Number Package With Many Long 
;        Streams and Substreams. 2001.
;        To appear in Operations Research.
;
; Arithmetics
; ===========
;
; The MRG32k3a generator produces values in {0..2^32-209-1}. All 
; subexpressions of the actual generator fit into {-2^53..2^53-1}. 
; The code below assumes that Scheme's "integer" covers this range.
; In addition, it is assumed that floating point literals can be
; read and there is some arithmetics with inexact numbers.
;
; However, for advancing the state of the generator by more than
; one step at a time, the full range {0..2^32-209-1} is needed.


; Required: Backbone Generator
; ============================
;
; At this point in the code, the following procedures are assumed
; to be defined to execute the core generator:
;
;   (mrg32k3a-pack-state unpacked-state) -> packed-state
;   (mrg32k3a-unpack-state packed-state) -> unpacked-state
;      pack/unpack a state of the generator. The core generator works
;      on packed states, passed as an explicit argument, only. This
;      allows native code implementations to store their state in a
;      suitable form. Unpacked states are #(x10 x11 x12 x20 x21 x22) 
;      with integer x_ij. Pack/unpack need not allocate new objects
;      in case packed and unpacked states are identical.
;
;   (mrg32k3a-random-range) -> m-max
;   (mrg32k3a-random-integer packed-state range) -> x in {0..range-1}
;      advance the state of the generator and return the next random
;      range-limited integer. 
;        Note that the state is not necessarily advanced by just one 
;      step because we use the rejection method to avoid any problems 
;      with distribution anomalies.
;        The range argument must be an exact integer in {1..m-max}.
;      It can be assumed that range is a fixnum if the Scheme system
;      has such a number representation.
;
;   (mrg32k3a-random-real packed-state) -> x in (0,1)
;      advance the state of the generator and return the next random
;      real number between zero and one (both excluded). The type of
;      the result should be a flonum if possible.

; Required: Record Data Type
; ==========================
;
; At this point in the code, the following procedures are assumed
; to be defined to create and access a new record data type:
;
;   (:random-source-make a0 a1 a2 a3 a4 a5) -> s
;     constructs a new random source object s consisting of the 
;     objects a0 .. a5 in this order.
;
;   (:random-source? obj) -> bool
;     tests if a Scheme object is a :random-source.
;
;   (:random-source-state-ref         s) -> a0
;   (:random-source-state-set!        s) -> a1
;   (:random-source-randomize!        s) -> a2
;   (:random-source-pseudo-randomize! s) -> a3
;   (:random-source-make-integers     s) -> a4
;   (:random-source-make-reals        s) -> a5
;     retrieve the values in the fields of the object s.

; Required: Current Time as an Integer
; ====================================
;
; At this point in the code, the following procedure is assumed
; to be defined to obtain a value that is likely to be different
; for each invokation of the Scheme system:
;
;   (:random-source-current-time) -> x
;     an integer that depends on the system clock. It is desired
;     that the integer changes as fast as possible.


; Accessing the State
; ===================

(define (mrg32k3a-state-ref packed-state)
  (cons 'lecuyer-mrg32k3a 
        (vector->list (mrg32k3a-unpack-state packed-state))))

(define (mrg32k3a-state-set external-state)

  (define (check-value x m)
    (if (and (integer? x)
             (exact? x)
             (<= 0 x (- m 1)))
        #t
        (error "illegal value" x)))

  (if (and (list? external-state)
           (= (length external-state) 7)
           (eq? (car external-state) 'lecuyer-mrg32k3a))
      (let ((s (cdr external-state)))
        (check-value (list-ref s 0) mrg32k3a-m1)
        (check-value (list-ref s 1) mrg32k3a-m1)
        (check-value (list-ref s 2) mrg32k3a-m1)
        (check-value (list-ref s 3) mrg32k3a-m2)
        (check-value (list-ref s 4) mrg32k3a-m2)
        (check-value (list-ref s 5) mrg32k3a-m2)
        (if (or (zero? (+ (list-ref s 0) (list-ref s 1) (list-ref s 2)))
                (zero? (+ (list-ref s 3) (list-ref s 4) (list-ref s 5))))
            (error "illegal degenerate state" external-state))
        (mrg32k3a-pack-state (list->vector s)))
      (error "malformed state" external-state)))


; Pseudo-Randomization
; ====================
;
; Reference [1] above shows how to obtain many long streams and 
; substream from the backbone generator.
;
; The idea is that the generator is a linear operation on the state.
; Hence, we can express this operation as a 3x3-matrix acting on the
; three most recent states. Raising the matrix to the k-th power, we
; obtain the operation to advance the state by k steps at once. The
; virtual streams and substreams are now simply parts of the entire
; periodic sequence (which has period around 2^191).
;
; For the implementation it is necessary to compute with matrices in
; the ring (Z/(m1*m1)*Z)^(3x3). By the Chinese-Remainder Theorem, this
; is isomorphic to ((Z/m1*Z) x (Z/m2*Z))^(3x3). We represent such a pair
; of matrices 
;   [ [[x00 x01 x02],
;      [x10 x11 x12],
;      [x20 x21 x22]], mod m1
;     [[y00 y01 y02],
;      [y10 y11 y12],
;      [y20 y21 y22]]  mod m2]
; as a vector of length 18 of the integers as writen above:
;   #(x00 x01 x02 x10 x11 x12 x20 x21 x22
;     y00 y01 y02 y10 y11 y12 y20 y21 y22)
;
; As the implementation should only use the range {-2^53..2^53-1}, the
; fundamental operation (x*y) mod m, where x, y, m are nearly 2^32, 
; is computed by breaking up x and y as x = x1*w + x0 and y = y1*w + y0 
; where w = 2^16. In this case, all operations fit the range because 
; w^2 mod m is a small number. If proper multiprecision integers are
; available this is not necessary, but pseudo-randomize! is an expected
; to be called only occasionally so we do not provide this implementation.

(define mrg32k3a-m1 4294967087) ; modulus of component 1
(define mrg32k3a-m2 4294944443) ; modulus of component 2

(define mrg32k3a-initial-state ; 0 3 6 9 12 15 of A^16, see below
  '#( 1062452522
      2961816100 
       342112271 
      2854655037 
      3321940838 
      3542344109))

(define mrg32k3a-generators #f) ; computed when needed

(define (mrg32k3a-pseudo-randomize-state i j)

  (define (product A B) ; A*B in ((Z/m1*Z) x (Z/m2*Z))^(3x3)

    (define w      65536) ; wordsize to split {0..2^32-1}
    (define w-sqr1 209)   ; w^2 mod m1
    (define w-sqr2 22853) ; w^2 mod m2

    (define (lc i0 i1 i2 j0 j1 j2 m w-sqr) ; linear combination
      (let ((a0h (quotient (vector-ref A i0) w))
            (a0l (modulo   (vector-ref A i0) w))
            (a1h (quotient (vector-ref A i1) w))
            (a1l (modulo   (vector-ref A i1) w))
            (a2h (quotient (vector-ref A i2) w))
            (a2l (modulo   (vector-ref A i2) w))
            (b0h (quotient (vector-ref B j0) w))
            (b0l (modulo   (vector-ref B j0) w))
            (b1h (quotient (vector-ref B j1) w))
            (b1l (modulo   (vector-ref B j1) w))
            (b2h (quotient (vector-ref B j2) w))
            (b2l (modulo   (vector-ref B j2) w)))
        (modulo
         (+ (* (+ (* a0h b0h) 
                  (* a1h b1h) 
                  (* a2h b2h)) 
               w-sqr)
            (* (+ (* a0h b0l) 
                  (* a0l b0h)
                  (* a1h b1l) 
                  (* a1l b1h)
                  (* a2h b2l) 
                  (* a2l b2h))
               w)
            (* a0l b0l)
            (* a1l b1l)
            (* a2l b2l))
         m)))
    
    (vector
     (lc  0  1  2   0  3  6  mrg32k3a-m1 w-sqr1) ; (A*B)_00 mod m1
     (lc  0  1  2   1  4  7  mrg32k3a-m1 w-sqr1) ; (A*B)_01
     (lc  0  1  2   2  5  8  mrg32k3a-m1 w-sqr1)
     (lc  3  4  5   0  3  6  mrg32k3a-m1 w-sqr1) ; (A*B)_10
     (lc  3  4  5   1  4  7  mrg32k3a-m1 w-sqr1)
     (lc  3  4  5   2  5  8  mrg32k3a-m1 w-sqr1)
     (lc  6  7  8   0  3  6  mrg32k3a-m1 w-sqr1)
     (lc  6  7  8   1  4  7  mrg32k3a-m1 w-sqr1)
     (lc  6  7  8   2  5  8  mrg32k3a-m1 w-sqr1)
     (lc  9 10 11   9 12 15  mrg32k3a-m2 w-sqr2) ; (A*B)_00 mod m2
     (lc  9 10 11  10 13 16  mrg32k3a-m2 w-sqr2)
     (lc  9 10 11  11 14 17  mrg32k3a-m2 w-sqr2)
     (lc 12 13 14   9 12 15  mrg32k3a-m2 w-sqr2)
     (lc 12 13 14  10 13 16  mrg32k3a-m2 w-sqr2)
     (lc 12 13 14  11 14 17  mrg32k3a-m2 w-sqr2)
     (lc 15 16 17   9 12 15  mrg32k3a-m2 w-sqr2)
     (lc 15 16 17  10 13 16  mrg32k3a-m2 w-sqr2)
     (lc 15 16 17  11 14 17  mrg32k3a-m2 w-sqr2)))

  (define (power A e) ; A^e
    (cond
     ((zero? e)
      '#(1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1))
     ((= e 1)
      A)
     ((even? e)
      (power (product A A) (quotient e 2)))
     (else
      (product (power A (- e 1)) A))))

  (define (power-power A b) ; A^(2^b)
    (if (zero? b)
        A
        (power-power (product A A) (- b 1))))

  (define A                        ; the MRG32k3a recursion
    '#(     0 1403580 4294156359
            1       0          0
            0       1          0
       527612       0 4293573854
            1       0          0
            0       1          0))

  ; check arguments
  (if (not (and (integer? i) 
                (exact? i)
                (integer? j)
                (exact? j)))
      (error "i j must be exact integer" i j))

  ; precompute A^(2^127) and A^(2^76) only once

  (if (not mrg32k3a-generators)
      (set! mrg32k3a-generators
            (list (power-power A 127)
                  (power-power A  76)
                  (power A 16))))

  ; compute M = A^(16 + i*2^127 + j*2^76)
  (let ((M (product 
            (list-ref mrg32k3a-generators 2)
            (product
             (power (list-ref mrg32k3a-generators 0)
                    (modulo i (expt 2 28)))
             (power (list-ref mrg32k3a-generators 1) 
                    (modulo j (expt 2 28)))))))
    (mrg32k3a-pack-state
     (vector
      (vector-ref M 0)
      (vector-ref M 3)
      (vector-ref M 6)
      (vector-ref M 9)
      (vector-ref M 12)
      (vector-ref M 15)))))

; True Randomization
; ==================
;
; The value obtained from the system time is feed into a very
; simple pseudo random number generator. This in turn is used
; to obtain numbers to randomize the state of the MRG32k3a
; generator, avoiding period degeneration.

(define (mrg32k3a-randomize-state state)
  ;; G. Marsaglia's simple 16-bit generator with carry
  (let* ((m 65536)
         (x (modulo (:random-source-current-time) m)))
    (define (random-m)
      (let ((y (modulo x m)))
        (set! x (+ (* 30903 y) (quotient x m)))
        y))
    (define (random n)                        ; m < n < m^2
      (modulo (+ (* (random-m) m) (random-m)) n))

                                        ; modify the state
    (let ((m1 mrg32k3a-m1)
          (m2 mrg32k3a-m2)
          (s (mrg32k3a-unpack-state state)))
      (mrg32k3a-pack-state
       (vector
        (+ 1 (modulo (+ (vector-ref s 0) (random (- m1 1))) (- m1 1)))
        (modulo (+ (vector-ref s 1) (random m1)) m1)
        (modulo (+ (vector-ref s 2) (random m1)) m1)
        (+ 1 (modulo (+ (vector-ref s 3) (random (- m2 1))) (- m2 1)))
        (modulo (+ (vector-ref s 4) (random m2)) m2)
        (modulo (+ (vector-ref s 5) (random m2)) m2))))))


; Large Integers
; ==============
;
; To produce large integer random deviates, for n > m-max, we first 
; construct large random numbers in the range {0..m-max^k-1} for some 
; k such that m-max^k >= n and then use the rejection method to choose
; uniformly from the range {0..n-1}.

(define mrg32k3a-m-max
  (mrg32k3a-random-range))

(define (mrg32k3a-random-power state k) ; n = m-max^k, k >= 1
  (if (= k 1)
      (mrg32k3a-random-integer state mrg32k3a-m-max)
      (+ (* (mrg32k3a-random-power state (- k 1)) mrg32k3a-m-max)
         (mrg32k3a-random-integer state mrg32k3a-m-max))))

(define (mrg32k3a-random-large state n) ; n > m-max
  (do ((k 2 (+ k 1))
       (mk (* mrg32k3a-m-max mrg32k3a-m-max) (* mk mrg32k3a-m-max)))
      ((>= mk n)
       (let* ((mk-by-n (quotient mk n))
              (a (* mk-by-n n)))
         (do ((x (mrg32k3a-random-power state k)
                 (mrg32k3a-random-power state k)))
             ((< x a) (quotient x mk-by-n)))))))


; Multiple Precision Reals
; ========================
;
; To produce multiple precision reals we produce a large integer value
; and convert it into a real value. This value is then normalized.
; The precision goal is unit <= 1/(m^k + 1), or 1/unit - 1 <= m^k.
; If you know more about the floating point number types of the
; Scheme system, this can be improved.

(define (mrg32k3a-random-real-mp state unit)
  (do ((k 1 (+ k 1))
       (u (- (/ 1 unit) 1) (/ u mrg32k3a-m1)))
      ((<= u 1)
       (/ (exact->inexact (+ (mrg32k3a-random-power state k) 1))
          (exact->inexact (+ (expt mrg32k3a-m-max k) 1))))))


; Provide the Interface as Specified in the SRFI
; ==============================================
;
; An object of type random-source is a record containing the procedures
; as components. The actual state of the generator is stored in the
; binding-time environment of make-random-source.

(define (make-random-source)
  (let ((state (mrg32k3a-pack-state ; make a new copy
                (list->vector (vector->list mrg32k3a-initial-state)))))
    (:random-source-make
     (lambda ()
       (mrg32k3a-state-ref state))
     (lambda (new-state)
       (set! state (mrg32k3a-state-set new-state)))
     (lambda ()
       (set! state (mrg32k3a-randomize-state state)))
     (lambda (i j)
       (set! state (mrg32k3a-pseudo-randomize-state i j)))
     (lambda ()
       (lambda (n)
         (cond
          ((not (and (integer? n) (exact? n) (positive? n)))
           (error "range must be exact positive integer" n))           
          ((<= n mrg32k3a-m-max)
           (mrg32k3a-random-integer state n))
          (else
           (mrg32k3a-random-large state n)))))
     (lambda args
       (cond
        ((null? args)
         (lambda () 
           (mrg32k3a-random-real state)))
        ((null? (cdr args))
         (let ((unit (car args)))
           (cond
            ((not (and (real? unit) (< 0 unit 1)))
             (error "unit must be real in (0,1)" unit))
            ((<= (- (/ 1 unit) 1) mrg32k3a-m1)
             (lambda () 
               (mrg32k3a-random-real state)))
            (else
             (lambda () 
               (mrg32k3a-random-real-mp state unit))))))
        (else
         (error "illegal arguments" args)))))))

(define random-source? 
  :random-source?)

(define (random-source-state-ref s)
  ((:random-source-state-ref s)))

(define (random-source-state-set! s state)
  ((:random-source-state-set! s) state))

(define (random-source-randomize! s)
  ((:random-source-randomize! s)))

(define (random-source-pseudo-randomize! s i j)
  ((:random-source-pseudo-randomize! s) i j))

; ---

(define (random-source-make-integers s)
  ((:random-source-make-integers s)))

(define (random-source-make-reals s . unit)
  (apply (:random-source-make-reals s) unit))

; ---

(define default-random-source 
  (make-random-source))

(define random-integer
  (random-source-make-integers default-random-source))

(define random-real
  (random-source-make-reals default-random-source))

)

(library (srfi :27)

  (export random-integer random-real default-random-source
          make-random-source random-source?
          random-source-state-ref random-source-state-set!
          random-source-randomize! random-source-pseudo-randomize!
          random-source-make-integers random-source-make-reals)

  (import (srfi :27 random-bits)))

; eof
