; random number test

(define (digits n . args)
  (let ((b (if (null? args) 10 (car args))))
    (let loop ((n n) (d '()))
      (if (zero? n) d
          (loop (quotient n b)
                (cons (modulo n b) d))))))

(define seed 20170707)

(define (lcg) (set! seed (modulo (* 16807 seed) 2147483647)) seed)

(define (zero-pad n len)
  (let* ((ds (digits n 2)))
    (append (make-list (- len (length ds)) 0) ds)))

(define (make-bits n)
  (let loop ((n n) (bits (list)))
    (if (zero? n) bits
      (loop (- n 1) (append (zero-pad (lcg) 31) bits)))))

(define (one? n) (= n 1))

(define (count-ones n)
  (let* ((bits (make-bits n))
         (expected (* n 31 1/2))
         (std-dev (sqrt (/ n 4)))
         (lo (- expected (* 2 std-dev)))
         (hi (+ expected (* 2 std-dev))))
    (display "Number of ones: ")
    (display (length (filter one? bits)))
    (newline)
    (display "Expected: ")
    (display (inexact->exact (round lo)))
    (display " to ")
    (display (inexact->exact (round hi)))
    (newline)))

(count-ones 10000)

(define (max-run n)
  (let* ((bits (make-bits n))
         (expected (- (/ (log (* 0.5 31 n)) (log 0.5))))
         (std-dev (/ (log 2)))
         (lo (- expected (* 2 std-dev)))
         (hi (+ expected (* 2 std-dev))))
    (let loop ((bits bits) (current-run 0) (max-run 0))
      (if (null? bits)
          (begin
            (display "Maximum run length: ")
            (display max-run)
            (newline)
            (display "Expected: ")
            (display (inexact->exact (round lo)))
            (display " to ")
            (display (inexact->exact (round hi)))
            (newline))
          (if (zero? (car bits))
              (loop (cdr bits) 0 (max current-run max-run))
              (loop (cdr bits) (+ current-run 1) max-run))))))

(max-run 10000)