Ignore:
Timestamp:
Dec 22, 2009, 7:10:27 AM (10 years ago)
Author:
rme
Message:

Improve CL:RANDOM.

The new generator is the MRG321k3p generator described in

  1. L'Ecuyer and R. Touzin, "Fast Combined Multiple Recursive

Generators with Multipliers of the form a = +/- 2d +/- 2e"",
Proceedings of the 2000 Winter Simulation Conference, Dec. 2000,
683--689.

It has a period of about 2185 and produces output of much higher
statistical quality than the previous generator.

Performance of the new generator should generally be comparable to that
of the old generator, despite the fact that the new generator does
a lot more work.

Location:
trunk/source
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/source

  • trunk/source/level-1/l1-numbers.lisp

    r13067 r13327  
    422422      nil)))
    423423
    424 (defun %cons-random-state (seed-1 seed-2)
    425   #+32-bit-target
    426   (%istruct 'random-state seed-1 seed-2)
    427   #+64-bit-target
    428   (%istruct 'random-state (the fixnum (+ (the fixnum seed-2)
    429                           (the fixnum (ash (the fixnum seed-1) 16))))))
    430 
    431424;;; random associated stuff except for the print-object method which
    432425;;; is still in "lib;numbers.lisp"
    433 (defun initialize-random-state (seed-1 seed-2)
    434   (unless (and (fixnump seed-1) (%i<= 0 seed-1) (%i< seed-1 #x10000))
    435     (report-bad-arg seed-1 '(unsigned-byte 16)))
    436   (unless (and (fixnump seed-2) (%i<= 0 seed-2) (%i< seed-2 #x10000))
    437     (report-bad-arg seed-2 '(unsigned-byte 16)))
    438     (%cons-random-state seed-1 seed-2))
     426
     427(defun init-random-state-seeds ()
     428  (let* ((ticks (ldb (byte 32 0)
     429                     (+ (mixup-hash-code (%current-tcr))
     430                        (let* ((iface (primary-ip-interface)))
     431                          (or (and iface (ip-interface-addr iface))
     432                              0))
     433                        (mixup-hash-code
     434                         (logand (get-internal-real-time)
     435                                 (1- target::target-most-positive-fixnum))))))
     436         (high (ldb (byte 16 16) (if (zerop ticks) #x10000 ticks)))
     437         (low (ldb (byte 16 0) ticks)))
     438    (declare (fixnum high low))
     439    (values high low)))
     440
     441(defun %cons-mrg31k3p-state (x0 x1 x2 x3 x4 x5)
     442  (let ((array (make-array 6 :element-type '(unsigned-byte 32)
     443                           :initial-contents (list x0 x1 x2 x3 x4 x5))))
     444    (%istruct 'random-state array)))
     445
     446(defun initialize-mrg31k3p-state (x0 x1 x2 x3 x4 x5)
     447  (let ((args (list x0 x1 x2 x3 x4 x5)))
     448    (declare (dynamic-extent args))
     449    (dolist (a args)
     450      (unless (and (fixnump a) (%i<= 0 a) (< a mrg31k3p-limit))
     451        (report-bad-arg a `(integer 0 (,mrg31k3p-limit)))))
     452    (when (and (zerop x0) (zerop x1) (zerop x2))
     453      (error "The first three arguments must not all be zero."))
     454    (when (and (zerop x3) (zerop x4) (zerop x5))
     455      (error "The second three arguments must not all be zero."))
     456    (%cons-mrg31k3p-state x0 x1 x2 x3 x4 x5)))
     457
     458(defun random-mrg31k3p-state ()
     459  (loop repeat 6
     460        for n = (init-random-state-seeds)
     461        ;; The first three seed elements must not be all zero, and
     462        ;; likewise for the second three.  Avoid the issue by
     463        ;; excluding zero values.
     464        collect (1+ (mod n (1- mrg31k3p-limit))) into seed
     465        finally (return (apply #'%cons-mrg31k3p-state seed))))
     466
     467(defun initial-random-state ()
     468  (initialize-mrg31k3p-state 12345 12345 12345 12345 12345 12345))
    439469
    440470(defun make-random-state (&optional state)
    441   "Make a random state object. If STATE is not supplied, return a copy
    442   of the default random state. If STATE is a random state, then return a
    443   copy of it. If STATE is T then return a random state generated from
    444   the universal time."
    445   (let* ((seed-1 0)
    446          (seed-2 0))
    447     (if (eq state t)
    448       (multiple-value-setq (seed-1 seed-2) (init-random-state-seeds))
    449       (progn
    450         (setq state (require-type (or state *random-state*) 'random-state))
    451         #+32-bit-target
    452         (setq seed-1 (random.seed-1 state) seed-2 (random.seed-2 state))
    453         #+64-bit-target
    454         (let* ((seed (random.seed-1 state)))
    455           (declare (type (unsigned-byte 32) seed))
    456           (setq seed-1 (ldb (byte 16 16) seed)
    457                 seed-2 (ldb (byte 16 0) seed)))))
    458     (%cons-random-state seed-1 seed-2)))
     471  "Make a new random state object. If STATE is not supplied, return a
     472  copy of the current random state. If STATE is a random state, then
     473  return a copy of it. If STATE is T then return a randomly
     474  initialized random state."
     475  (if (eq state t)
     476    (random-mrg31k3p-state)
     477    (progn
     478      (setq state (require-type (or state *random-state*) 'random-state))
     479      (let ((seed (coerce (random.mrg31k3p-state state) 'list)))
     480        (apply #'%cons-mrg31k3p-state seed)))))
    459481
    460482(defun random-state-p (thing) (istruct-typep thing 'random-state))
     483
     484(defun %random-state-equalp (x y)
     485  ;; x and y are both random-state objects
     486  (equalp (random.mrg31k3p-state x) (random.mrg31k3p-state y)))
    461487
    462488;;; transcendental stuff.  Should go in level-0;l0-float
Note: See TracChangeset for help on using the changeset viewer.