;;; A uniform [0,1] random number generator; is Pierre L'Ecuyer's generator from his paper
;;; "Good parameters and implementations for combined multiple recursive random number generators"
;;; available at his web site http://www.iro.umontreal.ca/~lecuyer

(define seed-set! #f)
(define seed-ref #f)
(define random #f)

(let ((norm 2.328306549295728e-10)
      (m1 4294967087.0)
      (m2 4294944443.0)
      (a12 1403580.0)
      (a13n 810728.0)
      (a21 527612.0)
      (a23n 1370589.0)
      (seed (f64vector 1.0 0.0 0.0 1.0 0.0 0.0)))  ; cannot be a constant, since we mutate it.
  
  (set! random (if (##fixnum? (- (expt 2 32) 1)) ; at least 32 bit fixnums
		   (lambda ()    ; twice as fast as the next version on alpha 21264
		     (declare (flonum))
		     (let ((seed seed))  ; make it local
		       (let ((p1 (- (* a12 (f64vector-ref seed 1)) (* a13n (f64vector-ref seed 0))))
			     (p2 (- (* a21 (f64vector-ref seed 5)) (* a23n (f64vector-ref seed 3)))))
			 (let ((k1 (##flonum.<-fixnum (##flonum.->fixnum (/ p1 m1))))
			       (k2 (##flonum.<-fixnum (##flonum.->fixnum (/ p2 m2))))
			       (ignore1 (f64vector-set! seed 0 (f64vector-ref seed 1)))
			       (ignore3 (f64vector-set! seed 3 (f64vector-ref seed 4))))
			   (let ((p1 (- p1 (* k1 m1)))
				 (p2 (- p2 (* k2 m2)))
				 (ignore2 (f64vector-set! seed 1 (f64vector-ref seed 2)))
				 (ignore4 (f64vector-set! seed 4 (f64vector-ref seed 5))))
			     (let ((p1 (if (< p1 0.0) (+ p1 m1) p1))
				   (p2 (if (< p2 0.0) (+ p2 m2) p2)))
			       (f64vector-set! seed 2 p1)
			       (f64vector-set! seed 5 p2)
			       (if (<= p1 p2)
				   (* norm (+ (- p1 p2) m1))
				   (* norm (- p1 p2)))))))))
		   (lambda ()    ; uses no conversions between flonums and fixnums.
		     (declare (flonum))
		     (let ((seed seed))  ; make it local
		       (let ((p1 (- (* a12 (f64vector-ref seed 1)) (* a13n (f64vector-ref seed 0))))
			     (p2 (- (* a21 (f64vector-ref seed 5)) (* a23n (f64vector-ref seed 3)))))
			 (let ((k1 (truncate (/ p1 m1)))
			       (k2 (truncate (/ p2 m2)))
			       (ignore1 (f64vector-set! seed 0 (f64vector-ref seed 1)))
			       (ignore3 (f64vector-set! seed 3 (f64vector-ref seed 4))))
			   (let ((p1 (- p1 (* k1 m1)))
				 (p2 (- p2 (* k2 m2)))
				 (ignore2 (f64vector-set! seed 1 (f64vector-ref seed 2)))
				 (ignore4 (f64vector-set! seed 4 (f64vector-ref seed 5))))
			     (let ((p1 (if (< p1 0.0) (+ p1 m1) p1))
				   (p2 (if (< p2 0.0) (+ p2 m2) p2)))
			       (f64vector-set! seed 2 p1)
			       (f64vector-set! seed 5 p2)
			       (if (<= p1 p2)
				   (* norm (+ (- p1 p2) m1))
				   (* norm (- p1 p2)))))))))))
	
  (set! seed-ref (lambda () (f64vector->list seed)))
  
  (set! seed-set! (lambda l
		    (set! seed (list->f64vector l)))))

(define (random32)
  (inexact->exact (truncate (* 4294967296. (random)))))

(define (write-diehard filename s bytes-per-call calls)
  (let ((port (open-output-file filename))
        (n (expt 256 bytes-per-call)))
    (do ((i 0 (+ i 1))) ((= i calls) (close-output-port port))
      (let ((x (s)))
        (do ((k 0 (+ k 1))) ((= k bytes-per-call))
          (write-char (integer->char (modulo x 256)) port)
          (set! x (quotient x 256)))))))