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