From bae0952bf7b66a174f23492b95b1183d8aff48fa Mon Sep 17 00:00:00 2001 From: Matthew Flatt Date: Wed, 22 Aug 2007 17:01:20 +0000 Subject: [PATCH] Chongkai's SRFI-27 re-implementation in terms of the native MzScheme random-number procedures svn: r7145 --- collects/srfi/27/random-bits.ss | 693 +++++--------------------------- 1 file changed, 97 insertions(+), 596 deletions(-) diff --git a/collects/srfi/27/random-bits.ss b/collects/srfi/27/random-bits.ss index ff0c473298..a62d614c79 100644 --- a/collects/srfi/27/random-bits.ss +++ b/collects/srfi/27/random-bits.ss @@ -1,598 +1,99 @@ -; MODULE DEFINITION FOR SRFI-27 -; ============================= -; -; Sebastian.Egner@philips.com, Mar-2002, in PLT 204 -; -; This file contains the top-level definition for the 54-bit integer-only -; implementation of SRFI 27 for the PLT 204 DrScheme system. -; -; 1. The core generator is implemented in 'mrg32k3a-a.scm'. -; 2. The generic parts of the interface are in 'mrg32k3a.scm'. -; 3. The non-generic parts (record type, time, error) are here. -; -; load the module with -; (require (lib "random-bits.ss" "srfi")) -; -; history of this file: -; SE, 17-May-2003: initial version - (module random-bits mzscheme - (require (lib "9.ss" "srfi")) - (require (lib "23.ss" "srfi")) - - (provide - 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) - - (define-record-type :random-source - (:random-source-make - state-ref - state-set! - randomize! - pseudo-randomize! - make-integers - make-reals) - :random-source? - (state-ref :random-source-state-ref) - (state-set! :random-source-state-set!) - (randomize! :random-source-randomize!) - (pseudo-randomize! :random-source-pseudo-randomize!) - (make-integers :random-source-make-integers) - (make-reals :random-source-make-reals)) - (define :random-source-current-time - current-milliseconds) - -; implementation begins here - -; 54-BIT INTEGER IMPLEMENTATION OF THE "MRG32K3A"-GENERATOR -; ========================================================= -; -; Sebastian.Egner@philips.com, Mar-2002. -; -; This file is an implementation of Pierre L'Ecuyer's MRG32k3a -; pseudo random number generator. Please refer to 'mrg32k3a.scm' -; for more information. -; -; compliance: -; Scheme R5RS with integers covering at least {-2^53..2^53-1}. -; -; history of this file: -; SE, 18-Mar-2002: initial version -; SE, 22-Mar-2002: comments adjusted, range added -; SE, 25-Mar-2002: pack/unpack just return their argument - -; the actual generator - -(define (mrg32k3a-random-m1 state) - (let ((x11 (vector-ref state 0)) - (x12 (vector-ref state 1)) - (x13 (vector-ref state 2)) - (x21 (vector-ref state 3)) - (x22 (vector-ref state 4)) - (x23 (vector-ref state 5))) - (let ((x10 (modulo (- (* 1403580 x12) (* 810728 x13)) 4294967087)) - (x20 (modulo (- (* 527612 x21) (* 1370589 x23)) 4294944443))) - (vector-set! state 0 x10) - (vector-set! state 1 x11) - (vector-set! state 2 x12) - (vector-set! state 3 x20) - (vector-set! state 4 x21) - (vector-set! state 5 x22) - (modulo (- x10 x20) 4294967087)))) - -; interface to the generic parts of the generator - -(define (mrg32k3a-pack-state unpacked-state) - unpacked-state) - -(define (mrg32k3a-unpack-state state) - state) - -(define (mrg32k3a-random-range) ; m1 - 4294967087) - -(define (mrg32k3a-random-integer state range) ; rejection method - (let* ((q (quotient 4294967087 range)) - (qn (* q range))) - (do ((x (mrg32k3a-random-m1 state) (mrg32k3a-random-m1 state))) - ((< x qn) (quotient x q))))) - -(define (mrg32k3a-random-real state) ; normalization is 1/(m1+1) - (* 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 - (define m 65536) - (define 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)) - - -) ; module ends + (require (lib "contract.ss")) + + (provide random-real + default-random-source + (rename s:make-random-source make-random-source) + random-source? + random-source-state-ref + random-source-state-set! + random-source-randomize! + random-source-pseudo-randomize!) + + (define positive-integer/c + (and/c integer? positive? exact?)) + + (define-struct random-source (generator)) + + (provide/contract + (random-integer (-> positive-integer/c any)) + (random-source-make-integers (-> positive-integer/c any)) + (random-source-make-reals (case-> (-> random-source? any) + (-> random-source? (and/c (>/c 0) (vector (random-source-generator s))) + + (define-syntax random-source-state-set! + (syntax-rules (default-random-source) + ((_ default-random-source state) + (current-pseudo-random-generator (vector->pseudo-random-generator state))) + ((_ s state) + (set-random-source-generator! s + (vector->pseudo-random-generator state))))) + + (define-syntax random-source-randomize! + (syntax-rules (default-random-source) + ((_ default-random-source) + (current-pseudo-random-generator (make-pseudo-random-generator))) + ((_ s) + (set-random-source-generator! s + (make-pseudo-random-generator))))) + + (define-syntax random-source-pseudo-randomize! + (syntax-rules (default-random-source) + ((_ default-random-source . ij) + (random-seed (modulo (equal-hash-code 'ij) 2147483648))) + ((_ s . ij) + (parameterize ((current-pseudo-random-generator + (random-source-generator s))) + (random-seed (modulo (equal-hash-code 'ij) 2147483648)))))) + + (define (random-source-make-integers s) + (lambda (n) + (parameterize ((current-pseudo-random-generator + (random-source-generator s))) + (random-integer n)))) + + (define random-source-make-reals + (case-lambda + ((s) + (lambda () + (parameterize ((current-pseudo-random-generator + (random-source-generator s))) + (random)))) + ((s unit) + (let ((n (- (inexact->exact (floor (/ unit))) 2))) + (lambda () + (parameterize ((current-pseudo-random-generator + (random-source-generator s))) + (* (add1 (random n)) unit))))))) + + (random-seed 0))