Changeset 13327


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:
12 edited

Legend:

Unmodified
Added
Removed
  • trunk/source

  • trunk/source/compiler/PPC/ppc-lapmacros.lisp

    r13067 r13327  
    10741074    (and. ,mask ,was ,mask)))
    10751075                                           
     1076(defppclapmacro u32-ref (dest index vector)
     1077  `(lwz ,dest (+ (* 4 ,index) target::misc-data-offset) ,vector))
     1078
     1079(defppclapmacro u32-set (new-value index vector)
     1080  `(stw ,new-value (+ (* 4 ,index) target::misc-data-offset) ,vector))
    10761081
    10771082(provide "PPC-LAPMACROS")
  • trunk/source/level-0/PPC/ppc-numbers.lisp

    r13067 r13327  
    438438    (beq cr0 @shiftv)
    439439    (b @loop)))
     440
     441(defppclapfunction %mrg31k3p ((state arg_z))
     442  (let ((seed temp0))
     443    (svref seed 1 state)
     444    (u32-ref imm0 1 seed)
     445    (u32-ref imm3 2 seed)
     446    (rlwinm imm1 imm0 22 1 9)
     447    (srwi imm2 imm0 9)
     448    (add imm0 imm1 imm2)
    440449   
    441 
    442 
     450    ;; construct m1 (1- (expt 2 31))
     451    (lis imm1 #x7fff)
     452    (ori imm1 imm1 #xffff)
     453
     454    (rlwinm imm4 imm3 7 1 24)
     455    (srwi imm5 imm3 24)
     456    (add imm0 imm0 imm4)
     457    (add imm0 imm0 imm5)
     458
     459    ;; reduce mod m1
     460    (cmplw cr7 imm0 imm1)
     461    (blt cr7 @ok1)
     462    (sub imm0 imm0 imm1)
     463    @ok1
     464
     465    (add imm0 imm0 imm3)
     466
     467    ;; reduce mod m1
     468    (cmplw cr7 imm0 imm1)
     469    (blt cr7 @ok2)
     470    (sub imm0 imm0 imm1)
     471    @ok2
     472
     473    ;; update state
     474    (u32-ref imm1 1 seed)
     475    (u32-set imm1 2 seed)
     476    (u32-ref imm1 0 seed)
     477    (u32-set imm1 1 seed)
     478    (u32-set imm0 0 seed)
     479
     480    ;; construct m2 (- (expt 2 31) 21069))
     481    (lis imm5 #x7fff)
     482    (ori imm5 imm5 44467)
     483
     484    ;; second component
     485    (u32-ref imm0 3 seed)
     486    (rlwinm imm1 imm0 15 1 16)
     487    (srwi imm2 imm0 16)
     488    (mulli imm2 imm2 21069)
     489    (add imm0 imm1 imm2)
     490
     491    ;; reduce mod m2
     492    (cmplw cr7 imm0 imm5)
     493    (blt cr7 @ok3)
     494    (sub imm0 imm0 imm5)
     495    @ok3
     496
     497    (u32-ref imm1 5 seed)
     498    (rlwinm imm2 imm1 15 1 16)
     499    (srwi imm3 imm1 16)
     500    (mulli imm3 imm3 21069)
     501    (add imm2 imm2 imm3)
     502
     503    ;; reduce mod m2
     504    (cmplw cr7 imm2 imm5)
     505    (blt cr7 @ok4)
     506    (sub imm2 imm2 imm5)
     507    @ok4
     508
     509    (add imm2 imm1 imm2)
     510    (cmplw cr7 imm2 imm5)
     511    (blt cr7 @ok5)
     512    (sub imm2 imm2 imm5)
     513    @ok5
     514
     515    (add imm2 imm2 imm0)
     516    (cmplw cr7 imm2 imm5)
     517    (blt cr7 @ok6)
     518    (sub imm2 imm2 imm5)
     519    @ok6
     520
     521    ;; update state
     522    (u32-ref imm0 4 seed)
     523    (u32-set imm0 5 seed)
     524    (u32-ref imm0 3 seed)
     525    (u32-set imm0 4 seed)
     526    (u32-set imm2 3 seed)
     527
     528    ;; construct m1 (1- (expt 2 31))
     529    (lis imm5 #x7fff)
     530    (ori imm5 imm5 #xffff)
     531
     532    ;; combination
     533    (u32-ref imm0 0 seed)
     534    (cmplw cr7 imm0 imm2)
     535    (sub imm0 imm0 imm2)
     536    (bgt cr7 @finish)
     537    (add imm0 imm0 imm5)
     538    @finish
     539    #+ppc32-target
     540    (clrlwi imm0 imm0 3)                ;don't want negative fixnums
     541    (box-fixnum arg_z imm0)
     542    (blr)))   
    443543
    444544; End of ppc-numbers.lisp
  • trunk/source/level-0/X86/X8632/x8632-numbers.lisp

    r13067 r13327  
    219219  (single-value-return))
    220220
     221(defx8632lapfunction %mrg31k3p ((state arg_z))
     222  (let ((seed temp0)
     223        (m1 #x7fffffff)
     224        (m2 #x7fffadb3)
     225        (negative-m1 #x80000001)
     226        (negative-m2 #x8000524d)
     227        (imm1 temp1))
     228    (mark-as-imm temp1)
     229    (svref state 1 seed)
     230    (movl (@ (+ x8632::misc-data-offset (* 4 1)) (% seed)) (% imm0))
     231    (andl ($ #x1ff) (% imm0))
     232    (shll ($ 22) (% imm0))
     233    (movl (@ (+ x8632::misc-data-offset (* 4 1)) (% seed)) (% imm1))
     234    (shrl ($ 9) (% imm1))
     235    (addl (% imm1) (% imm0))
     236
     237    (movl (@ (+ x8632::misc-data-offset (* 4 2)) (% seed)) (% imm1))
     238    (andl ($ #xffffff) (% imm1))
     239    (shll ($ 7) (% imm1))
     240    (addl (% imm1) (% imm0))
     241    (movl (@ (+ x8632::misc-data-offset (* 4 2)) (% seed)) (% imm1))
     242    (shrl ($ 24) (% imm1))
     243
     244    (addl (% imm1) (% imm0))
     245    (leal (@ negative-m1 (% imm0)) (% imm1))
     246    (cmpl ($ m1) (% imm0))
     247    (cmovael (% imm1) (% imm0))
     248
     249    (addl (@ (+ x8632::misc-data-offset (* 4 2)) (% seed)) (% imm0))
     250    (leal (@ negative-m1 (% imm0)) (% imm1))
     251    (cmpl ($ m1) (% imm0))
     252    (cmovael (% imm1) (% imm0))
     253
     254    ;; update state
     255    (movl (@ (+ x8632::misc-data-offset (* 4 1)) (% seed)) (% imm1))
     256    (movl (% imm1) (@ (+ x8632::misc-data-offset (* 4 2)) (% seed)))
     257    (movl (@ (+ x8632::misc-data-offset (* 4 0)) (% seed)) (% imm1))
     258    (movl (% imm1) (@ (+ x8632::misc-data-offset (* 4 1)) (% seed)))
     259    (movl (% imm0) (@ (+ x8632::misc-data-offset (* 4 0)) (% seed)))
     260
     261    ;; second component
     262    (movzwl (@ (+ x8632::misc-data-offset (* 4 3)) (% seed)) (% imm0))
     263    ;(andl ($ #xffff) (% imm0))
     264    (shll ($ 15) (% imm0))
     265    (movl (@ (+ x8632::misc-data-offset (* 4 3)) (% seed)) (% imm1))
     266    (shrl ($ 16) (% imm1))
     267    (imull ($ 21069) (% imm1) (% imm1))
     268
     269    (addl (% imm1) (% imm0))
     270    (leal (@ negative-m2 (% imm0)) (% imm1))
     271    (cmpl ($ m2) (% imm0))
     272    (cmovael (% imm1) (% imm0))
     273    (movl (% imm0) (:rcontext x8632::tcr.unboxed0))     ;stash t1
     274
     275    (movzwl (@ (+ x8632::misc-data-offset (* 4 5)) (% seed)) (% imm0))
     276    ;(andl ($ #xffff) (% imm0))
     277    (shll ($ 15) (% imm0))
     278    (movl (@ (+ x8632::misc-data-offset (* 4 5)) (% seed)) (% imm1))
     279    (shrl ($ 16) (% imm1))
     280    (imull ($ 21069) (% imm1) (% imm1))
     281
     282    (addl (% imm1) (% imm0))
     283    (leal (@ negative-m2 (% imm0)) (% imm1))
     284    (cmpl ($ m2) (% imm0))
     285    (cmovael (% imm1) (% imm0))
     286
     287    (addl (@ (+ x8632::misc-data-offset (* 4 5)) (% seed)) (% imm0))
     288    (leal (@ negative-m2 (% imm0)) (% imm1))
     289    (cmpl ($ m2) (% imm0))
     290    (cmovael (% imm1) (% imm0))
     291
     292    (addl (:rcontext x8632::tcr.unboxed0) (% imm0))     ;add in t1
     293    (leal (@ negative-m2 (% imm0)) (% imm1))
     294    (cmpl ($ m2) (% imm0))
     295    (cmovael (% imm1) (% imm0))
     296
     297    ;; update state
     298    (movl (@ (+ x8632::misc-data-offset (* 4 4)) (% seed)) (% imm1))
     299    (movl (% imm1) (@ (+ x8632::misc-data-offset (* 4 5)) (% seed)))
     300    (movl (@ (+ x8632::misc-data-offset (* 4 3)) (% seed)) (% imm1))
     301    (movl (% imm1) (@ (+ x8632::misc-data-offset (* 4 4)) (% seed)))
     302    (movl (% imm0) (@ (+ x8632::misc-data-offset (* 4 3)) (% seed)))
     303
     304    ;; combination
     305    (movl (@ (+ x8632::misc-data-offset (* 4 0)) (% seed)) (% imm1))
     306    (xchgl (% imm1) (% imm0))           ;for sanity
     307    (rcmpl (% imm0) (% imm1))
     308    (ja @ok)
     309    (subl (% imm1) (% imm0))
     310    (mark-as-node temp1)
     311    (addl ($ m1) (% imm0))
     312    (box-fixnum imm0 arg_z)
     313    (andl ($ #x7fffffff) (% arg_z))
     314    (single-value-return)
     315    @ok
     316    (subl (% imm1) (% imm0))
     317    (mark-as-node temp1)
     318    (box-fixnum imm0 arg_z)
     319    (andl ($ #x7fffffff) (% arg_z))
     320    (single-value-return)))
  • trunk/source/level-0/X86/x86-numbers.lisp

    r13067 r13327  
    190190    (single-value-return)))
    191191
    192 
     192(defx86lapfunction %mrg31k3p ((state arg_z))
     193  (let ((seed temp0)
     194        (m1 #x7fffffff)
     195        (m2 #x7fffadb3)
     196        (negative-m1 #x80000001)
     197        (negative-m2 #x8000524d))
     198    (svref state 1 seed)
     199    (movl (@ (+ x8664::misc-data-offset (* 4 1)) (% seed)) (% imm0.l))
     200    (andl ($ #x1ff) (% imm0.l))
     201    (shll ($ 22) (% imm0.l))
     202    (movl (@ (+ x8664::misc-data-offset (* 4 1)) (% seed)) (% imm1.l))
     203    (shrl ($ 9) (% imm1.l))
     204    (addl (% imm1.l) (% imm0.l))
     205
     206    (movl (@ (+ x8664::misc-data-offset (* 4 2)) (% seed)) (% imm1.l))
     207    (andl ($ #xffffff) (% imm1.l))
     208    (shll ($ 7) (% imm1.l))
     209    (addl (% imm1.l) (% imm0.l))
     210    (movl (@ (+ x8664::misc-data-offset (* 4 2)) (% seed)) (% imm1.l))
     211    (shrl ($ 24) (% imm1.l))
     212
     213    (addl (% imm1.l) (% imm0.l))
     214    (leal (@ negative-m1 (% imm0.l)) (% imm1.l))
     215    (cmpl ($ m1) (% imm0.l))
     216    (cmovael (% imm1.l) (% imm0.l))
     217
     218    (addl (@ (+ x8664::misc-data-offset (* 4 2)) (% seed)) (% imm0.l))
     219    (leal (@ negative-m1 (% imm0.l)) (% imm1.l))
     220    (cmpl ($ m1) (% imm0.l))
     221    (cmovael (% imm1.l) (% imm0.l))
     222
     223    ;; update state
     224    (movl (@ (+ x8664::misc-data-offset (* 4 1)) (% seed)) (% imm1.l))
     225    (movl (% imm1.l) (@ (+ x8664::misc-data-offset (* 4 2)) (% seed)))
     226    (movl (@ (+ x8664::misc-data-offset (* 4 0)) (% seed)) (% imm1.l))
     227    (movl (% imm1.l) (@ (+ x8664::misc-data-offset (* 4 1)) (% seed)))
     228    (movl (% imm0.l) (@ (+ x8664::misc-data-offset (* 4 0)) (% seed)))
     229
     230    ;; second component
     231    (movl (@ (+ x8664::misc-data-offset (* 4 3)) (% seed)) (% imm0.l))
     232    (andl ($ #xffff) (% imm0.l))
     233    (shll ($ 15) (% imm0.l))
     234    (movl (@ (+ x8664::misc-data-offset (* 4 3)) (% seed)) (% imm1.l))
     235    (shrl ($ 16) (% imm1.l))
     236    (imull ($ 21069) (% imm1.l) (% imm1.l))
     237
     238    (addl (% imm1.l) (% imm0.l))
     239    (leal (@ negative-m2 (% imm0.l)) (% imm1.l))
     240    (cmpl ($ m2) (% imm0.l))
     241    (cmovael (% imm1.l) (% imm0.l))
     242
     243    (movl (% imm0.l) (% imm2.l))        ;stash t1
     244
     245    (movl (@ (+ x8664::misc-data-offset (* 4 5)) (% seed)) (% imm0.l))
     246    (andl ($ #xffff) (% imm0.l))
     247    (shll ($ 15) (% imm0.l))
     248    (movl (@ (+ x8664::misc-data-offset (* 4 5)) (% seed)) (% imm1.l))
     249    (shrl ($ 16) (% imm1.l))
     250    (imull ($ 21069) (% imm1.l) (% imm1.l))
     251
     252    (addl (% imm1.l) (% imm0.l))
     253    (leal (@ negative-m2 (% imm0.l)) (% imm1.l))
     254    (cmpl ($ m2) (% imm0.l))
     255    (cmovael (% imm1.l) (% imm0.l))
     256
     257    (addl (@ (+ x8664::misc-data-offset (* 4 5)) (% seed)) (% imm0.l))
     258    (leal (@ negative-m2 (% imm0.l)) (% imm1.l))
     259    (cmpl ($ m2) (% imm0.l))
     260    (cmovael (% imm1.l) (% imm0.l))
     261
     262    (addl (% imm2.l) (% imm0.l))        ;add in t1
     263    (leal (@ negative-m2 (% imm0.l)) (% imm1.l))
     264    (cmpl ($ m2) (% imm0.l))
     265    (cmovael (% imm1.l) (% imm0.l))
     266
     267    ;; update state
     268    (movl (@ (+ x8664::misc-data-offset (* 4 4)) (% seed)) (% imm1.l))
     269    (movl (% imm1.l) (@ (+ x8664::misc-data-offset (* 4 5)) (% seed)))
     270    (movl (@ (+ x8664::misc-data-offset (* 4 3)) (% seed)) (% imm1.l))
     271    (movl (% imm1.l) (@ (+ x8664::misc-data-offset (* 4 4)) (% seed)))
     272    (movl (% imm0.l) (@ (+ x8664::misc-data-offset (* 4 3)) (% seed)))
     273
     274    ;; combination
     275    (movl (@ (+ x8664::misc-data-offset (* 4 0)) (% seed)) (% imm1.l))
     276    (xchgl (% imm1.l) (% imm0.l))               ;for sanity
     277    (rcmpl (% imm0.l) (% imm1.l))
     278    (ja @ok)
     279    (subl (% imm1.l) (% imm0.l))
     280    (addl ($ m1) (% imm0.l))
     281    (box-fixnum imm0 arg_z)
     282    (single-value-return)
     283    @ok
     284    (subl (% imm1.l) (% imm0.l))
     285    (box-fixnum imm0 arg_z)
     286    (single-value-return)))
    193287
    194288;;; End of x86-numbers.lisp
  • trunk/source/level-0/l0-bignum64.lisp

    r13067 r13327  
    22272227    (with-bignum-buffers ((bignum ndigits))
    22282228      (dotimes (i sign-index)
    2229         (setf (bignum-ref bignum i) (%next-random-seed state)))
     2229        (setf (bignum-ref bignum i) (random (expt 2 digit-size) state)))
    22302230      (setf (bignum-ref bignum sign-index)
    22312231            (logand #x7fffffff (the (unsigned-byte 32)
    2232                                  (%next-random-seed state))))
     2232                                 (random (expt 2 (1- digit-size)) state))))
    22332233      (let* ((result (mod bignum number)))
    22342234        (if (eq result bignum)
  • 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)
  • trunk/source/level-1/l1-aprims.lisp

    r13067 r13327  
    5656(def-standard-initial-binding *package*)
    5757(def-standard-initial-binding *gensym-counter* 0)
    58 (def-standard-initial-binding *random-state* (initialize-random-state #xFBF1 9))
     58(def-standard-initial-binding *random-state* (initial-random-state))
    5959(def-standard-initial-binding *whostate* "Reset")
    6060(setq *whostate* "Reset")
  • 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
  • trunk/source/level-1/sysutils.lisp

    r13184 r13327  
    512512        ((and (hash-table-p x) (hash-table-p y))
    513513         (%hash-table-equalp x y))
     514        ((and (random-state-p x) (random-state-p y))
     515         (%random-state-equalp x y))
    514516        (t nil)))
    515517
  • trunk/source/lib/numbers.lisp

    r13067 r13327  
    347347(defparameter a-short-float 1.0s0)
    348348
    349 #+32-bit-target
    350349(defmethod print-object ((rs random-state) stream)
    351   (format stream "#.(~S ~S ~S)"         ;>> #.GAG!!!
    352           'ccl::initialize-random-state
    353           (random.seed-1 rs)
    354           (random.seed-2 rs)))
    355 
    356 #+64-bit-target
    357 (defmethod print-object ((rs random-state) stream)
    358   (let* ((s1 (random.seed-1 rs)))
    359     (format stream "#.(~S ~S ~S)"       ;>> #.GAG!!!
    360             'ccl::initialize-random-state
    361             (ldb (byte 16 16) s1)
    362             (ldb (byte 16 0) s1))))
    363 
    364 
     350  (let* ((s1 (random.mrg31k3p-state rs)))
     351    (format stream "#.(~s~{ ~s~})"       ;>> #.GAG!!!
     352            'ccl::initialize-mrg31k3p-state
     353            (coerce s1 'list))))
    365354
    366355(defun float-radix (float)
  • trunk/source/library/lispequ.lisp

    r13279 r13327  
    824824(def-accessors (random-state) %svref
    825825  ()
    826   random.seed-1
    827   random.seed-2)
     826  random.mrg31k3p-state)
    828827
    829828;;; IEEE-floating-point constants.
Note: See TracChangeset for help on using the changeset viewer.