home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / bsd_srcs / usr.bin / lisp / lispnews / text0137.txt < prev    next >
Encoding:
Text File  |  1985-11-10  |  2.7 KB  |  76 lines
  1. Here is a random number generator i use.  It seems to work fairly well, but i have
  2. not looked to closely at the statistics.  Since it will occasionally require
  3. bignums, it is probably not the fastest either.  I was just looking for something
  4. that was likely to be portable between LISP's.
  5. I would be very interested in hearing your evaluation of it.
  6.  
  7. k
  8.  
  9. ;;; RANDOM NUMBERS
  10. (declare (macros t))
  11. (include math.h)
  12.  
  13. (defvar $uniform-a 16807) ; = 7^5
  14. (defvar $mersenne-prime 2147483647) ; = 2^31 - 1
  15. (defvar $mersenne-prime-1 (- $mersenne-prime 1))
  16.  
  17. (defmacro with-seed (s-newseed . body)
  18.   ; evaluates body with the seed of the random numbers set to s-newseed.
  19.   ; the value of s-newseed is updated.  Thus this is a way of
  20.   ; Keepining several sequences of random numbers with their own seeds
  21.   `(let (($uniform-seed ,s-newseed))
  22.     (prog1 (progn ,@body) 
  23.            (setq ,s-newseed $uniform-seed))))
  24.  
  25. (defun uniform-basic (previous-fixnum)
  26.   ; -> a fixnum 0 < fixnum < 2^31 - 1
  27.   ; Repeated calls will generate fixnums in the range
  28.   ; 1 -> 2^31 - 2.
  29.  
  30.   ; The basic idea is new = A^k * old (mod p)
  31.   ; where A is a primitive root of p, k is not a factor of  p-1
  32.   ; and p is a large prime.
  33.  
  34.   ; This is a good random number generator but is not be the fastest!
  35.   ; On FRANZ LISP, and LISP MACHINE it will require bignums since
  36.   ; (* $uniform-a previous-fixnum) can have 46 bits in it. Also the remainder
  37.   ; can be done more efficiently.
  38.   ; See: Linus Schrage, A More Portable Fortran Random Number Generator,
  39.   ;      ACM Trans. Math. Soft., V5, No. 2, p 132-138, 1979.
  40.   (remainder (*$ $uniform-a previous-fixnum) $mersenne-prime))
  41.  
  42. (defvar $uniform-seed 53) ; 0 < fixnum < $mersenne-prime-1
  43.  
  44. (defun uniform ()
  45.   ; -> the next uniform random number in the sequence
  46.   ; To have your own sequence, rebind $uniform-seed.
  47.   (setq $uniform-seed (uniform-basic $uniform-seed)))
  48.  
  49. (defun uniform-between (low-num high-num)
  50.   ; -> a uniform random  number, x, low-num <= x <= high-num
  51.   ; If low-num and high-num are fixnums, a fixnum is returned.
  52.   (cond ((not (and (fixp low-num) (fixp high-num)))
  53.      (+$ low-num (*$ (//$ (uniform) (float $mersenne-prime-1))
  54.                (-$ high-num low-num))))
  55.     (t (+ low-num (// (uniform)
  56.               (// $mersenne-prime-1 (max 1 (- high-num low-num -1))))))))
  57.  
  58. (defun gaussian-random-1 ()
  59.   ; -> a gaussian random variable with mean 0.0 and
  60.   ; standard deviation 1.0.
  61.   ; Good tails.
  62.   (*$ (sqrt (*$ -2.0 (log (uniform-between 0.0 1.0))))
  63.      (sin (*$ $2pi (uniform-between 0.0 1.0)))))
  64.  
  65. (defun gaussian-random (mean standard-deviation)
  66.   (+$ mean (*$ (gaussian-random-1) standard-deviation)))
  67.  
  68. ;;(defun gaussian (x)
  69. ;;  (* (sqrt $2pi) 
  70. ;;     (exp (minus (// (square x) 2.0)))))
  71.  
  72. (defun random-yes-no (fraction-yes)
  73.     (<= (uniform-between 0.0 1.0) fraction-yes))
  74.  
  75.  
  76.