Raymond Toy pushed to branch rtoy-xoro at cmucl / cmucl
Commits: b119b34f by Raymond Toy at 2017-12-15T09:00:38-08:00 Initial implementation of xoroshiro rng
Not yet tested or integrated.
- - - - -
1 changed file:
- + src/code/rand-xoroshiro.lisp
Changes:
===================================== src/code/rand-xoroshiro.lisp ===================================== --- /dev/null +++ b/src/code/rand-xoroshiro.lisp @@ -0,0 +1,362 @@ +;;; -*- Mode: Lisp; Package: Kernel -*- +;;; +;;; ********************************************************************** +(ext:file-comment + "$Header: src/code/rand-xoroshiro.lisp $") + +;;; +;;; ********************************************************************** +;;; +;;; Support for the xoroshiro128+ random number generator by David +;;; Blackman and Sebastiano Vigna (vigna@acm.org) + +(in-package "LISP") +(intl:textdomain "cmucl") + +(export '(random-state random-state-p random *random-state* + make-random-state)) + +(in-package "KERNEL") +(export '(%random-single-float %random-double-float random-chunk init-random-state)) + +(sys:register-lisp-feature :random-xoroshiro) + +(defun int-init-random-state (&optional (seed 5772156649015328606) state) + (let ((state (or state (make-array 2 :element-type 'double-float))) + (splitmix-state (ldb (byte 64 0) seed))) + (flet ((splitmix64 () + (let ((z (setf splitmix-state + (ldb (byte 64 0) (+ splitmix-state #x9e3779b97f4a7c15))))) + (declare (type (unsigned-byte 64) z)) + (setf z (ldb (byte 64 0) + (* (logxor z (ash z -30)) + #xbf58476d1ce4e5b9))) + (setf z (ldb (byte 64 0) + (* (logxor z (ash z -27)) + #x94d049bb133111eb))) + (logxor z (ash z -31)))) + (make-double (x) + (let ((lo (ldb (byte 32 0) x)) + (hi (ldb (byte 32 32) x))) + (kernel:make-double-float + (if (< hi #x80000000) + hi + (- hi #x100000000)) + lo)))) + (let* ((s0 (splitmix64)) + (s1 (splitmix64))) + (setf (aref state 0) (make-double s0) + (aref state 1) (make-double s1)))))) + +(defun vec-init-random-state (key &optional state) + (declare (type (array (unsigned-byte 32) (4)) key) + (type (simple-array double-float (2)) state)) + (flet ((make-double (hi lo) + (kernel:make-double-float + (if (< hi #x80000000) + hi + (- hi #x100000000)) + lo))) + (setf (aref state 0) (make-double (aref key 0) (aref key 1)) + (aref state 1) (make-double (aref key 2) (aref key 3))))) + + +(defun init-random-state (&optional (seed 5772156649015328606) state) + "Generate an random state vector from the given SEED. The seed can be + either an integer or a vector of (unsigned-byte 32)" + (declare (type (or null integer + (array (unsigned-byte 32) (*))) + seed)) + (etypecase seed + (integer + (int-init-random-state (ldb (byte 64 0) seed) state)) + ((array (unsigned-byte 32) (4)) + (vec-init-random-state seed state)))) + +(defstruct (xoro-random-state + (:constructor make-xoroshiro-object) + (:make-load-form-fun :just-dump-it-normally)) + (state (init-random-state) + :type (simple-array double-float (2))) + (rand (make-array 1 :element-type '(unsigned-byte 32) :initial-element 0) + :type (simple-array (unsigned-byte 32) (1))) + (cached-p nil :type (member t nil))) + + + +;;;; Random entries: + +;;; Size of the chunks returned by random-chunk. +;;; +(defconstant random-chunk-length 32) + +;;; random-chunk -- Internal +;;; +;;; This function generaters a 32bit integer between 0 and #xffffffff +;;; inclusive. +;;; +(declaim (inline random-chunk)) + +(defun random-chunk (rng-state) + (declare (type xoro-state rng-state) + (optimize (speed 3) (safety 0))) + (let ((cached (xoro-state-cached-p rng-state))) + (cond (cached + (setf (xoro-state-cached-p rng-state) nil) + (aref (xoro-state-rand rng-state) 0)) + (t + (let ((s (xoro-state-state rng-state))) + (declare (type (simple-array double-float (2)) s)) + (multiple-value-bind (r1 r0) + (vm::xoroshiro-next s) + (setf (aref (xoro-state-rand rng-state) 0) r1) + (setf (xoro-state-cached-p rng-state) t) + r0)))))) + +#-x86 +(defun xoroshiro-next (state) + (declare (type (simple-array double-float (2)) state)) + (flet ((rotl-55 (x1 x0) + (declare (type (unsigned-byte 32) x0 x1) + (optimize (speed 3) (safety 0))) + ;; x << 55 + (let ((sl55-h (ldb (byte 32 0) (ash x0 (- 55 32)))) + (sl55-l 0)) + ;; x >> 9 + (let ((sr9-h (ash x1 -9)) + (sr9-l (ldb (byte 32 0) + (logior (ash x0 -9) + (ash x1 23))))) + (values (logior sl55-h sr9-h) + (logior sl55-l sr9-l))))) + (rotl-36 (x1 x0) + (declare (type (unsigned-byte 32) x0 x1) + (optimize (speed 3) (safety 0))) + ;; x << 36 + (let ((sl36-h (ldb (byte 32 0) (ash x0 4)))) + ;; x >> 28 + (let ((sr28-l (ldb (byte 32 0) + (logior (ash x0 -28) + (ash x1 4)))) + (sr28-h (ash x1 -28))) + (values (logior sl36-h sr28-h) + sr28-l)))) + (shl-14 (x1 x0) + (declare (type (unsigned-byte 32) x1 x0) + (optimize (speed 3) (safety 0))) + (values (ldb (byte 32 0) + (logior (ash x1 14) + (ash x0 (- 14 32)))) + (ldb (byte 32 0) + (ash x0 14)))) + (make-double (hi lo) + (kernel:make-double-float + (if (< hi #x80000000) + hi + (- hi #x100000000)) + lo))) + (let ((s0-1 0) + (s0-0 0) + (s1-1 0) + (s1-0 0) + (r1 0) + (r0 0)) + (declare (type (unsigned-byte 32)) s0-1 s0-0 s1-1 s1-0 r1 r0) + (multiple-value-bind (x1 x0) + (kernel:double-float-bits (aref state 0)) + (setf s0-1 (ldb (byte 32 0) x1) + s0-0 x0)) + (multiple-value-bind (x1 x0) + (kernel:double-float-bits (aref state 1)) + (setf s1-1 (ldb (byte 32 0) x1) + s1-0 x0)) + + (multiple-value-prog1 + (multiple-value-bind (sum-0 c) + (bignum::%add-with-carry s0-0 s1-0 0) + (values (bignum::%add-with-carry s0-1 s1-1 c) + sum-0)) + ;; s1 ^= s0 + (setf s1-1 (logxor s1-1 s0-1) + s1-0 (logxor s1-0 s0-0)) + ;; s[0] = rotl(s0,55) ^ s1 ^ (s1 << 14) + (multiple-value-setq (s0-1 s0-0) + (rotl-55 s0-1 s0-0)) + (setf s0-1 (logxor s0-1 s1-1) + s0-0 (logxor s0-0 s1-0)) + (multiple-value-bind (s14-1 s14-0) + (shl-14 s1-1 s1-0) + (setf s0-1 (logxor s0-1 s14-1) + s0-0 (logxor s0-0 s14-0))) + (setf (aref s 0) s0-0) + (setf (aref s 1) s0-1) + + (multiple-value-bind (r1 r0) + (rotl-36 s1-1 s1-0) + (setf (aref s 2) r0 + (aref s 3) r1)) + (setf (aref state 0) (make-double s0-1 s0-0) + (aref state 1) (make-double s1-1 s1-0)))))) + +;;; %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)) +(declaim (ftype (function ((single-float (0f0)) random-state) + (single-float 0f0)) + %random-single-float)) +;;; +(defun %random-single-float (arg state) + (declare (type (single-float (0f0)) arg) + (type random-state state)) + (* arg + (- (make-single-float + (dpb (ash (random-chunk state) + (- vm:single-float-digits random-chunk-length)) + vm:single-float-significand-byte + (single-float-bits 1.0))) + 1.0))) +;;; +(declaim (ftype (function ((double-float (0d0)) random-state) + (double-float 0d0)) + %random-double-float)) +;;; +;;; 53bit version. +;;; +#-x86 +(defun %random-double-float (arg state) + (declare (type (double-float (0d0)) arg) + (type random-state state)) + (* arg + (- (lisp::make-double-float + (dpb (ash (random-chunk state) + (- vm:double-float-digits random-chunk-length + vm:word-bits)) + vm:double-float-significand-byte + (lisp::double-float-high-bits 1d0)) + (random-chunk state)) + 1d0))) + +;;; Using a faster inline VOP. +#+x86 +(defun %random-double-float (arg state) + (declare (type (double-float (0d0)) arg) + (type random-state state)) + (let ((state-vector (random-state-state state))) + (* arg + (- (lisp::make-double-float + (dpb (ash (vm::random-mt19937 state-vector) + (- vm:double-float-digits random-chunk-length + vm:word-bits)) + vm:double-float-significand-byte + (lisp::double-float-high-bits 1d0)) + (vm::random-mt19937 state-vector)) + 1d0)))) + +#+long-float +(declaim (inline %random-long-float)) +#+long-float +(declaim (ftype (function ((long-float (0l0)) random-state) (long-float 0l0)) + %random-long-float)) + +;;; Using a faster inline VOP. +#+(and long-float x86) +(defun %random-long-float (arg state) + (declare (type (long-float (0l0)) arg) + (type random-state state)) + (let ((state-vector (random-state-state state))) + (* arg + (- (lisp::make-long-float + (lisp::long-float-exp-bits 1l0) + (logior (vm::random-mt19937 state-vector) vm:long-float-hidden-bit) + (vm::random-mt19937 state-vector)) + 1l0)))) + +#+(and long-float sparc) +(defun %random-long-float (arg state) + (declare (type (long-float (0l0)) arg) + (type random-state state)) + (* arg + (- (lisp::make-long-float + (lisp::long-float-exp-bits 1l0) ; X needs more work + (random-chunk state) (random-chunk state) (random-chunk state)) + 1l0))) +#+double-double +(defun %random-double-double-float (arg state) + (declare (type (double-double-float (0w0)) arg) + (type random-state state)) + ;; Generate a 31-bit integer, scale it and sum them up + (let* ((r 0w0) + (scale (scale-float 1d0 -31)) + (mult scale)) + (declare (double-float mult) + (type double-double-float r) + (optimize (speed 3) (inhibit-warnings 3))) + (dotimes (k 4) + (setf r (+ r (* mult (ldb (byte 31 0) (random-chunk state))))) + (setf mult (* mult scale))) + (* arg r))) + + +;;;; 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))))) + +(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." + (declare (inline %random-single-float %random-double-float + #+long-float %long-float)) + (cond + ((typep arg '(integer 1 #x100000000)) + ;; Let the compiler deftransform take care of this case. + (random arg state)) + ((and (typep arg 'single-float) (> arg 0.0F0)) + (%random-single-float arg state)) + ((and (typep arg 'double-float) (> arg 0.0D0)) + (%random-double-float arg state)) + #+long-float + ((and (typep arg 'long-float) (> arg 0.0L0)) + (%random-long-float arg state)) + #+double-double + ((and (typep arg 'double-double-float) (> arg 0.0w0)) + (%random-double-double-float arg state)) + ((and (integerp arg) (> arg 0)) + (%random-integer arg state)) + (t + (error 'simple-type-error + :expected-type '(or (integer 1) (float (0.0))) :datum arg + :format-control (intl:gettext "Argument is not a positive integer or a positive float: ~S") + :format-arguments (list arg))))) +
View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/commit/b119b34f82807a862a763e93e8...
--- View it on GitLab: https://gitlab.common-lisp.net/cmucl/cmucl/commit/b119b34f82807a862a763e93e8... You're receiving this email because of your account on gitlab.common-lisp.net.