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-0/l0-numbers.lisp

    r13132 r13327  
    17831783  'real)
    17841784
    1785 
    1786 
    1787 (defun init-random-state-seeds ()
    1788   (let* ((ticks (ldb (byte 32 0) (+ (mixup-hash-code (%current-tcr))
    1789                                     (let* ((iface(primary-ip-interface)))
    1790                                       (or (and iface (ip-interface-addr iface))
    1791                                           0))
    1792                                     (mixup-hash-code
    1793                                      (logand (get-internal-real-time)
    1794                                              (1- target::target-most-positive-fixnum))))))
    1795          (high (ldb (byte 16 16) (if (zerop ticks) #x10000 ticks)))
    1796          (low (ldb (byte 16 0) ticks)))
    1797     (declare (fixnum high low))
    1798     (values high low)))
    1799 
    1800 
    1801 #+32-bit-target
    1802 (defun random (number &optional (state *random-state*))
    1803   (if (not (typep state 'random-state)) (report-bad-arg state 'random-state))
    1804   (cond
    1805      ((and (fixnump number) (> (the fixnum number) 0))
    1806       (locally (declare (fixnum number))
    1807         (if (< number 65536)
    1808           (fast-mod (%next-random-seed state) number)
    1809           (let* ((n 0)
    1810                  (nhalf (ash (+ 15 (integer-length number)) -4)))
    1811             (declare (fixnum n nhalf))
    1812             (dotimes (i nhalf (fast-mod n number))
    1813               (setq n (logior (the fixnum (ash n 16))
    1814                               (the fixnum (%next-random-seed state)))))))))
    1815      ((and (typep number 'double-float) (> (the double-float number) 0.0))
    1816       (%float-random number state))
    1817      ((and (typep number 'short-float) (> (the short-float number) 0.0s0))
    1818       (%float-random number state))
    1819      ((and (bignump number) (> number 0))
    1820       (%bignum-random number state))
    1821      (t (report-bad-arg number '(or (integer (0)) (float (0.0)))))))
     1785;;; This is the MRG31k3p random number generator described in
     1786;;; P. L'Ecuyer and R. Touzin, "Fast Combined Multiple Recursive
     1787;;; Generators with Multipliers of the form a = +/- 2^d +/- 2^e"",
     1788;;; Proceedings of the 2000 Winter Simulation Conference, Dec. 2000,
     1789;;; 683--689.
     1790;;;
     1791;;; A link to the paper is available on L'Ecuyer's web site:
     1792;;; http://www.iro.umontreal.ca/~lecuyer/papers.html.
     1793;;;
     1794;;; This generator has a period of about 2^185.  It produces values in
     1795;;; in the half-open interval [0, 2^31 - 1).
     1796;;;
     1797;;; It uses 6 words of state.
     1798
     1799(defconstant mrg31k3p-m1 #.(- (expt 2 31) 1))
     1800(defconstant mrg31k3p-m2 #.(- (expt 2 31) 21069))
     1801(defconstant mrg31k3p-limit #.(1- (expt 2 31))
     1802             "Exclusive upper bound on values returned by %mrg31k3p.")
     1803
     1804
     1805;;; This is a portable version of the MRG31k3p generator.  It's not
     1806;;; too bad in a 64-bit CCL, but the generator pretty much has to be
     1807;;; in LAP for 32-bit ports.
     1808#-(or x8632-target ppc32-target x8664-target)
     1809(defun %mrg31k3p (state)
     1810  (let* ((v (random.mrg31k3p-state state)))
     1811    (declare (type (simple-array (unsigned-byte 32) (*)) v)
     1812             (optimize speed))
     1813    (let ((y1 (+ (+ (ash (logand (aref v 1) #x1ff) 22)
     1814                    (ash (aref v 1) -9))
     1815                 (+ (ash (logand (aref v 2) #xffffff) 7)
     1816                    (ash (aref v 2) -24)))))
     1817      (declare (type (unsigned-byte 32) y1))
     1818      (if (>= y1 mrg31k3p-m1) (decf y1 mrg31k3p-m1))
     1819      (incf y1 (aref v 2))
     1820      (if (>= y1 mrg31k3p-m1) (decf y1 mrg31k3p-m1))
     1821      (setf (aref v 2) (aref v 1)
     1822            (aref v 1) (aref v 0)
     1823            (aref v 0) y1))
     1824    (let ((y1 (+ (ash (logand (aref v 3) #xffff) 15)
     1825                 (* 21069 (ash (aref v 3) -16))))
     1826          (y2 (+ (ash (logand (aref v 5) #xffff) 15)
     1827                 (* 21069 (ash (aref v 5) -16)))))
     1828      (declare (type (unsigned-byte 32) y1 y2))
     1829      (if (>= y1 mrg31k3p-m2) (decf y1 mrg31k3p-m2))
     1830      (if (>= y2 mrg31k3p-m2) (decf y2 mrg31k3p-m2))
     1831      (incf y2 (aref v 5))
     1832      (if (>= y2 mrg31k3p-m2) (decf y2 mrg31k3p-m2))
     1833      (incf y2 y1)
     1834      (if (>= y2 mrg31k3p-m2) (decf y2 mrg31k3p-m2))
     1835      (setf (aref v 5) (aref v 4)
     1836            (aref v 4) (aref v 3)
     1837            (aref v 3) y2))
     1838    (let* ((x10 (aref v 0))
     1839           (x20 (aref v 3)))
     1840      (if (<= x10 x20)
     1841        (+ (- x10 x20) mrg31k3p-m1)
     1842        (- x10 x20)))))
     1843
     1844(eval-when (:compile-toplevel :execute)
     1845  (declaim (inline %16-random-bits)))
     1846
     1847(defun %16-random-bits (state)
     1848  (logand #xffff (the fixnum (%mrg31k3p state))))
    18221849
    18231850#+64-bit-target
    1824 (defun random (number &optional (state *random-state*))
    1825   (if (not (typep state 'random-state)) (report-bad-arg state 'random-state))
    1826   (cond
    1827     ((and (fixnump number) (> (the fixnum number) 0))
    1828      (locally (declare (fixnum number))
    1829        (let* ((n 0)
    1830               (n32 (ash (+ 31 (integer-length number)) -5)))
    1831          (declare (fixnum n n32))
    1832          (dotimes (i n32 (fast-mod n number))
    1833            (setq n (logior (the fixnum (ash n 32))
    1834                            (the fixnum (%next-random-seed state))))))))
    1835     ((and (typep number 'double-float) (> (the double-float number) 0.0))
    1836      (%float-random number state))
    1837     ((and (typep number 'short-float) (> (the short-float number) 0.0s0))
    1838      (%float-random number state))
    1839     ((and (bignump number) (> number 0))
    1840      (%bignum-random number state))
    1841     (t (report-bad-arg number '(or (integer (0)) (float (0.0)))))))
    1842 
    1843 
    1844 #|
    1845 Date: Mon, 3 Feb 1997 10:04:08 -0500
    1846 To: info-mcl@digitool.com, wineberg@franz.scs.carleton.ca
    1847 From: dds@flavors.com (Duncan Smith)
    1848 Subject: Re: More info on the random number generator
    1849 Sender: owner-info-mcl@digitool.com
    1850 Precedence: bulk
    1851 
    1852 The generator is a Linear Congruential Generator:
    1853 
    1854    X[n+1] = (aX[n] + c) mod m
    1855 
    1856 where: a = 16807  (Park&Miller recommend 48271)
    1857        c = 0
    1858        m = 2^31 - 1
    1859 
    1860 See: Knuth, Seminumerical Algorithms (Volume 2), Chapter 3.
    1861 
    1862 The period is: 2^31 - 2  (zero is excluded).
    1863 
    1864 What makes this generator so simple is that multiplication and addition mod
    1865 2^n-1 is easy.  See Knuth Ch. 4.3.2 (2nd Ed. p 272).
    1866 
    1867     ab mod m = ...
    1868 
    1869 If         m = 2^n-1
    1870            u = ab mod 2^n
    1871            v = floor( ab / 2^n )
    1872 
    1873     ab mod m = u + v                   :  u+v < 2^n
    1874     ab mod m = ((u + v) mod 2^n) + 1   :  u+v >= 2^n
    1875 
    1876 What we do is use 2b and 2^n so we can do arithemetic mod 2^32 instead of
    1877 2^31.  This reduces the whole generator to 5 instructions on the 680x0 or
    1878 80x86, and 8 on the 60x.
    1879 
    1880 -Duncan
    1881 
    1882 |#
    1883 
    1884 #+64-bit-target
    1885 (defun %next-random-seed (state)
    1886   (let* ((n (the fixnum (* (the fixnum (random.seed-1 state)) 48271))))
    1887     (declare (fixnum n))
    1888     (setf (random.seed-1 state) (fast-mod n (1- (expt 2 31))))
    1889     (logand n (1- (ash 1 32)))))
    1890 
    1891 #+32-bit-target
    1892 (defun %next-random-seed (state)
    1893   (multiple-value-bind (high low) (%next-random-pair (random.seed-1 state)
    1894                                                      (random.seed-2 state))
    1895     (declare (type (unsigned-byte 15) high)
    1896              (type (unsigned-byte 16) low))
    1897     (setf (random.seed-1 state) high
    1898           (random.seed-2 state) low)
    1899     (logior high (the fixnum (logand low (ash 1 15))))))
    1900 
     1851(defun %big-fixnum-random (number state)
     1852  (declare (fixnum number)
     1853           (ftype (function (random-state) fixnum) %mrg31k3p))
     1854  (let ((low (ldb (byte 30 0) (%mrg31k3p state)))
     1855        (high (ldb (byte 30 0) (%mrg31k3p state))))
     1856    (declare (fixnum low high))
     1857    (fast-mod (logior low (the fixnum (ash high 30)))
     1858              number)))
     1859
     1860;;; When using a dead simple random number generator, it's reasonable
     1861;;; to take 16 bits of the output and discard the rest.  With a more
     1862;;; expensive generator, however, it may be worthwhile to do more bit
     1863;;; fiddling here here so that we can use all of the random bits
     1864;;; produced by %mrg31k2p.
    19011865#+32-bit-target
    19021866(defun %bignum-random (number state)
     
    19141878       ;; This had better inline due to the lie above, or it will error
    19151879       #+big-endian-target
    1916        (setf (aref 16-bit-dividend index) (%next-random-seed state))
     1880       (setf (aref 16-bit-dividend index) (%16-random-bits state))
    19171881       #+little-endian-target
    19181882       (setf (aref 16-bit-dividend (the fixnum (1- index)))
    1919              (%next-random-seed state))
     1883             (%16-random-bits state))
    19201884       (decf half-words)
    19211885       (when (<= half-words 0) (return))
    19221886       #+big-endian-target
    19231887       (setf (aref 16-bit-dividend (the fixnum (1- index)))
    1924              (%next-random-seed state))
     1888             (%16-random-bits state))
    19251889       #+little-endian-target
    1926        (setf (aref 16-bit-dividend index) (%next-random-seed state))
     1890       (setf (aref 16-bit-dividend index) (%16-random-bits state))
    19271891       (decf half-words)
    19281892       (when (<= half-words 0) (return))
     
    19381902    (declare (dynamic-extent ratio))
    19391903    (* number ratio)))
     1904
     1905(defun random (number &optional (state *random-state*))
     1906  (if (not (typep state 'random-state)) (report-bad-arg state 'random-state))
     1907  (cond
     1908    ((and (fixnump number) (> (the fixnum number) 0))
     1909     #+32-bit-target
     1910     (fast-mod (%mrg31k3p state) number)
     1911     #+64-bit-target
     1912     (if (< number mrg31k3p-limit)
     1913       (fast-mod (%mrg31k3p state) number)
     1914       (%big-fixnum-random number state)))
     1915    ((and (typep number 'double-float) (> (the double-float number) 0.0))
     1916     (%float-random number state))
     1917    ((and (typep number 'short-float) (> (the short-float number) 0.0s0))
     1918     (%float-random number state))
     1919    ((and (bignump number) (> number 0))
     1920     (%bignum-random number state))
     1921    (t (report-bad-arg number '(or (integer (0)) (float (0.0)))))))
    19401922
    19411923(eval-when (:compile-toplevel :execute)
Note: See TracChangeset for help on using the changeset viewer.