Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Two declarations to quiet compiler when building on Allegro Common Lisp #3

Merged
merged 5 commits into from
Nov 4, 2013
Merged

Two declarations to quiet compiler when building on Allegro Common Lisp #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
61f5b17
Add a SPECIAL declaration to quiet ACL compiler.
rpgoldman Jun 20, 2013
0904584
Added IGNORE declaration.
rpgoldman Jun 20, 2013
125c812
Added a terse PRINT-METHOD for MT-RANDOM-STATE objects when *PRINT-RE…
rpgoldman Jul 24, 2013
66f9874
Document the argument to MAKE-MT-RANDOM-STATE.
rpgoldman Jul 25, 2013
0752d69
Added MAKE-LOAD-FORM method for MT-RANDOM-STATE.
rpgoldman Oct 11, 2013
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
183 changes: 101 additions & 82 deletions jmt.lisp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,21 +70,30 @@
(defconstant +mt-lower-mask+ #x7FFFFFFF "least significant r bits")

(defstruct (mt-random-state
(:constructor mt-internal-make-random-state))
(:constructor mt-internal-make-random-state))
;; Could have avoided MTI, which is an index into ARR, with a
;; fill pointer in ARR. MTI more closely follows the reference
;; implementation.
;; ARR corresponts to "mt[]" in the reference implementation.
;; Probably should have called it MT after all. Oh well.
mti ; index into ARR
arr) ; array of numbers
mti ; index into ARR
arr) ; array of numbers

(defmethod print-object ((obj mt-random-state) stream)
(if *print-readably*
(call-next-method)
(print-unreadable-object (obj stream :type t :identity t))))

(defmethod make-load-form ((obj mt-random-state) &optional environment)
(declare (ignore environment))
(make-load-form-saving-slots obj :slot-names '(mti arr)))

(labels
((next-seed (n) (mod (1+ (* 69069 n)) +mt-k2^32+))
(get-hi16 (n) (logand n #xFFFF0000))
(next-elt (n)
(logior (get-hi16 n)
(ash (get-hi16 (next-seed n)) -16))))
(logior (get-hi16 n)
(ash (get-hi16 (next-seed n)) -16))))
(defun mt-make-random-state-integer (n)
"Use the single integer to expand into a bunch of
integers to use as an MT-RANDOM-STATE.
Expand All @@ -94,12 +103,12 @@ MAKE-MT-RANDOM-STATE unless specific circumstances dictate otherwise."
(mt-internal-make-random-state
:mti +mt-n+
:arr (make-array
+mt-n+
:element-type 'integer
:initial-contents (do ((i 0 (1+ i))
(sd n (next-seed (next-seed sd)))
(lst () (cons (next-elt sd) lst)))
((>= i +mt-n+) (nreverse lst)))))))
+mt-n+
:element-type 'integer
:initial-contents (do ((i 0 (1+ i))
(sd n (next-seed (next-seed sd)))
(lst () (cons (next-elt sd) lst)))
((>= i +mt-n+) (nreverse lst)))))))

(let ((some-number 0))
(defun mt-make-random-state-random ()
Expand All @@ -108,42 +117,52 @@ random, value.
This is mostly an internal function. I recommend using
MAKE-MT-RANDOM-STATE unless specific circumstances dictate otherwise."
(mt-make-random-state-integer (+ (get-universal-time)
(incf some-number)))))
(incf some-number)))))

(defun make-mt-random-state (&optional state)
"Analogous to Common Lisp's MAKE-RANDOM-STATE except that this function
works on random states for JMT's Mersenne Twister implementation.

Optional state argument is interpreted as follows:
if T -- make a new random state.
if NIL -- return a copy of the current value of *MT-RANDOM-STATE*
if an integer -- return a random state created by expanding the
integer to a full random state
if a list or array -- create a new random state object using a copy of
that sequence. CF MT-INTERNAL-MAKE-RANDOM-STATE
if a random state -- return a copy of the input state.

AJR-comment: this version is better for statistical computations,
since it can be used by an externally defined seed. However, it
must be optimized, being the core of any pseudo-stochastic
algorithm."
(declare (special *mt-random-state*))
(cond ((eq state t) (mt-make-random-state-random))
((null state)
;; For NIL, return a copy of the current state.
(make-mt-random-state *mt-random-state*))
((integerp state)
;; Expand the integer STATE into controlled junk that is an
;; MT RANDOM STATE.
(mt-make-random-state-integer state))
((typep state 'sequence)
;; It's a list or an array. It must be of length +MT-N+, & it
;; must contain integers. We'll create a random state object
;; using a copy of that sequence.
(assert state) ; should have caught NIL earlier
(assert (eql (length state) +mt-n+))
(assert (not (find-if #'integerp state)))
(mt-internal-make-random-state
:mti 0
:arr (copy-seq (coerce state 'array))))
((mt-random-state-p state)
;; Return a copy of state. It is an instance of MT-RANDOM-STATE.
(mt-internal-make-random-state
:mti (mt-random-state-mti state)
:arr (copy-seq (mt-random-state-arr state))))
(t
;; For anything else, error.
(cerror "STATE should not have a value of ~A" state))))
((null state)
;; For NIL, return a copy of the current state.
(make-mt-random-state *mt-random-state*))
((integerp state)
;; Expand the integer STATE into controlled junk that is an
;; MT RANDOM STATE.
(mt-make-random-state-integer state))
((typep state 'sequence)
;; It's a list or an array. It must be of length +MT-N+, & it
;; must contain integers. We'll create a random state object
;; using a copy of that sequence.
(assert state) ; should have caught NIL earlier
(assert (eql (length state) +mt-n+))
(assert (not (find-if #'integerp state)))
(mt-internal-make-random-state
:mti 0
:arr (copy-seq (coerce state 'array))))
((mt-random-state-p state)
;; Return a copy of state. It is an instance of MT-RANDOM-STATE.
(mt-internal-make-random-state
:mti (mt-random-state-mti state)
:arr (copy-seq (mt-random-state-arr state))))
(t
;; For anything else, error.
(cerror "STATE should not have a value of ~A" state))))

(defvar *mt-random-state* (make-mt-random-state t)
"Unlike the reference implementation, we'll initialize the random
Expand All @@ -161,47 +180,47 @@ MT-GENRAND function for clarity."
(let (y kk)
(setq kk 0)
(do ()
((>= kk (- +mt-n+ +mt-m+)))
(setq y (logior
(logand (aref (mt-random-state-arr *mt-random-state*)
kk)
+mt-upper-mask+)
(logand (aref (mt-random-state-arr *mt-random-state*)
(1+ kk))
+mt-lower-mask+)))
(setf (aref (mt-random-state-arr *mt-random-state*) kk)
(logxor
(aref (mt-random-state-arr *mt-random-state*) (+ kk +mt-m+))
(ash y -1)
(aref mag01 (logand y 1))))
(incf kk))
((>= kk (- +mt-n+ +mt-m+)))
(setq y (logior
(logand (aref (mt-random-state-arr *mt-random-state*)
kk)
+mt-upper-mask+)
(logand (aref (mt-random-state-arr *mt-random-state*)
(1+ kk))
+mt-lower-mask+)))
(setf (aref (mt-random-state-arr *mt-random-state*) kk)
(logxor
(aref (mt-random-state-arr *mt-random-state*) (+ kk +mt-m+))
(ash y -1)
(aref mag01 (logand y 1))))
(incf kk))
(do ()
((>= kk (- +mt-n+ 1)))
(setq y (logior
(logand
(aref (mt-random-state-arr *mt-random-state*) kk)
+mt-upper-mask+)
(logand
(aref (mt-random-state-arr *mt-random-state*) (1+ kk))
+mt-lower-mask+)))
(setf (aref (mt-random-state-arr *mt-random-state*) kk)
(logxor (aref (mt-random-state-arr *mt-random-state*)
(+ kk (- +mt-m+ +mt-n+)))
(ash y -1)
(aref mag01 (logand y 1))))
(incf kk))
((>= kk (- +mt-n+ 1)))
(setq y (logior
(logand
(aref (mt-random-state-arr *mt-random-state*) kk)
+mt-upper-mask+)
(logand
(aref (mt-random-state-arr *mt-random-state*) (1+ kk))
+mt-lower-mask+)))
(setf (aref (mt-random-state-arr *mt-random-state*) kk)
(logxor (aref (mt-random-state-arr *mt-random-state*)
(+ kk (- +mt-m+ +mt-n+)))
(ash y -1)
(aref mag01 (logand y 1))))
(incf kk))
(setq y (logior
(logand
(aref (mt-random-state-arr *mt-random-state*) (- +mt-n+ 1))
+mt-upper-mask+)
(logand
(aref (mt-random-state-arr *mt-random-state*) 0)
+mt-lower-mask+)))
(logand
(aref (mt-random-state-arr *mt-random-state*) (- +mt-n+ 1))
+mt-upper-mask+)
(logand
(aref (mt-random-state-arr *mt-random-state*) 0)
+mt-lower-mask+)))
(setf (aref (mt-random-state-arr *mt-random-state*) (- +mt-n+ 1))
(logxor
(aref (mt-random-state-arr *mt-random-state*) (- +mt-m+ 1))
(ash y -1)
(aref mag01 (logand y 1))))
(logxor
(aref (mt-random-state-arr *mt-random-state*) (- +mt-m+ 1))
(ash y -1)
(aref mag01 (logand y 1))))
(setf (mt-random-state-mti *mt-random-state*) 0))
'mt-refill))

Expand All @@ -223,7 +242,7 @@ MT-GENRAND function for clarity."
(when (>= (mt-random-state-mti *mt-random-state*) +mt-n+)
(mt-refill))
(let ((y (aref (mt-random-state-arr *mt-random-state*)
(mt-random-state-mti *mt-random-state*))))
(mt-random-state-mti *mt-random-state*))))
(incf (mt-random-state-mti *mt-random-state*))
;; The following separate, explicit SETQ & other expressions
;; could be compacted/optimized into a single arithmetic expression
Expand All @@ -232,9 +251,9 @@ MT-GENRAND function for clarity."
;; the statements in the C program, mt19937int.c.
(setq y (logxor y (mt-tempering-shift-u y)))
(setq y (logxor y (logand (mt-tempering-shift-s y)
mt-tempering-mask-b)))
mt-tempering-mask-b)))
(setq y (logxor y (logand (mt-tempering-shift-t y)
mt-tempering-mask-c)))
mt-tempering-mask-c)))
(setq y (logxor y (mt-tempering-shift-l y)))
y)))

Expand All @@ -247,10 +266,10 @@ MT-GENRAND function for clarity."
(setq *mt-random-state* (make-mt-random-state state)))
(if (integerp n)
(mod (do ((bits-needed (log n 2) )
(bit-count 0 (+ 32 bit-count))
(r 0 (+ (ash r 32) (mt-genrand))))
((>= bit-count bits-needed) r))
n)
(bit-count 0 (+ 32 bit-count))
(r 0 (+ (ash r 32) (mt-genrand))))
((>= bit-count bits-needed) r))
n)
(* (mt-genrand) +mt-k-inverse-2^32f+ n)))

;;; --- end of file ---
11 changes: 6 additions & 5 deletions tests.lisp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,24 +21,25 @@
(defun test-dist (f &optional (n 1000))
(let ((lst (genlist f n)))
(values (float (mean lst))
(float (var lst)))))
(float (var lst)))))

(defvar *n* (* 100 1000 1000))
(defun test-gamma-speed (&optional (n *n*))
(time (dotimes (i n)
(random-gamma1 2.4d0 0.5d0))))
(random-gamma1 2.4d0 0.5d0))))

(defun test-gamma-mt-speed (&optional (n *n*))
(time (dotimes (i n)
(random-gamma-mt 2.4d0 0.5d0))))
(random-gamma-mt 2.4d0 0.5d0))))

(defun test-zigg-speed (&optional (n *n*))
(time (dotimes (i n)
(random-normal-ziggurat 0d0 1d0))))
(random-normal-ziggurat 0d0 1d0))))

(defun test-binomial (&optional (n *n*))
(test-dist #'(lambda (x)
(random-binomial 0.3d0 10)) n))
(declare (ignore x))
(random-binomial 0.3d0 10)) n))

#+sbcl
(defun profile-gamma (&optional (n (* 100 1000)))
Expand Down