;;;; srfi.27.normals.scm
;;;; Kon Lovett, Sep '23
;;;; Kon Lovett, Mar '21
;;;; Kon Lovett, Dec '17
;;;; Kon Lovett, May '06

(module srfi.27.normals

(;export
  *make-random-normals
  make-random-normals)

(import scheme
  (chicken base)
  (chicken type)
  (only (type-checks-numbers scheme) check-real)
  (only type-checks-structured check-procedure)
  (srfi 27)
  srfi-27-distributions-support)

;;;

(include-relative "srfi-27-common-types")

(: *make-random-normals (number number random-real-function -> number-function))
(: make-random-normals (#!key (mu number) (sigma number) (randoms procedure) -> number-function procedure))

;;; Normal distribution

;; Knuth's "The Art of Computer Programming", Vol. II, 2nd ed.,
;; Algorithm P of Section 3.4.1.C.

(define (*make-random-normals mu sigma randoms)
  (let ((next (the maybe-number #f)))
    (lambda ()
      (if next
          (let ((result next))
            (set! next #f)
            (+ mu (* sigma result)) )
          (let loop ()
            (let* ((v1 (- (* 2 (randoms)) 1) )
                   (v2 (- (* 2 (randoms)) 1) )
                   (s (+ (* v1 v1) (* v2 v2)) ) )
              (if (<= 1 s)
                  (loop)
                  (let ((scale (sqrt (/ (* -2 (log s)) s))))
                    (set! next (* scale v2))
                    (+ mu (* sigma scale v1)) ) ) ) ) ) ) ) )

(define (make-random-normals #!key mu sigma randoms)
  (if (not mu) (set! mu 0)
      (check-real 'make-random-normals mu 'mu))
  (if (not sigma) (set! sigma 1)
      (check-nonzero-real 'make-random-normals sigma 'sigma))
  (if (not randoms) (set! randoms (current-random-real))
      (check-procedure 'make-random-normals randoms 'randoms))
  (values (*make-random-normals mu sigma randoms)
          (lambda () (values mu sigma randoms))) )

) ;module srfi.27.normals