Ignore:
Timestamp:
Dec 22, 2009, 12:00:06 AM (10 years ago)
Author:
rme
Message:

New random number generator %mrg31k3p.

Remove function init-random-state-seeds (moved to l1-numbers.lisp).

Update cl:random to use new generator.

In %bignum-random, use new inline function %16-random-bits.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/new-random/level-0/l0-numbers.lisp

    r13132 r13314  
    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 
     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#-x86-target
     1809(defun %mrg31k3p (state)
     1810  (let* ((v 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#-x86-target
     1848(defun %16-random-bits (state)
     1849  (%next-random-seed state))
     1850
     1851#+x86-target
     1852(defun %16-random-bits (state)
     1853  (logand #xffff (the fixnum (%mrg31kp3 state))))
     1854
     1855#+64-bit-target
     1856(defun %big-fixnum-random (number state)
     1857  (declare (fixnum number)
     1858           (ftype (function (random-state) fixnum) %mrg31k3p))
     1859  (let ((low (ldb (byte 30 0) (%mrg31k3p state)))
     1860        (high (ldb (byte 30 0) (%mrg31k3p state))))
     1861    (declare (fixnum low high))
     1862    (fast-mod (logior low (the fixnum (ash high 30)))
     1863              number)))
     1864
     1865#|
     1866Date: Mon, 3 Feb 1997 10:04:08 -0500
     1867To: info-mcl@digitool.com, wineberg@franz.scs.carleton.ca
     1868From: dds@flavors.com (Duncan Smith)
     1869Subject: Re: More info on the random number generator
     1870Sender: owner-info-mcl@digitool.com
     1871Precedence: bulk
     1872
     1873The generator is a Linear Congruential Generator:
     1874
     1875   X[n+1] = (aX[n] + c) mod m
     1876
     1877where: a = 16807  (Park&Miller recommend 48271)
     1878       c = 0
     1879       m = 2^31 - 1
     1880
     1881See: Knuth, Seminumerical Algorithms (Volume 2), Chapter 3.
     1882
     1883The period is: 2^31 - 2  (zero is excluded).
     1884
     1885What makes this generator so simple is that multiplication and addition mod
     18862^n-1 is easy.  See Knuth Ch. 4.3.2 (2nd Ed. p 272).
     1887
     1888    ab mod m = ...
     1889
     1890If         m = 2^n-1
     1891           u = ab mod 2^n
     1892           v = floor( ab / 2^n )
     1893
     1894    ab mod m = u + v                   :  u+v < 2^n
     1895    ab mod m = ((u + v) mod 2^n) + 1   :  u+v >= 2^n
     1896
     1897What we do is use 2b and 2^n so we can do arithemetic mod 2^32 instead of
     18982^31.  This reduces the whole generator to 5 instructions on the 680x0 or
     189980x86, and 8 on the 60x.
     1900
     1901-Duncan
     1902
     1903|#
     1904
     1905#+(and 32-bit-target (not x86-target))
     1906(defun %next-random-seed (state)
     1907  (multiple-value-bind (high low) (%next-random-pair (random.seed-1 state)
     1908                                                     (random.seed-2 state))
     1909    (declare (type (unsigned-byte 15) high)
     1910             (type (unsigned-byte 16) low))
     1911    (setf (random.seed-1 state) high
     1912          (random.seed-2 state) low)
     1913    (logior high (thes fixnum (logand low (ash 1 15))))))
     1914
     1915;;; When using a dead simple random number generator, it's reasonable
     1916;;; to take 16 bits of the output and discard the rest.  With a more
     1917;;; expensive generator, however, it may be worthwhile to do more bit
     1918;;; fiddling here here so that we can use all of the random bits
     1919;;; produced by %mrg31k2p.
    18011920#+32-bit-target
     1921(defun %bignum-random (number state)
     1922  (let* ((bits (+ (integer-length number) 8))
     1923         (half-words (ash (the fixnum (+ bits 15)) -4))
     1924         (long-words (ash (+ half-words 1) -1))
     1925         (dividend (%alloc-misc long-words target::subtag-bignum))
     1926         (16-bit-dividend dividend)
     1927         (index 1))
     1928    (declare (fixnum long-words index bits)
     1929             (dynamic-extent dividend)
     1930             (type (simple-array (unsigned-byte 16) (*)) 16-bit-dividend) ;lie
     1931             (optimize (speed 3) (safety 0)))
     1932    (loop
     1933       ;; This had better inline due to the lie above, or it will error
     1934       #+big-endian-target
     1935       (setf (aref 16-bit-dividend index) (%16-random-bits state))
     1936       #+little-endian-target
     1937       (setf (aref 16-bit-dividend (the fixnum (1- index)))
     1938             (%16-random-bits state))
     1939       (decf half-words)
     1940       (when (<= half-words 0) (return))
     1941       #+big-endian-target
     1942       (setf (aref 16-bit-dividend (the fixnum (1- index)))
     1943             (%16-random-bits state))
     1944       #+little-endian-target
     1945       (setf (aref 16-bit-dividend index) (%16-random-bits state))
     1946       (decf half-words)
     1947       (when (<= half-words 0) (return))
     1948       (incf index 2))
     1949    ;; The bignum code expects normalized bignums
     1950    (let* ((result (mod dividend number)))
     1951      (if (eq dividend result)
     1952        (copy-uvector result)
     1953        result))))
     1954
     1955(defun %float-random (number state)
     1956  (let ((ratio (gvector :ratio (random target::target-most-positive-fixnum state) target::target-most-positive-fixnum)))
     1957    (declare (dynamic-extent ratio))
     1958    (* number ratio)))
     1959
     1960#+(and 32-bit-target (not x86-target))
    18021961(defun random (number &optional (state *random-state*))
    18031962  (if (not (typep state 'random-state)) (report-bad-arg state 'random-state))
     
    18211980     (t (report-bad-arg number '(or (integer (0)) (float (0.0)))))))
    18221981
    1823 #+64-bit-target
     1982#+x86-target
    18241983(defun random (number &optional (state *random-state*))
    18251984  (if (not (typep state 'random-state)) (report-bad-arg state 'random-state))
    18261985  (cond
    18271986    ((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))))))))
     1987     #+32-bit-target
     1988     (fast-mod (%mrg31k3p state) number)
     1989     #+64-bit-target
     1990     (if (< number mrg31k3p-limit)
     1991       (fast-mod (%mrg31k3p state) number)
     1992       (%big-fixnum-random number state)))
    18351993    ((and (typep number 'double-float) (> (the double-float number) 0.0))
    18361994     (%float-random number state))
     
    18401998     (%bignum-random number state))
    18411999    (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 
    1901 #+32-bit-target
    1902 (defun %bignum-random (number state)
    1903   (let* ((bits (+ (integer-length number) 8))
    1904          (half-words (ash (the fixnum (+ bits 15)) -4))
    1905          (long-words (ash (+ half-words 1) -1))
    1906          (dividend (%alloc-misc long-words target::subtag-bignum))
    1907          (16-bit-dividend dividend)
    1908          (index 1))
    1909     (declare (fixnum long-words index bits)
    1910              (dynamic-extent dividend)
    1911              (type (simple-array (unsigned-byte 16) (*)) 16-bit-dividend) ;lie
    1912              (optimize (speed 3) (safety 0)))
    1913     (loop
    1914        ;; This had better inline due to the lie above, or it will error
    1915        #+big-endian-target
    1916        (setf (aref 16-bit-dividend index) (%next-random-seed state))
    1917        #+little-endian-target
    1918        (setf (aref 16-bit-dividend (the fixnum (1- index)))
    1919              (%next-random-seed state))
    1920        (decf half-words)
    1921        (when (<= half-words 0) (return))
    1922        #+big-endian-target
    1923        (setf (aref 16-bit-dividend (the fixnum (1- index)))
    1924              (%next-random-seed state))
    1925        #+little-endian-target
    1926        (setf (aref 16-bit-dividend index) (%next-random-seed state))
    1927        (decf half-words)
    1928        (when (<= half-words 0) (return))
    1929        (incf index 2))
    1930     ;; The bignum code expects normalized bignums
    1931     (let* ((result (mod dividend number)))
    1932       (if (eq dividend result)
    1933         (copy-uvector result)
    1934         result))))
    1935 
    1936 (defun %float-random (number state)
    1937   (let ((ratio (gvector :ratio (random target::target-most-positive-fixnum state) target::target-most-positive-fixnum)))
    1938     (declare (dynamic-extent ratio))
    1939     (* number ratio)))
    19402000
    19412001(eval-when (:compile-toplevel :execute)
Note: See TracChangeset for help on using the changeset viewer.