use CHICKEN SRFI-27 generators whenever possible

1 files changed, 20 insertions, 320 deletions

diff --git a/194-impl.scm b/194-impl.scm
index 8936d35..ae98046 100644
--- a/194-impl.scm
+++ b/194-impl.scm
@@ -40,18 +40,8 @@
 ;;
 
 (define (make-random-integer-generator low-bound up-bound)
-  (unless (and (integer? low-bound)
-               (exact? low-bound))
-    (error "expected exact integer for lower bound"))
-  (unless (and (integer? up-bound)
-               (exact? up-bound))
-    (error "expected exact integer for upper bound"))
-  (unless (< low-bound up-bound)
-    (error "upper bound should be greater than lower bound"))
-  (let ((rand-int-proc (random-source-make-integers (current-random-source)))
-        (range (- up-bound low-bound)))
-    (lambda ()
-      (+ low-bound (rand-int-proc range)))))
+  (make-uniform-random-integers high: (- up-bound 1)
+                                low: low-bound))
 
 (define (make-random-u1-generator)
   (make-random-integer-generator 0 2))
@@ -142,10 +132,7 @@
                 (r (sqrt (+ (* m t) b))))
           (+ origin (make-polar r phi))))))))
 
-(define (make-random-boolean-generator)
-  (define u1 (make-random-u1-generator))
-  (lambda ()
-    (zero? (u1))))
+(define (make-random-boolean-generator) (make-random-bernoullis))
 
 (define (make-random-char-generator str)
   (when (not (string? str))
@@ -168,16 +155,17 @@
 
 (define PI (* 4 (atan 1.0)))
 
-(define (make-bernoulli-generator p)
-  (unless (real? p)
-    (error "expected p to be real"))
-  (unless (<= 0 p 1)
-    (error "expected 0 <= p <= 1"))
-  (let ((rand-real-proc (random-source-make-reals (current-random-source))))
-   (lambda ()
-     (if (<= (rand-real-proc) p)
-         1
-         0))))
+(define make-bernoulli-generator
+  (case-lambda
+    (() (make-bernoulli-generator 1/2 (current-random-real)))
+    ((p) (make-bernoulli-generator p (current-random-real)))
+    ((p source)
+     (gmap (lambda (x) (if x 1 0))
+           (make-random-bernoullis p: p randoms: source)))))
+
+(define (make-binomial-generator n p)
+  ;; Already in CHICKEN's SRFI-27
+  (make-random-binomials t: n p: p randoms: (current-random-real)))
 
 (define (make-categorical-generator weights-vec)
   (define weight-sum
@@ -203,40 +191,13 @@
            (it newsum
                (+ i 1)))))))
 
-;; Normal distribution (continuous - generates real numbers)
-;; Box-Muller algorithm
-;; NB: We tested Ziggurat method, too,
-;; only to find out Box-Muller is faster about 12% - presumably
-;; the overhead of each ops is larger in Gauche than C/C++, and
-;; so the difference of cost of log or sin from the primitive
-;; addition/multiplication are negligible.
-
-;; NOTE: this implementation is not thread safe
 (define make-normal-generator
   (case-lambda
-    (()
-     (make-normal-generator 0.0 1.0))
-    ((mean)
-     (make-normal-generator mean 1.0))
+    (() (make-normal-generator #i0 #i1))
+    ((mean) (make-normal-generator mean #i1))
     ((mean deviation)
-     (let ((rand-real-proc (random-source-make-reals (current-random-source)))
-           (state #f))
-       (unless (and (real? mean)
-                    (finite? mean))
-         (error "expected mean to be finite real number"))
-       (unless (and (real? deviation)
-                    (finite? deviation)
-                    (> deviation 0))
-         (error "expected deviation to be positive finite real number"))
-       (lambda ()
-         (if state
-             (let ((result state))
-              (set! state #f)
-              result)
-             (let* ((r (sqrt (* -2 (log (rand-real-proc)))))
-                    (theta (* 2 PI (rand-real-proc))))
-               (set! state (+ mean (* deviation r (cos theta))))
-               (+ mean (* deviation r (sin theta))))))))))
+     (make-random-normals mu: mean sigma: deviation
+                          randoms: (current-random-real)))))
 
 (define (make-exponential-generator mean)
   (unless (and (real? mean)
@@ -248,114 +209,10 @@
      (- (* mean (log (rand-real-proc)))))))
 
 (define (make-geometric-generator p)
+  (make-random-geometrics p: p randoms: (current-random-real)))
 
-  (define (log1p x)
-    ;; Adapted from SRFI 144
-    (let ((u (+ 1.0 x)))
-      (cond ((= u 1.0)
-             x) ;; gets sign of zero result correct
-            ((= u x)
-             (log u)) ;; large arguments and infinities
-            (else
-             (* (log u) (/ x (- u 1.0)))))))
-
-  (unless (and (real? p)
-               (> p 0)
-               (<= p 1))
-          (error "expected p to be real number, 0 < p <= 1"))
-  (if (zero? (- p 1.))
-      ;; p is indistinguishable from 1.
-      (lambda () 1)
-      (let ((c (/ (log1p (- p))))
-            (rand-real-proc (random-source-make-reals (current-random-source))))
-        (lambda ()
-          (exact (ceiling (* c (log (rand-real-proc)))))))))
-
-;; Draw from poisson distribution with mean L, variance L.
-;; For small L, we use Knuth's method.  For larger L, we use rejection
-;; method by Atkinson, The Computer Generation of Poisson Random Variables,
-;; J. of the Royal Statistical Society Series C (Applied Statistics), 28(1),
-;; pp29-35, 1979.  The code here is a port by John D Cook's C++ implementation
-;; (http://www.johndcook.com/stand_alone_code.html )
-
-;; NOTE: this implementation calculates and stores a table of log(n!) on first invocation of L >= 36
-;; and therefore is not entirely thread safe (should still produce correct result, but with performance hit if table
-;; is recalculated multiple times)
 (define (make-poisson-generator L)
-  (unless (and (real? L)
-               (finite? L)
-               (> L 0))
-    (error "expected L to be finite positive real number"))
-  (let ((rand-real-proc (random-source-make-reals (current-random-source))))
-   (if (< L 30)
-       (make-poisson/small rand-real-proc L)
-       (make-poisson/large rand-real-proc L))))
-
-;; private
-(define (make-poisson/small rand-real-proc L)
-  (lambda ()
-    (do ((exp-L (exp (- L)))
-         (k 0 (+ k 1))
-         (p 1.0 (* p (rand-real-proc))))
-        ((<= p exp-L) (- k 1)))))
-
-;; private
-(define (make-poisson/large rand-real-proc L)
-  (let* ((c (- 0.767 (/ 3.36 L)))
-         (beta (/ PI (sqrt (* 3 L))))
-         (alpha (* beta L))
-         (k (- (log c) L (log beta))))
-    (define (loop)
-      (let* ((u (rand-real-proc))
-             (x (/ (- alpha (log (/ (- 1.0 u) u))) beta))
-             (n (exact (floor (+ x 0.5)))))
-        (if (< n 0)
-            (loop)
-            (let* ((v (rand-real-proc))
-                   (y (- alpha (* beta x)))
-                   (t (+ 1.0 (exp y)))
-                   (lhs (+ y (log (/ v (* t t)))))
-                   (rhs (+ k (* n (log L)) (- (log-of-fact n)))))
-              (if (<= lhs rhs)
-                  n
-                  (loop))))))
-    loop))
-
-;; private
-;; log(n!) table for n 1 to 256. Vector, where nth index corresponds to log((n+1)!)
-;; Computed on first invocation of `log-of-fact`
-(define log-fact-table #f)
-
-;; private
-;; computes log-fact-table
-;; log(n!) = log((n-1)!) + log(n)
-(define (make-log-fact-table!)
-  (define table (make-vector 256))
-  (vector-set! table 0 0)
-  (do ((i 1 (+ i 1)))
-      ((> i 255) #t)
-      (vector-set! table i (+ (vector-ref table (- i 1))
-                              (log (+ i 1)))))
-  (set! log-fact-table table))
-
-;; private
-;; returns log(n!)
-;; adapted from https://www.johndcook.com/blog/2010/08/16/how-to-compute-log-factorial/
-(define (log-of-fact n)
-  (when (not log-fact-table)
-    (make-log-fact-table!))
-  (cond
-    ((<= n 1) 0)
-    ((<= n 256) (vector-ref log-fact-table (- n 1)))
-    (else (let ((x (+ n 1)))
-           (+ (* (- x 0.5)
-                 (log x))
-              (- x)
-              (* 0.5
-                 (log (* 2 PI)))
-              (/ 1.0 (* x 12.0)))))))
-
-
+  (make-random-poissons mu: L random: (current-random-real)))
 
 (define (gsampling . generators-lst)
   (let ((gen-vec (list->vector generators-lst))
@@ -392,160 +249,3 @@
           (pick)))))
 
 
-;;; Code for binomial random variable generation.
-;;; Written by Brad Lucier, lucier@math.purdue.edu
-
-;;; binomial-geometric is somewhat classical, the
-;;; "First waiting time algorithm" from page 525 of
-;;; Devroye, L. (1986), Non-Uniform Random Variate
-;;; Generation, Springer-Verlag, New York.
-
-;;; binomial-rejection is algorithm BTRS from
-;;; Hormann, W. (1993), The generation of binomial
-;;; random variates, Journal of Statistical Computation
-;;; and Simulation, 46:1-2, 101-110,
-;;; DOI: https://doi.org/10.1080/00949659308811496
-;;; stirling-tail and BTRD (mentioned in the comments)
-;;; are also from that paper.
-
-;;; Another implementation of the same algorithm is at
-;;; https://github.com/tensorflow/tensorflow/blob/master/tensorflow/core/kernels/random_binomial_op.cc
-;;; That implementation pointed out at least two bugs in the
-;;; BTRS paper.
-
-(define (stirling-tail k)
-  ;; Computes
-  ;;
-  ;; \log(k!)-[\log(\sqrt{2\pi})+(k+\frac12)\log(k+1)-(k+1)]
-  ;;
-  (let ((small-k-table
-         ;; Computed using computable reals package
-         ;; Matches values in paper, which are given
-         ;; for 0\leq k < 10
-         '#(.08106146679532726
-            .0413406959554093
-            .02767792568499834
-            .020790672103765093
-            .016644691189821193
-            .013876128823070748
-            .01189670994589177
-            .010411265261972096
-            .009255462182712733
-            .00833056343336287
-            .007573675487951841
-            .00694284010720953
-            .006408994188004207
-            .0059513701127588475
-            .005554733551962801
-            .0052076559196096404
-            .004901395948434738
-            .004629153749334028
-            .004385560249232324
-            .004166319691996922)))
-    (if (< k 20)
-        (vector-ref small-k-table k)
-        ;; the correction term (+ (/ (* 12 (+ k 1))) ...)
-        ;; in Stirling's approximation to log(k!)
-        (let* ((inexact-k+1 (inexact (+ k 1)))
-               (inexact-k+1^2 (square inexact-k+1)))
-          (/ (- #i1/12
-                (/ (- #i1/360
-                      (/ #i1/1260 inexact-k+1^2))
-                   inexact-k+1^2))
-             inexact-k+1)))))
-
-(define (make-binomial-generator n p)
-  (if (not (and (real? p)
-                (<= 0 p 1)
-                (exact-integer? n)
-                (positive? n)))
-       (error "make-binomial-generator: Bad parameters: " n p)
-       (cond ((< 1/2 p)
-              (let ((complement (make-binomial-generator n (- 1 p))))
-                (lambda ()
-                  (- n (complement)))))
-             ((zero? p)
-              (lambda () 0))
-             ((< (* n p) 10)
-              (binomial-geometric n p))
-             (else
-              (binomial-rejection n p)))))
-
-(define (binomial-geometric n p)
-  (let ((geom (make-geometric-generator p)))
-    (lambda ()
-      (let loop ((X -1)
-                 (sum 0))
-        (if (< n sum)
-            X
-            (loop (+ X 1)
-                  (+ sum (geom))))))))
-
-(define (binomial-rejection n p)
-  ;; call when p <= 1/2 and np >= 10
-  ;; Use notation from the paper
-  (let* ((spq
-          (inexact (sqrt (* n p (- 1 p)))))
-         (b
-          (+ 1.15 (* 2.53 spq)))
-         (a
-          (+ -0.0873
-             (* 0.0248 b)
-             (* 0.01 p)))
-         (c
-          (+ (* n p) 0.5))
-         (v_r
-          (- 0.92
-             (/ 4.2 b)))
-         (alpha
-          ;; The formula in BTRS has 1.5 instead of 5.1;
-          ;; The formula for alpha in algorithm BTRD and Table 1
-          ;; and the tensorflow code uses 5.1, so we use 5.1
-          (* (+ 2.83 (/ 5.1 b)) spq))
-         (lpq
-          (log (/ p (- 1 p))))
-         (m
-          (exact (floor (* (+ n 1) p))))
-         (rand-real-proc
-          (random-source-make-reals (current-random-source))))
-    (lambda ()
-      (let loop ()
-        (let* ((u (rand-real-proc))
-               (v (rand-real-proc))
-               (u
-                (- u 0.5))
-               (us
-                (- 0.5 (abs u)))
-               (k
-                (exact
-                 (floor
-                  (+ (* (+ (* 2. (/ a us)) b) u) c)))))
-          (cond ((or (< k 0)
-                     (< n k))
-                 (loop))
-                ((and (<= 0.07 us)
-                      (<= v v_r))
-                 k)
-                (else
-                 (let ((v
-                        ;; The tensorflow code notes that BTRS doesn't have
-                        ;; this logarithm; BTRS is incorrect (see BTRD, step 3.2)
-                        (log (* v (/ alpha
-                                     (+ (/ a (square us)) b))))))
-                   (if (<=  v
-                            (+ (* (+ m 0.5)
-                                  (log (* (/ (+ m 1.)
-                                             (- n m -1.)))))
-                               (* (+ n 1.)
-                                  (log (/ (- n m -1.)
-                                          (- n k -1.))))
-                               (* (+ k 0.5)
-                                  (log (* (/ (- n k -1.)
-                                             (+ k 1.)))))
-                               (* (- k m) lpq)
-                               (- (+ (stirling-tail m)
-                                     (stirling-tail (- n m)))
-                                  (+ (stirling-tail k)
-                                     (stirling-tail (- n k))))))
-                       k
-                       (loop))))))))))