; gb-flip

(define two31 #x80000000)

(define a (make-vector 56 -1))

(define fptr #f)

(define (mod-diff x y)
  (modulo (- x y) two31))

; use this if you have logand
; (define (mod-diff x y)
;   (logand (- x y) #x7FFFFFFF))

(define (flip-cycle)
  (do ((ii 1 (+ ii 1)) (jj 32 (+ jj 1))) ((< 55 jj))
    (vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj))))
  (do ((ii 25 (+ ii 1)) (jj 1 (+ jj 1))) ((< 55 ii))
    (vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj))))
  (set! fptr 54)
  (vector-ref a 55))

(define (init-rand seed)
  (let* ((seed (mod-diff seed 0)) (prev seed) (next 1))
    (vector-set! a 55 prev)
    (do ((i 21 (modulo (+ i 21) 55))) ((zero? i))
      (vector-set! a i next)
      (set! next (mod-diff prev next))
      (set! seed (+ (quotient seed 2) (if (odd? seed) #x40000000 0)))
      (set! next (mod-diff next seed))
      (set! prev (vector-ref a i)))
    (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle)))

(define (next-rand)
  (if (negative? (vector-ref a fptr)) (flip-cycle)
    (let ((next (vector-ref a fptr)))
      (set! fptr (- fptr 1)) next)))

(define (unif-rand m)
  (let ((t (- two31 (modulo two31 m))))
    (let loop ((r (next-rand)))
      (if (<= t r) (loop (next-rand)) (modulo r m)))))

(define (test-flip)
  (init-rand -314159)
  (when (not (= (next-rand) 119318998))
    (error 'test-flip "Failure on the first try!"))
  (do ((j 1 (+ j 1))) ((< 133 j)) (next-rand))
  (when (not (= (unif-rand #x55555555) 748103812))
    (error 'test-flip "Failure on the second try!"))
  (display "OK, the gb-flip routines seem to work!") (newline))

(test-flip)

OK, the gb-flip routines seem to work!