random-number-array.lisp

;; based on https://marui.hatenablog.com/entry/2022/12/13/222030

(defun mapfn-into (arr fn &rest args)
  (loop for i below (array-total-size arr)
        do (setf (row-major-aref arr i) (apply fn args))
        finally (return arr)))

(defun box-muller-random-sampling ()
  "Return a pair of independent, standard, normally distributed
 (zero expectation, unit variance) random number, using a random
 number sampling method known as the Box–Muller transform."
  (flet ((rndm () (loop for x = (random 1.0) when (plusp x) return x))
         (bmrs (u1 u2)
           (values (* (sqrt (* -2 (log u1))) (cos (* 2 pi u2)))
                   (* (sqrt (* -2 (log u1))) (sin (* 2 pi u2))))))
    (bmrs (rndm) (rndm))))

(let ((stock nil))
  (defun randn ()
    "Return a random number following uniform and normal distribion."
    (if (null stock)
      (multiple-value-bind (r1 r2) (box-muller-random-sampling)
        (progn (setq stock r2) r1))
      (prog1 stock (setq stock nil)))))

(defun make-randn-array (dims)
  (mapfn-into (make-array dims) #'randn))