diff --git a/code/rand.lisp b/code/rand.lisp index 306279f2dfe91fa4eceb33ba07517d2922c80fce..f7d00d80b9e9916c94e8c061d2257e2ead5d56d0 100644 --- a/code/rand.lisp +++ b/code/rand.lisp @@ -1,4 +1,4 @@ -;;; -*- Mode: Lisp; Package: Lisp; Log: code.log -*- +;;; -*- Mode: Lisp; Package: Kernel -*- ;;; ;;; ********************************************************************** ;;; This code was written as part of the CMU Common Lisp project at @@ -7,70 +7,61 @@ ;;; Scott Fahlman or slisp-group@cs.cmu.edu. ;;; (ext:file-comment - "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/rand.lisp,v 1.4 1991/12/14 09:01:01 wlott Exp $") + "$Header: /Volumes/share2/src/cmucl/cvs2git/cvsroot/src/code/rand.lisp,v 1.5 1993/05/07 19:01:51 ram Exp $") ;;; ;;; ********************************************************************** ;;; -;;; Functions to random number functions for Spice Lisp -;;; Written by David Adam. +;;; Functions to random number functions for Spice Lisp ;;; -;;; The random number functions are part of the standard Spicelisp environment. +;;; Originally written by David Adam. Python tuning, better large integer +;;; randomness and efficient IEEE float support by Rob MacLachlan. ;;; -;;; ********************************************************************** -;;; -(in-package 'lisp) +(in-package "LISP") (export '(random-state random-state-p random *random-state* make-random-state)) +(in-package "KERNEL") +(export '(%random-single-float %random-double-float random-chunk + random-fixnum-max)) + + +;;;; Random state hackery: + (defconstant random-const-a 8373) (defconstant random-const-c 101010101) -(defconstant random-upper-bound (1- most-positive-fixnum)) (defconstant random-max 54) -(defconstant %fixnum-length (integer-length most-positive-fixnum)) -(defvar rand-seed 0) +;;; Inclusive upper bound on the size of fixnum kept in the state (and returned +;;; by random-chunk.) Must be even. +;;; +(defconstant random-upper-bound (- most-positive-fixnum 3)) +(defconstant random-chunk-length (integer-length random-upper-bound)) +(deftype random-chunk () `(integer 0 ,random-upper-bound)) + +(defvar rand-seed 0) (defstruct (random-state (:constructor make-random-object) (:make-load-form-fun :just-dump-it-normally)) - (j 24 :type integer) - (k 0 :type integer) + (j 24 :type index) + (k 0 :type index) (seed (make-array (1+ random-max) :initial-contents (do ((list-rands () (cons (rand1) list-rands)) (i 0 (1+ i))) - ((> i random-max) list-rands))) + ((> i random-max) list-rands) + (declare (fixnum i)))) :type simple-vector)) ;;; Generates a random number from rand-seed. (defun rand1 () - (setq rand-seed (mod (+ (* rand-seed random-const-a) random-const-c) - (1+ random-upper-bound)))) + (declare (optimize (inhibit-warnings 3))) + (setq rand-seed + (mod (+ (* rand-seed random-const-a) random-const-c) + (1+ random-upper-bound)))) (defvar *random-state* (make-random-object)) - -;;; rand3 -- Internal -;;; -;;; This function generates fixnums between 0 and random-upper-bound, -;;; inclusive For the algorithm to work random-upper-bound must be an -;;; even positive fixnum. State is the random state to use. -;;; -(defun rand3 (state) - (let ((seed (random-state-seed state)) - (j (random-state-j state)) - (k (random-state-k state))) - (declare (fixnum j k) (simple-vector seed)) - (setf (svref seed k) - (let ((a (- random-upper-bound - (svref seed - (setf (random-state-j state) - (if (= j 0) random-max (1- j)))) - (svref seed - (setf (random-state-k state) - (if (= k 0) random-max (1- k))))))) - (if (minusp a) (- a) (- random-upper-bound a)))))) - (defun copy-state (cur-state) (let ((state (make-random-object @@ -92,28 +83,121 @@ ((random-state-p state) (copy-state state)) ((eq state t) (setq rand-seed (get-universal-time)) (make-random-object)) - (t (error "Bad argument, ~A, for RANDOM-STATE." state)))) + (t (error "Argument is not a RANDOM-STATE, T or NIL: ~S" state)))) + + +;;;; Random entries: + +(declaim (start-block random %random-single-float %random-double-float + random-chunk)) + +;;; random-chunk -- Internal +;;; +;;; This function generates fixnums between 0 and random-upper-bound, +;;; inclusive. For the algorithm to work random-upper-bound must be an +;;; even positive fixnum. State is the random state to use. +;;; +(declaim (ftype (function (random-state) random-chunk) random-chunk)) +(defun random-chunk (state) + (let* ((seed (random-state-seed state)) + (j (random-state-j state)) + (k (random-state-k state)) + (a (- (- random-upper-bound + (the random-chunk + (svref seed + (setf (random-state-j state) + (if (= j 0) random-max (1- j)))))) + (the random-chunk + (svref seed + (setf (random-state-k state) + (if (= k 0) random-max (1- k)))))))) + (declare (fixnum a)) + (setf (svref seed k) + (the random-chunk (if (minusp a) (- a) (- random-upper-bound a)))))) + + +;;; %RANDOM-SINGLE-FLOAT, %RANDOM-DOUBLE-FLOAT -- Interface +;;; +;;; Handle the single or double float case of RANDOM. We generate a float +;;; between 0.0 and 1.0 by clobbering the significand of 1.0 with random bits, +;;; then subtracting 1.0. This hides the fact that we have a hidden bit. +;;; +(declaim (inline %random-single-float %random-double-float)) +(defun %random-single-float (arg state) + (declare (type (single-float (0f0)) arg) + (type (or random-state null) state)) + (* arg + (- (make-single-float + (dpb (ash (random-chunk (or state *random-state*)) + (- vm:single-float-digits random-chunk-length)) + vm:single-float-significand-byte + (single-float-bits 1.0))) + 1.0))) +;;; +(defun %random-double-float (arg state) + (declare (type (double-float (0d0)) arg) + (type (or random-state null) state)) + (let ((state (or state *random-state*))) + (* arg + (- (make-double-float + (dpb (ash (random-chunk state) + (- vm:double-float-digits random-chunk-length + vm:word-bits)) + vm:double-float-significand-byte + (double-float-high-bits 1d0)) + (logxor (ash (random-chunk state) + (- vm:word-bits random-chunk-length)) + (ash (random-chunk state) + (- random-chunk-length vm:word-bits)))) + 1d0)))) + + +;;;; Random integers: + +;;; Amount we overlap chunks by when building a large integer to make up for +;;; the loss of randomness in the low bits. +;;; +(defconstant random-integer-overlap 3) + +;;; Extra bits of randomness that we generate before taking the value MOD the +;;; limit, to avoid loss of randomness near the limit. +;;; +(defconstant random-integer-extra-bits 10) + +;;; Largest fixnum we can compute from one chunk of bits. +;;; +(defconstant random-fixnum-max + (1- (ash 1 (- random-chunk-length random-integer-extra-bits)))) + + +;;; %RANDOM-INTEGER -- Internal +;;; +(defun %random-integer (arg state) + (declare (type (integer 1) arg) (type random-state state)) + (let ((shift (- random-chunk-length random-integer-overlap))) + (do ((bits (random-chunk state) + (logxor (ash bits shift) (random-chunk state))) + (count (+ (integer-length arg) + (- random-integer-extra-bits shift)) + (- count shift))) + ((minusp count) + (rem bits arg)) + (declare (fixnum count))))) -(proclaim '(ftype (function (t) fixnum) rand3)) (defun random (arg &optional (state *random-state*)) "Generate a uniformly distributed pseudo-random number between zero and Arg. State, if supplied, is the random state to use." - (typecase arg - (fixnum - (unless (plusp (the fixnum arg)) - (error "Non-positive argument, ~A, to RANDOM." arg)) - (rem (the fixnum (rand3 state)) (the fixnum arg))) - (float - (unless (plusp arg) - (error "Non-positive argument, ~A, to RANDOM." arg)) - (let ((arg-length (float-digits arg))) - (* arg (/ (float (random (ash 2 arg-length) state)) - (float (ash 2 arg-length)))))) - (integer - (unless (plusp arg) - (error "Non-positive argument, ~A, to RANDOM." arg)) - (do ((tot (rand3 state) (+ (ash tot %fixnum-length) (rand3 state))) - (end (ash arg (- %fixnum-length)) - (ash end (- %fixnum-length)))) - ((zerop end) (mod tot arg)))) - (t (error "Wrong type argument, ~A, to RANDOM." arg)))) + (declare (inline %random-single-float %random-double-float)) + (cond + ((and (fixnump arg) (<= arg random-fixnum-max)) + (rem (random-chunk state) arg)) + ((typep arg 'single-float) + (%random-single-float arg state)) + ((typep arg 'double-float) + (%random-double-float arg state)) + ((integerp arg) + (%random-integer arg state)) + (t + (error 'simple-type-error :expected-type '(real (0)) :datum arg + :format-string "Argument is not a positive real number: ~S" + :format-arguments (list arg)))))