110 lines
2.9 KiB
Scheme
110 lines
2.9 KiB
Scheme
; This is a purely functional random generator, based on Lennart
|
|
; Augustsson's generator shipped with Hugs.
|
|
|
|
; Its comment says:
|
|
|
|
; This implementation uses the Portable Combined Generator of L'Ecuyer
|
|
; ["System.Random\#LEcuyer"] for 32-bit computers, transliterated by
|
|
; Lennart Augustsson. It has a period of roughly 2.30584e18.
|
|
|
|
; This makes it not as good as Sebastian Egner's reference
|
|
; implementation of SRFI 27, but much faster for applications that
|
|
; need a pure random-number generator with a `split' operation.
|
|
|
|
(define-record-type :random-generator
|
|
(really-make-random-generator s1 s2)
|
|
random-generator?
|
|
(s1 random-generator-s1)
|
|
(s2 random-generator-s2))
|
|
|
|
(define-record-discloser :random-generator
|
|
(lambda (r)
|
|
(list 'random-generator
|
|
(random-generator-s1 r)
|
|
(random-generator-s2 r))))
|
|
|
|
(define min-bound (- (expt 2 31)))
|
|
(define max-bound (- (expt 2 31) 1))
|
|
(define int-range (- max-bound min-bound))
|
|
|
|
(define (make-random-generator s)
|
|
(if (negative? s)
|
|
(make-random-generator (- s))
|
|
(let ((q (quotient s 2147483562))
|
|
(s1 (remainder s 2147483562)))
|
|
(let ((s2 (remainder q 2147483398)))
|
|
(really-make-random-generator (+ 1 s1) (+ 1 s2))))))
|
|
|
|
(define (random-generator-next rg)
|
|
(let ((s1 (random-generator-s1 rg))
|
|
(s2 (random-generator-s2 rg)))
|
|
|
|
(let ((k (quotient s1 53668))
|
|
(k* (quotient s2 52774)))
|
|
(let ((s1* (- (* 40014 (- s1 (* k 53668)))
|
|
(* k 12211)))
|
|
(s2* (- (* 40692 (- s2 (* k* 52774)))
|
|
(* k* 3791))))
|
|
(let ((s1** (if (negative? s1*)
|
|
(+ s1* 2147483563)
|
|
s1*))
|
|
(s2** (if (negative? s2*)
|
|
(+ s2* 2147483399)
|
|
s2*)))
|
|
(let* ((z (- s1** s2**))
|
|
(z* (if (< z 1)
|
|
(+ z 2147483562)
|
|
z)))
|
|
(values z* (really-make-random-generator s1** s2**))))))))
|
|
|
|
(define (random-generator-split rg)
|
|
(let ((s1 (random-generator-s1 rg))
|
|
(s2 (random-generator-s2 rg)))
|
|
(let ((new-s1 (if (= s1 2147483562)
|
|
1
|
|
(+ s1 1)))
|
|
(new-s2 (if (= s2 1)
|
|
2147483398
|
|
(- s2 1))))
|
|
(call-with-values
|
|
(lambda ()
|
|
(random-generator-next rg))
|
|
(lambda (_ nrg)
|
|
(values (really-make-random-generator new-s1
|
|
(random-generator-s2 nrg))
|
|
(really-make-random-generator (random-generator-s1 nrg)
|
|
new-s2)))))))
|
|
|
|
|
|
; The intervals are inclusive.
|
|
|
|
(define (random-integer rg low high)
|
|
(let ((b 2147483561)
|
|
(k (+ (- high low) 1)))
|
|
(let loop ((n (ilogbase b k))
|
|
(acc low)
|
|
(rg rg))
|
|
(if (zero? n)
|
|
(values (+ low (modulo acc k))
|
|
rg)
|
|
(call-with-values
|
|
(lambda () (random-generator-next rg))
|
|
(lambda (x rgn)
|
|
(loop (- n 1) (+ x (* acc b)) rgn)))))))
|
|
|
|
(define (random-real rg low high)
|
|
(call-with-values
|
|
(lambda ()
|
|
(random-integer rg min-bound max-bound))
|
|
(lambda (x nrg)
|
|
(let ((scaled-x (+ (/ (+ low high) 2)
|
|
(* (/ (- high low) int-range)
|
|
x))))
|
|
(values scaled-x nrg)))))
|
|
|
|
(define (ilogbase b i)
|
|
(if (< i b)
|
|
1
|
|
(+ 1 (ilogbase b (quotient i b)))))
|
|
|