Initial implementation of xoroshiro rng

;;; -*- 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)))))