Changeset 13327 for trunk/source/level0/l0numbers.lisp
 Timestamp:
 Dec 22, 2009, 7:10:27 AM (10 years ago)
 Location:
 trunk/source
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/source
 Property svn:mergeinfo changed
/branches/newrandom (added) merged: 1331013326
 Property svn:mergeinfo changed

trunk/source/level0/l0numbers.lisp
r13132 r13327 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 1801 #+32bittarget 1802 (defun random (number &optional (state *randomstate*)) 1803 (if (not (typep state 'randomstate)) (reportbadarg state 'randomstate)) 1804 (cond 1805 ((and (fixnump number) (> (the fixnum number) 0)) 1806 (locally (declare (fixnum number)) 1807 (if (< number 65536) 1808 (fastmod (%nextrandomseed state) number) 1809 (let* ((n 0) 1810 (nhalf (ash (+ 15 (integerlength number)) 4))) 1811 (declare (fixnum n nhalf)) 1812 (dotimes (i nhalf (fastmod n number)) 1813 (setq n (logior (the fixnum (ash n 16)) 1814 (the fixnum (%nextrandomseed state))))))))) 1815 ((and (typep number 'doublefloat) (> (the doublefloat number) 0.0)) 1816 (%floatrandom number state)) 1817 ((and (typep number 'shortfloat) (> (the shortfloat number) 0.0s0)) 1818 (%floatrandom number state)) 1819 ((and (bignump number) (> number 0)) 1820 (%bignumrandom number state)) 1821 (t (reportbadarg 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 ;;; 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 #(or x8632target ppc32target x8664target) 1809 (defun %mrg31k3p (state) 1810 (let* ((v (random.mrg31k3pstate 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 (defun %16randombits (state) 1848 (logand #xffff (the fixnum (%mrg31k3p state)))) 1822 1849 1823 1850 #+64bittarget 1824 (defun random (number &optional (state *randomstate*)) 1825 (if (not (typep state 'randomstate)) (reportbadarg state 'randomstate)) 1826 (cond 1827 ((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)))))))) 1835 ((and (typep number 'doublefloat) (> (the doublefloat number) 0.0)) 1836 (%floatrandom number state)) 1837 ((and (typep number 'shortfloat) (> (the shortfloat number) 0.0s0)) 1838 (%floatrandom number state)) 1839 ((and (bignump number) (> number 0)) 1840 (%bignumrandom number state)) 1841 (t (reportbadarg number '(or (integer (0)) (float (0.0))))))) 1842 1843 1844 # 1845 Date: Mon, 3 Feb 1997 10:04:08 0500 1846 To: infomcl@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: ownerinfomcl@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^n1 is easy. See Knuth Ch. 4.3.2 (2nd Ed. p 272). 1866 1867 ab mod m = ... 1868 1869 If m = 2^n1 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 #+64bittarget 1885 (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 #+32bittarget 1892 (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) high 1898 (random.seed2 state) low) 1899 (logior high (the fixnum (logand low (ash 1 15)))))) 1900 1851 (defun %bigfixnumrandom (number state) 1852 (declare (fixnum number) 1853 (ftype (function (randomstate) 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 (fastmod (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. 1901 1865 #+32bittarget 1902 1866 (defun %bignumrandom (number state) … … 1914 1878 ;; This had better inline due to the lie above, or it will error 1915 1879 #+bigendiantarget 1916 (setf (aref 16bitdividend index) (% nextrandomseedstate))1880 (setf (aref 16bitdividend index) (%16randombits state)) 1917 1881 #+littleendiantarget 1918 1882 (setf (aref 16bitdividend (the fixnum (1 index))) 1919 (% nextrandomseedstate))1883 (%16randombits state)) 1920 1884 (decf halfwords) 1921 1885 (when (<= halfwords 0) (return)) 1922 1886 #+bigendiantarget 1923 1887 (setf (aref 16bitdividend (the fixnum (1 index))) 1924 (% nextrandomseedstate))1888 (%16randombits state)) 1925 1889 #+littleendiantarget 1926 (setf (aref 16bitdividend index) (% nextrandomseedstate))1890 (setf (aref 16bitdividend index) (%16randombits state)) 1927 1891 (decf halfwords) 1928 1892 (when (<= halfwords 0) (return)) … … 1938 1902 (declare (dynamicextent ratio)) 1939 1903 (* number ratio))) 1904 1905 (defun random (number &optional (state *randomstate*)) 1906 (if (not (typep state 'randomstate)) (reportbadarg state 'randomstate)) 1907 (cond 1908 ((and (fixnump number) (> (the fixnum number) 0)) 1909 #+32bittarget 1910 (fastmod (%mrg31k3p state) number) 1911 #+64bittarget 1912 (if (< number mrg31k3plimit) 1913 (fastmod (%mrg31k3p state) number) 1914 (%bigfixnumrandom number state))) 1915 ((and (typep number 'doublefloat) (> (the doublefloat number) 0.0)) 1916 (%floatrandom number state)) 1917 ((and (typep number 'shortfloat) (> (the shortfloat number) 0.0s0)) 1918 (%floatrandom number state)) 1919 ((and (bignump number) (> number 0)) 1920 (%bignumrandom number state)) 1921 (t (reportbadarg number '(or (integer (0)) (float (0.0))))))) 1940 1922 1941 1923 (evalwhen (:compiletoplevel :execute)
Note: See TracChangeset
for help on using the changeset viewer.