3
\$\begingroup\$

For reasons that are too dreary to detail, I need to reproduce several pseudo random number generators in different languages. It mostly involves translating ancient C and FORTRAN code. No problems with C, Java, Zig, Julia, etc. But for Lisp-derived languages, I would like to do something more idiomatic.

For example, here is a Scheme interpretation of ACM algorithm 647 for a portable PRNG to generate samples from a uniform distribution.

;; Return a random uniform deviate in the range 0 <= x < 1.0.
(define algo-647-uniform
  (let ([seed 12345])
    (lambda ()
      (let* ([k (floor (/ seed 127773))]
             [partial (- (* 16807 (- seed (* k 127773))) (* k 2836))])
        (if (< partial 0)
            (set! seed (+ partial 2147483647))
            (set! seed partial))
        (* seed 4.656612875e-10)))))

The above works and I'm satisfied with the implementation (though eager to hear suggestions for improvements).

The implementation below of a function to generate normally distributed PRNGs also works, but looks kind of ugly -- there is so much state to maintain. Is there a better way to do the following?

;; The following is an implementation of the Marsaglia polar
;; method for generating random standard normal deviates. See
;; https://en.wikipedia.org/wiki/Marsaglia_polar_method.

;; Return a random standard normal deviate.
(define gaussian-deviate-marsaglia
  (let* ([has-spare #f]
         [spare 0.0]
         [rsq 0.0] [r1 0.0] [r2 0.0]
         [nxt-rsq (lambda ()
                    (set! r1 (- (* 2.0 (algo-647-uniform)) 1.0))
                    (set! r2 (- (* 2.0 (algo-647-uniform)) 1.0))
                    (set! rsq (+ (* r1 r1) (* r2 r2)))
                    rsq)])
    (lambda ()
      (if has-spare
          (begin
            (set! has-spare #f)
            spare)
          (let loop ([_ (nxt-rsq)])
            (if (or (>= rsq 1.0) (= rsq 0.0))
                (loop (nxt-rsq))
                (let ([fac (sqrt (/ (* -2.0 (log rsq)) rsq))])
                  (set! spare (* r1 fac))
                  (set! has-spare #t)
                  (* r2 fac))))))))

In Chez 9.5.8, the first 10 results should be:

-1.0580380669115383
-0.5790254729247644
0.8434589541668004
-0.6708443571574382
-0.22644041228981196
-0.07818079860601053
-0.7443285279492631
-0.5232388154010481
-0.3300334725931046
0.41341813121639936
\$\endgroup\$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.