Artifact
af55930880b3b7452f9424331a05c4c77d4d95c4:
#!r6rs
(import (rnrs) (rnrs r5rs) (surfage s27 random-bits))
; CONFIDENCE TESTS FOR SRFI-27 "Sources of Random Bits"
; =====================================================
;
; Sebastian.Egner@philips.com, 2002.
;
; This file contains a small collection of checks for the
; implementation of SRFI-27. It is not meant to be complete
; or to test the actual properties of the underlying generator.
; It is merely meant to run the code and to check some of the
; assumptions made by specification. There is an interface to
; G. Marsaglia's DIEHARD battery of tests for random number
; generators, though.
; History of this file:
; SE, 19-Mar-2002: initial version, based on earlier tests
; SE, 22-Mar-2002: adapted to new procedure names
; SE, 25-Mar-2002: more descriptive output
; SE, 04-Apr-2002: some quick timings; check up
; (check expr)
; evals expr and issues an error if it is not #t.
#;
(define (check expr)
(if (not (eq? (eval expr (interaction-environment)) #t))
(error "check failed" expr)))
; Basic Tests of the Interface
; ============================
(define (my-random-integer n)
(let ((x (random-integer n)))
(if (<= 0 x (- n 1))
x
(error "(random-integer n) returned illegal value" x))))
(define (my-random-real)
(let ((x (random-real)))
(if (< 0 x 1)
x
(error "(random-real) returned illegal value" x))))
(define (check-basics-1)
; generate increasingly large numbers
(display "; generating large numbers [bits]: ")
(do ((k 0 (+ k 1))
(n 1 (* n 2)))
((> k 1024))
(display k)
(display " ")
(my-random-integer n))
(display "ok")
(newline)
; generate some reals
(display "; generating reals [1000 times]: ")
(do ((k 0 (+ k 1))
(x (my-random-real) (+ x (my-random-real))))
((= k 1000)
x))
(display "ok")
(newline)
; get/set the state
(display "; get/set state: ")
(let* ((state1 (random-source-state-ref default-random-source))
(x1 (my-random-integer (expt 2 32)))
(state2 (random-source-state-ref default-random-source))
(x2 (my-random-integer (expt 2 32))))
(random-source-state-set! default-random-source state1)
(let ((y1 (my-random-integer (expt 2 32))))
(if (not (= x1 y1))
(error "state get/set doesn't work" x1 y1 state1)))
(random-source-state-set! default-random-source state2)
(let ((y2 (my-random-integer (expt 2 32))))
(if (not (= x2 y2))
(error "state get/set doesn't work" x2 y2 state2))))
(display "ok")
(newline)
; randomize!
(display "; randomize!: ")
(let* ((state1 (random-source-state-ref default-random-source))
(x1 (my-random-integer (expt 2 32))))
(random-source-state-set! default-random-source state1)
(random-source-randomize! default-random-source)
(let ((y1 (my-random-integer (expt 2 32))))
(if (= x1 y1)
(error "random-source-randomize! didn't work" x1 state1))))
(display "ok")
(newline)
; pseudo-randomize!
(display "; pseudo-randomize!: ")
(let* ((state1 (random-source-state-ref default-random-source))
(x1 (my-random-integer (expt 2 32))))
(random-source-state-set! default-random-source state1)
(random-source-pseudo-randomize! default-random-source 0 1)
(let ((y1 (my-random-integer (expt 2 32))))
(if (= x1 y1)
(error "random-source-pseudo-randomize! didn't work" x1 state1)))
(random-source-state-set! default-random-source state1)
(random-source-pseudo-randomize! default-random-source 1 0)
(let ((y1 (my-random-integer (expt 2 32))))
(if (= x1 y1)
(error "random-source-pseudo-randomize! didn't work" x1 state1))))
(display "ok")
(newline)
(newline))
; Testing the MRG32k3a Generator (if implemented)
; ===============================================
; (check-mrg32k3a)
; tests if the underlying generator is the MRG32k3a generator
; as implemented in the reference implementation. This function
; is useful to check whether the reference implementation computes
; the right numbers.
(define (check-mrg32k3a)
; check if the initial state is A^16 * (1 0 0 1 0 0)
(display "; check A^16 * (1 0 0 1 0 0)")
(let* ((s (make-random-source))
(state1 (random-source-state-ref s))
(rand (random-source-make-reals s)))
(random-source-state-set! s '(lecuyer-mrg32k3a 1 0 0 1 0 0))
(do ((k 0 (+ k 1)))
((= k 16)
(let ((state2 (random-source-state-ref s)))
(if (not (equal? state1 state2))
(error "16-th state after (1 0 0 1 0 0) is wrong"))))
(rand)))
(display "ok")
(newline)
; check if pseudo-randomize! advances properly
(display "; checking (random-source-pseudo-randomize! s 1 2)")
(let ((s (make-random-source)))
(random-source-pseudo-randomize! s 1 2)
(if (not (equal? (random-source-state-ref s)
'(lecuyer-mrg32k3a
1250826159
3004357423
431373563
3322526864
623307378
2983662421)))
(error "pseudo-randomize! gives wrong result")))
(display "ok")
(newline)
; run the check published by Pierre L'Ecuyer:
; Note that the reference implementation deals slightly different
; with reals mapping m1-1 into 1-1/(m1+1) and not into 0 as in
; L'Ecuyer's original proposal. However, for the first 10^7 reals
; that makes no difference as m1-1 is not generated.
(display "; checking (random-source-pseudo-randomize! s 1 2)...")
(let* ((x 0.0)
(s (make-random-source))
(rand (random-source-make-reals s)))
(random-source-state-set!
s
'(lecuyer-mrg32k3a 12345 12345 12345 12345 12345 12345))
(do ((k 0 (+ k 1)))
((= k 10000000)
(if (not (< (abs (- x 5001090.95)) 0.01))
(error "bad sum over 10^7 reals" x)))
(set! x (+ x (rand)))))
(display "ok")
(newline))
; Writing Data to DIEHARD
; =======================
; (write-diehard filename s bytes-per-call calls)
; creates a binary file to which bytes-per-call * calls bytes are
; written. The bytes are obtained from the random source s using
; the range n = (expt 256 bytes-per-call).
; The intention of write-diehard is to give implementors a
; '15 min.'-way of running their favourite random number generator
; through a pretty tough testsuite.
;
; try: For the reference implementation, the call
;
; (write-diehard "outfile" (make-random-source) 4 2867200)
;
; should create a file that looks as follows (od -A x -t x1 outfile):
;
; 0000000 92 bb 7e db 1b 14 f6 bb bb 54 a1 55 c2 3e cd ca
; 0000010 23 01 20 35 06 47 65 b0 52 4c b8 c0 21 48 af 67
; 0000020 63 a9 8c 78 50 73 29 08 62 d1 22 7f a6 89 96 77
; 0000030 98 28 65 2d 2d 8b f9 52 41 be 8e 3f c5 84 0f ca
; 0000040 c0 fa 03 d6 f0 65 9d 3a 9b ab 6f fe d1 aa 5f 92
; 0000050 0f ea f6 3b 78 b9 fe ad 63 5e 49 f1 9d c9 8e 2f
; 0000060 53 a9 5d 32 d4 20 51 1d 1c 2e 82 f0 8b 26 40 c0
; ...total length is 11468800 bytes.
;
; The message digest is md5sum = 4df554f56cb5ed251bd04b0d50767443.
;
; try: For the reference implementation, the call
;
; (write-diehard "outfile" (make-random-source) 3 3822934)
;
; should create a file that looks as follows (od -A x -t x1 outfile):
;
; 000000 bb 7e db 30 a3 49 14 f6 bb d0 f2 d0 54 a1 55 8b
; 000010 8c 03 3e cd ca a3 88 1d 01 20 35 e8 50 c8 47 65
; 000020 b0 e7 d9 28 4c b8 c0 f2 82 35 48 af 67 42 3e 8a
; 000030 a9 8c 78 12 ef b6 73 29 08 ff e9 71 d1 22 7f 52
; 000040 b8 f0 89 96 77 dc 71 86 28 65 2d c2 82 fc 8b f9
; 000050 52 d7 23 2a be 8e 3f 61 a8 99 84 0f ca 44 83 65
; 000060 fa 03 d6 c2 11 c0 65 9d 3a c2 7a dd ab 6f fe 1c
; ...total length is 11468802 bytes.
;
; The message digest is md5sum = 750ac219ff40c50bb2d04ff5eff9b24c.
(define (write-diehard filename s bytes-per-call calls)
(let ((port (open-output-file filename))
(rand (random-source-make-integers s))
(n (expt 256 bytes-per-call)))
(do ((i 0 (+ i 1)))
((= i calls)
(close-output-port port))
(let ((x (rand n)))
(do ((k 0 (+ k 1))) ((= k bytes-per-call))
(put-u8 port (modulo x 256))
(set! x (quotient x 256)))))))
; run some tests
(check-basics-1)
(display "passed (check-basics-1)")
(newline)
(check-mrg32k3a)
(display "passed (check-mrg32k3a)")
(newline)
; (display "Generating diehard1 with expected MD5=4df554f56cb5ed251bd04b0d50767443\n")
; (write-diehard "diehard1" (make-random-source) 4 2867200)
;(display "Generating diehard2 with expected MD5=750ac219ff40c50bb2d04ff5eff9b24c\n")
; (display "Generating diehard2 with expected MD5=9c4cb1f6251efa301a98f226a76de5b9")
; (write-diehard "diehard2" (make-random-source) 3 3822934)