[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

*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))

- Prev by Date:
**random numbers** - Next by Date:
**Franz & SUNs** - Previous by thread:
**random numbers** - Next by thread:
**Franz & SUNs** - Index(es):