; 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)