ccm@warhol.ads.com (Chris McConnell) (07/03/90)
This random number generator is seedable, has a periord of 10^30
passes all of the statistical tests and is four times faster than the
built-in ones provided in Allegro and Lucid.
;;; -*- Mode: LISP; Package: random; Syntax: Common-lisp; Base: 10.; -*-
;;;%Header
;;;----------------------------------------------------------------------------
;;;
;;; Pseudo-random number generator
;;;
;;; Author: Chris McConnell, ccm@cs.cmu.edu
;;;
;;; This file implements a portable pseudo-random number generator for
;;; Common LISP. It has been converted from a C program that was
;;; converted from a FORTRAN program. I did not pick the variable
;;; names or pretend to have figured out how it works. The
;;; correctness of the generator can be verified by the TEST function
;;; at the end of the file.
;;;
;;; Original C header:
;;;
;;; This is the random number generator proposed by George Marsaglia in
;;; Florida State University Report: FSU-SCRI-87-50
;;;
;;; This Random Number Generator is based on the algorithm in a FORTRAN
;;; version published by George Marsaglia and Arif Zaman, Florida State
;;; University; ref.: see original comments below.
;;;
;;; At the fhw (Fachhochschule Wiesbaden, W.Germany), Dept. of Computer
;;; Science, we have written sources in further languages (C, Modula-2
;;; Turbo-Pascal(3.0, 5.0), Basic and Ada) to get exactly the same test
;;; results compared with the original FORTRAN version.
;;; April 1989
;;;
;;; Karl-L. Noell <NOELL@DWIFH1.BITNET>
;;; and Helmut Weber <WEBER@DWIFH1.BITNET>
;;;
(in-package 'random :nicknames '(r))
(export '(*state*
seed-state random-one random-float random-integer random-range))
;;;%Globals
(defvar *STATE* nil "The default state structure.")
;;;%%State
(defstruct (STATE (:print-function print-state))
"This contains random state for a state. The names of slots are the
same as the variable names in the orginal program."
(ij 1802 :type fixnum)
(kl 9373 :type fixnum)
(u (make-array 97) :type (simple-vector 97))
(c (/ 362436.0 16777216.0) :type single-float)
(cd (/ 7654321.0 16777216.0) :type single-float)
(cm (/ 16777213.0 16777216.0) :type single-float)
(i97 96 :type fixnum)
(j97 32 :type fixnum))
;;;
(defun PRINT-STATE (state stream level)
(declare (ignore level))
(format stream "#<State ~A ~A>" (state-ij state) (state-kl state)))
;;; %Interface
;;; This is the initialization routine for the random number generator
;;; NOTE: The seed variables can have values between: 0 <= IJ <= 31328
;;; 0 <= KL <= 30081
;;; The random number sequences created by these two seeds are of sufficient
;;; length to complete an entire calculation with. For example, if sveral
;;; different groups are working on different parts of the same calculation,
;;; each group could be assigned its own IJ seed. This would leave each group
;;; with 30000 choices for the second seed. That is to say, this random
;;; number generator can create 900 million different subsequences -- with
;;; each subsequence having a length of approximately 10^30.
;;;
(defun SEED-STATE (&optional (ij (mod (get-internal-real-time) 31329))
(kl (mod (get-internal-real-time) 30082)))
"Given the seed values 0 <= IJ <= 31328 and 0 <= KL <= 30081,
generate a state structure, set it as the default state and return it."
(declare (optimize (speed 3) (safety 0)))
(let* ((state (make-state :ij ij :kl kl))
(vector (state-u state))
(i (+ (mod (truncate ij 177) 177) 2))
(j (+ (mod ij 177) 2))
(k (1+ (mod (truncate kl 169) 178)))
(l (mod kl 169)))
(declare (type fixnum i j k l)
(type (simple-vector 97) vector))
(dotimes (ii 97)
(let ((s 0.0)
(t1 0.5))
(declare (type single-float s t1))
(dotimes (jj 24)
(let ((m (mod (* (mod (* i j) 179) k) 179)))
(declare (type fixnum m))
(setq i j
j k
k m
l (mod (1+ (* 53 l)) 169))
(when (>= (mod (* l m) 64) 32)
(incf s t1))
(setq t1 (* 0.5 t1))))
(setf (svref vector ii) s)))
(setq *state* state)))
;;;
(defun RANDOM-ONE (&optional (state *state*))
"Return a random value between 0.0 and 1.0 from STATE."
(declare (optimize (speed 3) (safety 0)))
(let* ((u (state-u state))
(uni (- (svref u (state-i97 state)) (svref u (state-j97 state)))))
(declare (type simple-vector u)
(type single-float uni))
(when (minusp uni)
(incf uni))
(setf (svref u (state-i97 state)) uni)
(when (minusp (decf (state-i97 state)))
(setf (state-i97 state) 96))
(when (minusp (decf (state-j97 state)))
(setf (state-j97 state) 96))
(when (minusp (decf (state-c state) (state-cd state)))
(incf (state-c state) (state-cm state)))
(when (minusp (decf uni (state-c state)))
(incf uni))
uni))
;;;
(defun RANDOM-FLOAT (n &optional (state *state*))
"Return a random float between 0 and N.
Use seed-state to initialize a pseudo-random sequence."
(declare (optimize (speed 3) (safety 0)))
(* n (random-one state)))
;;;
(defun RANDOM-INTEGER (n &optional (state *state*))
"Return a random integer between 0 and N.
Use seed-state to initialize a pseudo-random sequence."
(declare (optimize (speed 3) (safety 0)))
(truncate (* n (random-one state))))
;;;
(defun RANDOM-NUMBER (n &optional (state *state*))
"Return a random number between 0 and N with the same TYPE as N.
Use seed-state to initialize a pseudo-random sequence."
(declare (optimize (speed 3) (safety 0)))
(if (floatp n)
(random-float n state)
(random-integer n state)))
;;;
(defun RANDOM-RANGE (min max &optional (state *state*))
"Return a number between MIN and MAX with the same type.
Use seed-state to initialize a pseudo-random sequence."
(declare (optimize (speed 3) (safety 0)))
(let ((number (+ min (* (- max min) (random-one state)))))
(if (floatp min)
number
(truncate number))))
;;;%Test code
(defun TEST ()
"Test the random number generator. It should print out n = n T six
times if it works correctly."
(let ((state (seed-state 1802 9373)))
(time (dotimes (x 20000)
(random-one state)))
(dotimes (x 6)
(let ((value (round (* 4096 4096 (random-one state))))
(should-be (svref '#(6533892 14220222 7275067
6172232 8354498 10633180) x)))
(format t "~&~9A = ~9A ~A"
value should-be (= value should-be))))))
;;;%Initial seed
(setq *state* (seed-state))