Fwd: Random generator
Bernd Beuster <[email protected]> Sun, 9 Sep 2007 12:10:51 +0200
| Newsgroups | gmane.lisp.openmcl.bugs |
|---|---|
| Message-ID | <[email protected]> |
The typo in %NEXT-RANDOM-PAIR has been fixed in the current CVS
branch. But it looks like the division by 2^31-1 is not properly
implemented:
(defun %next-random-pair-1 (high low)
"Same function as in l0-numbers.lisp.
(NTH-VALUE 1 (%MULTIPLAY ...) has been replaced by *."
(let* ((n (* 48271 (dpb (ldb (byte 15 0) high)
(byte 16 16)
(ldb (byte 16 0) low )))))
(values (ldb (byte 16 16) n)
(ldb (byte 16 0) n))))
(defun %next-random-pair-2 (high low)
"Same as %NEXT-RANDOM-PAIR-1 but with modulo 2^31-1."
(let* ((n (rem (* 48271 (dpb (ldb (byte 15 0) high)
(byte 16 16)
(ldb (byte 16 0) low)))
(1- (expt 2 31)))))
(values (ldb (byte 16 16) n)
(ldb (byte 16 0) n))))
;;; (run-pairs #'ccl::%next-random-pair) ; wrong in snapshot 1.1-pre.
070722, partially fixed in CVS tree
;;; (run-pairs #'%next-random-pair-1) ; no modulo 2^31-1
;;; (run-pairs #'%next-random-pair-2) ; modulo 2^31-1 (correct
algorithm)
(defun run-pairs (fn &optional (n 5))
"Run the function FN N times and print every result."
(let ((hi 0)
(lo 1))
(loop repeat n do
(multiple-value-bind (high low) (funcall fn hi lo)
(setf hi high lo low)
(format t "~&~5D ~5D" high low)))))
The result of (RUN-PAIRS #'fn n) is:
1) (rem (expt 48271 n) (expt 2 32))
2) (rem (expt 48271 n) (1- (expt 2 31)))
Also, the high-low pairs are awkward (at least to me). Maybe we can
use a 64-bit number for *RANDOM-STATE* for the x86-64 implementation
and stay with the 32-bit pairs for the 32-bit implementation?
DPB/LDB is (maybe) more time consuming than a simple division by 2^31-
n, so the bit juggling with DPB/LDB is not worth the trouble and we
should simply discard it and use division for x86-64 -- in assembler
if you might.
--
Bernd