Changeset 13314 for branches/newrandom/level0
 Timestamp:
 Dec 22, 2009, 12:00:06 AM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/newrandom/level0/l0numbers.lisp
r13132 r13314 1783 1783 'real) 1784 1784 1785 1786 1787 (defun initrandomstateseeds () 1788 (let* ((ticks (ldb (byte 32 0) (+ (mixuphashcode (%currenttcr)) 1789 (let* ((iface(primaryipinterface))) 1790 (or (and iface (ipinterfaceaddr iface)) 1791 0)) 1792 (mixuphashcode 1793 (logand (getinternalrealtime) 1794 (1 target::targetmostpositivefixnum)))))) 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 ;;; 683689. 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 halfopen interval [0, 2^31  1). 1796 ;;; 1797 ;;; It uses 6 words of state. 1798 1799 (defconstant mrg31k3pm1 #.( (expt 2 31) 1)) 1800 (defconstant mrg31k3pm2 #.( (expt 2 31) 21069)) 1801 (defconstant mrg31k3plimit #.(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 64bit CCL, but the generator pretty much has to be 1807 ;;; in LAP for 32bit ports. 1808 #x86target 1809 (defun %mrg31k3p (state) 1810 (let* ((v state)) 1811 (declare (type (simplearray (unsignedbyte 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 (unsignedbyte 32) y1)) 1818 (if (>= y1 mrg31k3pm1) (decf y1 mrg31k3pm1)) 1819 (incf y1 (aref v 2)) 1820 (if (>= y1 mrg31k3pm1) (decf y1 mrg31k3pm1)) 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 (unsignedbyte 32) y1 y2)) 1829 (if (>= y1 mrg31k3pm2) (decf y1 mrg31k3pm2)) 1830 (if (>= y2 mrg31k3pm2) (decf y2 mrg31k3pm2)) 1831 (incf y2 (aref v 5)) 1832 (if (>= y2 mrg31k3pm2) (decf y2 mrg31k3pm2)) 1833 (incf y2 y1) 1834 (if (>= y2 mrg31k3pm2) (decf y2 mrg31k3pm2)) 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) mrg31k3pm1) 1842 ( x10 x20))))) 1843 1844 (evalwhen (:compiletoplevel :execute) 1845 (declaim (inline %16randombits))) 1846 1847 #x86target 1848 (defun %16randombits (state) 1849 (%nextrandomseed state)) 1850 1851 #+x86target 1852 (defun %16randombits (state) 1853 (logand #xffff (the fixnum (%mrg31kp3 state)))) 1854 1855 #+64bittarget 1856 (defun %bigfixnumrandom (number state) 1857 (declare (fixnum number) 1858 (ftype (function (randomstate) 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 (fastmod (logior low (the fixnum (ash high 30))) 1863 number))) 1864 1865 # 1866 Date: Mon, 3 Feb 1997 10:04:08 0500 1867 To: infomcl@digitool.com, wineberg@franz.scs.carleton.ca 1868 From: dds@flavors.com (Duncan Smith) 1869 Subject: Re: More info on the random number generator 1870 Sender: ownerinfomcl@digitool.com 1871 Precedence: bulk 1872 1873 The generator is a Linear Congruential Generator: 1874 1875 X[n+1] = (aX[n] + c) mod m 1876 1877 where: a = 16807 (Park&Miller recommend 48271) 1878 c = 0 1879 m = 2^31  1 1880 1881 See: Knuth, Seminumerical Algorithms (Volume 2), Chapter 3. 1882 1883 The period is: 2^31  2 (zero is excluded). 1884 1885 What makes this generator so simple is that multiplication and addition mod 1886 2^n1 is easy. See Knuth Ch. 4.3.2 (2nd Ed. p 272). 1887 1888 ab mod m = ... 1889 1890 If m = 2^n1 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 1897 What we do is use 2b and 2^n so we can do arithemetic mod 2^32 instead of 1898 2^31. This reduces the whole generator to 5 instructions on the 680x0 or 1899 80x86, and 8 on the 60x. 1900 1901 Duncan 1902 1903 # 1904 1905 #+(and 32bittarget (not x86target)) 1906 (defun %nextrandomseed (state) 1907 (multiplevaluebind (high low) (%nextrandompair (random.seed1 state) 1908 (random.seed2 state)) 1909 (declare (type (unsignedbyte 15) high) 1910 (type (unsignedbyte 16) low)) 1911 (setf (random.seed1 state) high 1912 (random.seed2 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. 1801 1920 #+32bittarget 1921 (defun %bignumrandom (number state) 1922 (let* ((bits (+ (integerlength number) 8)) 1923 (halfwords (ash (the fixnum (+ bits 15)) 4)) 1924 (longwords (ash (+ halfwords 1) 1)) 1925 (dividend (%allocmisc longwords target::subtagbignum)) 1926 (16bitdividend dividend) 1927 (index 1)) 1928 (declare (fixnum longwords index bits) 1929 (dynamicextent dividend) 1930 (type (simplearray (unsignedbyte 16) (*)) 16bitdividend) ;lie 1931 (optimize (speed 3) (safety 0))) 1932 (loop 1933 ;; This had better inline due to the lie above, or it will error 1934 #+bigendiantarget 1935 (setf (aref 16bitdividend index) (%16randombits state)) 1936 #+littleendiantarget 1937 (setf (aref 16bitdividend (the fixnum (1 index))) 1938 (%16randombits state)) 1939 (decf halfwords) 1940 (when (<= halfwords 0) (return)) 1941 #+bigendiantarget 1942 (setf (aref 16bitdividend (the fixnum (1 index))) 1943 (%16randombits state)) 1944 #+littleendiantarget 1945 (setf (aref 16bitdividend index) (%16randombits state)) 1946 (decf halfwords) 1947 (when (<= halfwords 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 (copyuvector result) 1953 result)))) 1954 1955 (defun %floatrandom (number state) 1956 (let ((ratio (gvector :ratio (random target::targetmostpositivefixnum state) target::targetmostpositivefixnum))) 1957 (declare (dynamicextent ratio)) 1958 (* number ratio))) 1959 1960 #+(and 32bittarget (not x86target)) 1802 1961 (defun random (number &optional (state *randomstate*)) 1803 1962 (if (not (typep state 'randomstate)) (reportbadarg state 'randomstate)) … … 1821 1980 (t (reportbadarg number '(or (integer (0)) (float (0.0))))))) 1822 1981 1823 #+ 64bittarget1982 #+x86target 1824 1983 (defun random (number &optional (state *randomstate*)) 1825 1984 (if (not (typep state 'randomstate)) (reportbadarg state 'randomstate)) 1826 1985 (cond 1827 1986 ((and (fixnump number) (> (the fixnum number) 0)) 1828 (locally (declare (fixnum number)) 1829 (let* ((n 0) 1830 (n32 (ash (+ 31 (integerlength number)) 5))) 1831 (declare (fixnum n n32)) 1832 (dotimes (i n32 (fastmod n number)) 1833 (setq n (logior (the fixnum (ash n 32)) 1834 (the fixnum (%nextrandomseed state)))))))) 1987 #+32bittarget 1988 (fastmod (%mrg31k3p state) number) 1989 #+64bittarget 1990 (if (< number mrg31k3plimit) 1991 (fastmod (%mrg31k3p state) number) 1992 (%bigfixnumrandom number state))) 1835 1993 ((and (typep number 'doublefloat) (> (the doublefloat number) 0.0)) 1836 1994 (%floatrandom number state)) … … 1840 1998 (%bignumrandom number state)) 1841 1999 (t (reportbadarg number '(or (integer (0)) (float (0.0))))))) 1842 1843 1844 #1845 Date: Mon, 3 Feb 1997 10:04:08 05001846 To: infomcl@digitool.com, wineberg@franz.scs.carleton.ca1847 From: dds@flavors.com (Duncan Smith)1848 Subject: Re: More info on the random number generator1849 Sender: ownerinfomcl@digitool.com1850 Precedence: bulk1851 1852 The generator is a Linear Congruential Generator:1853 1854 X[n+1] = (aX[n] + c) mod m1855 1856 where: a = 16807 (Park&Miller recommend 48271)1857 c = 01858 m = 2^31  11859 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 mod1865 2^n1 is easy. See Knuth Ch. 4.3.2 (2nd Ed. p 272).1866 1867 ab mod m = ...1868 1869 If m = 2^n11870 u = ab mod 2^n1871 v = floor( ab / 2^n )1872 1873 ab mod m = u + v : u+v < 2^n1874 ab mod m = ((u + v) mod 2^n) + 1 : u+v >= 2^n1875 1876 What we do is use 2b and 2^n so we can do arithemetic mod 2^32 instead of1877 2^31. This reduces the whole generator to 5 instructions on the 680x0 or1878 80x86, and 8 on the 60x.1879 1880 Duncan1881 1882 #1883 1884 #+64bittarget1885 (defun %nextrandomseed (state)1886 (let* ((n (the fixnum (* (the fixnum (random.seed1 state)) 48271))))1887 (declare (fixnum n))1888 (setf (random.seed1 state) (fastmod n (1 (expt 2 31))))1889 (logand n (1 (ash 1 32)))))1890 1891 #+32bittarget1892 (defun %nextrandomseed (state)1893 (multiplevaluebind (high low) (%nextrandompair (random.seed1 state)1894 (random.seed2 state))1895 (declare (type (unsignedbyte 15) high)1896 (type (unsignedbyte 16) low))1897 (setf (random.seed1 state) high1898 (random.seed2 state) low)1899 (logior high (the fixnum (logand low (ash 1 15))))))1900 1901 #+32bittarget1902 (defun %bignumrandom (number state)1903 (let* ((bits (+ (integerlength number) 8))1904 (halfwords (ash (the fixnum (+ bits 15)) 4))1905 (longwords (ash (+ halfwords 1) 1))1906 (dividend (%allocmisc longwords target::subtagbignum))1907 (16bitdividend dividend)1908 (index 1))1909 (declare (fixnum longwords index bits)1910 (dynamicextent dividend)1911 (type (simplearray (unsignedbyte 16) (*)) 16bitdividend) ;lie1912 (optimize (speed 3) (safety 0)))1913 (loop1914 ;; This had better inline due to the lie above, or it will error1915 #+bigendiantarget1916 (setf (aref 16bitdividend index) (%nextrandomseed state))1917 #+littleendiantarget1918 (setf (aref 16bitdividend (the fixnum (1 index)))1919 (%nextrandomseed state))1920 (decf halfwords)1921 (when (<= halfwords 0) (return))1922 #+bigendiantarget1923 (setf (aref 16bitdividend (the fixnum (1 index)))1924 (%nextrandomseed state))1925 #+littleendiantarget1926 (setf (aref 16bitdividend index) (%nextrandomseed state))1927 (decf halfwords)1928 (when (<= halfwords 0) (return))1929 (incf index 2))1930 ;; The bignum code expects normalized bignums1931 (let* ((result (mod dividend number)))1932 (if (eq dividend result)1933 (copyuvector result)1934 result))))1935 1936 (defun %floatrandom (number state)1937 (let ((ratio (gvector :ratio (random target::targetmostpositivefixnum state) target::targetmostpositivefixnum)))1938 (declare (dynamicextent ratio))1939 (* number ratio)))1940 2000 1941 2001 (evalwhen (:compiletoplevel :execute)
Note: See TracChangeset
for help on using the changeset viewer.