Random numbers

[Ken Anderson <kanderso@bbn-vax>]

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