# Random numbers

• To: franz-friends@Berkeley
• Subject: Random numbers
• From: Ken Anderson <kanderso@bbn-vax>
• Date: Thu, 14 Jul 83 19:47:19 GMT
• Cc: Malcolm.McRoberts@CMU-RI-ISL
• Original-date: 14 Jul 1983 15:47:19 EDT (Thursday)

```Here is a random number generator i use.  It seems to work fairly well, but i have
not looked to closely at the statistics.  Since it will occasionally require
bignums, it is probably not the fastest either.  I was just looking for something
that was likely to be portable between LISP's.
I would be very interested in hearing your evaluation of it.

k

;;; RANDOM NUMBERS
(declare (macros t))
(include math.h)

(defvar \$uniform-a 16807) ; = 7^5
(defvar \$mersenne-prime 2147483647) ; = 2^31 - 1
(defvar \$mersenne-prime-1 (- \$mersenne-prime 1))

(defmacro with-seed (s-newseed . body)
; evaluates body with the seed of the random numbers set to s-newseed.
; the value of s-newseed is updated.  Thus this is a way of
; Keepining several sequences of random numbers with their own seeds
`(let ((\$uniform-seed ,s-newseed))
(prog1 (progn ,@body)
(setq ,s-newseed \$uniform-seed))))

(defun uniform-basic (previous-fixnum)
; -> a fixnum 0 < fixnum < 2^31 - 1
; Repeated calls will generate fixnums in the range
; 1 -> 2^31 - 2.

; The basic idea is new = A^k * old (mod p)
; where A is a primitive root of p, k is not a factor of  p-1
; and p is a large prime.

; This is a good random number generator but is not be the fastest!
; On FRANZ LISP, and LISP MACHINE it will require bignums since
; (* \$uniform-a previous-fixnum) can have 46 bits in it. Also the remainder
; can be done more efficiently.
; See: Linus Schrage, A More Portable Fortran Random Number Generator,
;      ACM Trans. Math. Soft., V5, No. 2, p 132-138, 1979.
(remainder (*\$ \$uniform-a previous-fixnum) \$mersenne-prime))

(defvar \$uniform-seed 53) ; 0 < fixnum < \$mersenne-prime-1

(defun uniform ()
; -> the next uniform random number in the sequence
; To have your own sequence, rebind \$uniform-seed.
(setq \$uniform-seed (uniform-basic \$uniform-seed)))

(defun uniform-between (low-num high-num)
; -> a uniform random  number, x, low-num <= x <= high-num
; If low-num and high-num are fixnums, a fixnum is returned.
(cond ((not (and (fixp low-num) (fixp high-num)))
(+\$ low-num (*\$ (//\$ (uniform) (float \$mersenne-prime-1))
(-\$ high-num low-num))))
(t (+ low-num (// (uniform)
(// \$mersenne-prime-1 (max 1 (- high-num low-num -1))))))))

(defun gaussian-random-1 ()
; -> a gaussian random variable with mean 0.0 and
; standard deviation 1.0.
; Good tails.
(*\$ (sqrt (*\$ -2.0 (log (uniform-between 0.0 1.0))))
(sin (*\$ \$2pi (uniform-between 0.0 1.0)))))

(defun gaussian-random (mean standard-deviation)
(+\$ mean (*\$ (gaussian-random-1) standard-deviation)))

;;(defun gaussian (x)
;;  (* (sqrt \$2pi)
;;     (exp (minus (// (square x) 2.0)))))

(defun random-yes-no (fraction-yes)
(<= (uniform-between 0.0 1.0) fraction-yes))

```