[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
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))