; Checking the average complexity of our filters

(include "lib/myenv.scm")

(define less-counter 0)

(define (less? x y)
  (set! less-counter (+ 1 less-counter))
  (< x y))


(define (iota m n)
  (if (> m n) '()
    (cons m (iota (+ 1 m) n))))

(define (all-permutations lst)
  (cond
    ((null? lst) '(()))
    ((null? (cdr lst)) (list lst))
    (else
      (let ((head (car lst))
	    (res-perm (all-permutations (cdr lst))))
	(concat (map (lambda (l) (insert-all head l)) res-perm))))))

; [[a]] -> [a]
(define (concat llst)
  (let loop ((llst (reverse llst)) (result '()))
    (if (null? llst) result
      (loop (cdr llst) (append (car llst) result)))))

; a -> [a] -> [[a]]
; the result is a list of lists with x inserted in every position
(define (insert-all x lst)
  (let loop ((before (reverse lst)) (suffix '()) (result '()))
    (if (null? before)
      (cons (cons x suffix) result)
      (loop (cdr before) (cons (car before) suffix)
	(cons (append (reverse before) (cons x suffix)) result)))))

; Compute the median in the straightforward, unoptimized way
; as the floor((n-1)/2)-th min of the n-elem list
(define (unopt-median . args)
  (let* ((vec (apply vector args))
	 (n   (vector-length vec)))
    ; Sort vec from i (inclusive) through j (exclusive)
    ; Use the quicksort
    (define (sort i j)
      (let ((n (- j i)))
	(cond
	  ((= n 2)
	   (if (> (vector-ref vec i) (vector-ref vec (++ i)))
	     (swap! i (++ i))))
	  ((> n 2)
	   (let* ((pivot-pos (+ i (quotient n 2)))
		  (pivot-pos
		    (if (> (vector-ref vec i) (vector-ref vec pivot-pos))
		      i pivot-pos))
		  (pivot-pos
		    (if (< (vector-ref vec (-- j)) (vector-ref vec pivot-pos))
		      (-- j) pivot-pos))
		  (pivot (vector-ref vec pivot-pos)))
	     ;(cerr "pivot: " pivot ", pivot-pos " pivot-pos nl)
	     ; invariant: p <= pivot-pos <= q
	     (let loop ((p i) (q (-- j)) (pivot-pos pivot-pos))
	       (cond
		 ((= p q)
		   (sort i p)
		   (sort (++ p) j))
		 ((= p pivot-pos)
		  (if (>= (vector-ref vec q) pivot)
		    (loop p (-- q) pivot-pos)
		    (begin (swap! pivot-pos q) ; pivot-pos < q
		      (loop p q q))))
		 ((< (vector-ref vec p) pivot) (loop (++ p) q pivot-pos))
		 (else
		   (assert (< p q) (< p pivot-pos) (<= pivot-pos q))
		   (swap! p pivot-pos) (loop p q p))))))
	  )))
    (define (swap! i j)
      (let ((val (vector-ref vec i)))
	(vector-set! vec i (vector-ref vec j))
	(vector-set! vec j val)))

    (sort 0 n)
    (vector-ref vec (quotient (-- n) 2))))

		 
; Check the average number of comparisons to compute the median
(define (av-path n median-proc)
  (let ((test-set (all-permutations (iota 1 n))))
    (set! less-counter 0)
    (for-each (lambda (lst) (apply median-proc lst)) test-set)
    (display "total comparisons ") (display less-counter) (newline)
    (display "average: ") (display (/ less-counter (length test-set)))
    (newline) (newline)))

; Knuth (volume 3) says that the best average complexity
; of median3 is 2+2/3
(define (median3 x1 x2 x3)
  (if (less? x2 x3)
    (if (less? x2 x1) (if (less? x3 x1) x3 x1) x2)
    (if (less? x3 x1) (if (less? x2 x1) x2 x1) x3))
)

(define-macro (exhaustive-loop n body)
  (let loop ((i 1) (counters '()))
    (if (> i n) `(,body ,(cons 'list counters))
      (let ((counter (gensym)))
	`(do ((,counter 1 (++ ,counter))) ((> ,counter ,n))
	   ,(loop (+ i 1) (cons counter counters)))))))
(let ()
  (cerr nl "Exhaustive check of median3" nl)
  (exhaustive-loop 3
    (lambda (l)
      (let ((found (apply median3 l))
	    (need  (apply unopt-median l)))
	(assert (= found need)))))
  (cerr nl "done" nl)
)
(av-path 3 median3)

; This gives 96/4! comparisons, or 4
(define (median4 x1 x2 x3 x4)
(if (less? x3 x4)
  (if (less? x2 x1)
    (if (less? x3 x2) (if (less? x4 x2) x4 x2) (if (less? x1 x3) x1 x3))
    (if (less? x3 x1) (if (less? x4 x1) x4 x1) (if (less? x2 x3) x2 x3)))
  (if (less? x2 x1)
    (if (less? x4 x2) (if (less? x3 x2) x3 x2) (if (less? x1 x4) x1 x4))
    (if (less? x4 x1) (if (less? x3 x1) x3 x1) (if (less? x2 x4) x2 x4))))
)

(let ()
  (cerr nl "Exhaustive check of median4" nl)
  (exhaustive-loop 4
    (lambda (l)
      (let ((found (apply median4 l))
	    (need  (apply unopt-median l)))
	(assert (= found need)))))
  (cerr nl "done" nl)
)
(assert (= 1 (median4 1 3 1 1)))
(av-path 4 median4)

; This gives 716/5! =179/30.
; Knuth however shows that the absolute min is 5+13/15 = 176/30
(define (median5 x1 x2 x3 x4 x5)
  (if (less? x3 x4)
    (if (less? x5 x3)
      (if (less? x2 x1)
	(if (less? x3 x1)
	  (if (less? x3 x2) (if (less? x4 x2) x4 x2) x3)
	  (if (less? x5 x1) x1 x5))
	(if (less? x3 x2)
	  (if (less? x3 x1) (if (less? x4 x1) x4 x1) x3)
	  (if (less? x5 x2) x2 x5)))
      (if (less? x5 x4)
	(if (less? x2 x1)
	  (if (less? x5 x1)
	    (if (less? x5 x2) (if (less? x4 x2) x4 x2) x5)
	    (if (less? x3 x1) x1 x3))
	  (if (less? x5 x2)
	    (if (less? x5 x1) (if (less? x4 x1) x4 x1) x5)
	    (if (less? x3 x2) x2 x3)))
	(if (less? x2 x1)
	  (if (less? x4 x1)
	    (if (less? x4 x2) (if (less? x5 x2) x5 x2) x4)
	    (if (less? x3 x1) x1 x3))
	  (if (less? x4 x2)
	    (if (less? x4 x1) (if (less? x5 x1) x5 x1) x4)
	    (if (less? x3 x2) x2 x3)))))
    (if (less? x5 x4)
      (if (less? x2 x1)
	(if (less? x4 x1)
	  (if (less? x4 x2) (if (less? x3 x2) x3 x2) x4)
	  (if (less? x5 x1) x1 x5))
	(if (less? x4 x2)
	  (if (less? x4 x1) (if (less? x3 x1) x3 x1) x4)
	  (if (less? x5 x2) x2 x5)))
      (if (less? x5 x3)
	(if (less? x2 x1)
	  (if (less? x5 x1)
	    (if (less? x5 x2) (if (less? x3 x2) x3 x2) x5)
	    (if (less? x4 x1) x1 x4))
	  (if (less? x5 x2)
	    (if (less? x5 x1) (if (less? x3 x1) x3 x1) x5)
	    (if (less? x4 x2) x2 x4)))
	(if (less? x2 x1)
	  (if (less? x3 x1)
	    (if (less? x3 x2) (if (less? x5 x2) x5 x2) x3)
	    (if (less? x4 x1) x1 x4))
	  (if (less? x3 x2)
	    (if (less? x3 x1) (if (less? x5 x1) x5 x1) x3)
	    (if (less? x4 x2) x2 x4))))))
)
  

(let ()
  (cerr nl "Exhaustive check of median5" nl)
  (exhaustive-loop 5
    (lambda (l)
      (let ((found (apply median5 l))
	    (need  (apply unopt-median l)))
	(assert (= found need)))))
  (cerr nl "done" nl)
)
(assert (= 3 (median5 1 2 3 4 5)))
(assert (= 1 (median5 1 3 3 1 1)))
(av-path 5 median5)

; This gives 5568/6! = 116/15
(define (median6 x1 x2 x3 x4 x5 x6)
 (if (less? x4 x5)
  (if (less? x6 x4)
    (if (less? x3 x2)
      (if (less? x1 x3)
        (if (less? x4 x3)
          (if (less? x4 x1) (if (less? x5 x1) x5 x1) x4)
          (if (less? x6 x3) x3 (if (less? x2 x6) x2 x6)))
        (if (less? x1 x2)
          (if (less? x4 x1)
            (if (less? x4 x3) (if (less? x5 x3) x5 x3) x4)
            (if (less? x6 x1) x1 (if (less? x2 x6) x2 x6)))
          (if (less? x4 x2)
            (if (less? x4 x3) (if (less? x5 x3) x5 x3) x4)
            (if (less? x6 x2) x2 (if (less? x1 x6) x1 x6)))))
      (if (less? x1 x2)
        (if (less? x4 x2)
          (if (less? x4 x1) (if (less? x5 x1) x5 x1) x4)
          (if (less? x6 x2) x2 (if (less? x3 x6) x3 x6)))
        (if (less? x1 x3)
          (if (less? x4 x1)
            (if (less? x4 x2) (if (less? x5 x2) x5 x2) x4)
            (if (less? x6 x1) x1 (if (less? x3 x6) x3 x6)))
          (if (less? x4 x3)
            (if (less? x4 x2) (if (less? x5 x2) x5 x2) x4)
            (if (less? x6 x3) x3 (if (less? x1 x6) x1 x6))))))
    (if (less? x6 x5)
      (if (less? x3 x2)
        (if (less? x1 x3)
          (if (less? x6 x3)
            (if (less? x6 x1) (if (less? x5 x1) x5 x1) x6)
            (if (less? x4 x3) x3 (if (less? x2 x4) x2 x4)))
          (if (less? x1 x2)
            (if (less? x6 x1)
              (if (less? x6 x3) (if (less? x5 x3) x5 x3) x6)
              (if (less? x4 x1) x1 (if (less? x2 x4) x2 x4)))
            (if (less? x6 x2)
              (if (less? x6 x3) (if (less? x5 x3) x5 x3) x6)
              (if (less? x4 x2) x2 (if (less? x1 x4) x1 x4)))))
        (if (less? x1 x2)
          (if (less? x6 x2)
            (if (less? x6 x1) (if (less? x5 x1) x5 x1) x6)
            (if (less? x4 x2) x2 (if (less? x3 x4) x3 x4)))
          (if (less? x1 x3)
            (if (less? x6 x1)
              (if (less? x6 x2) (if (less? x5 x2) x5 x2) x6)
              (if (less? x4 x1) x1 (if (less? x3 x4) x3 x4)))
            (if (less? x6 x3)
              (if (less? x6 x2) (if (less? x5 x2) x5 x2) x6)
              (if (less? x4 x3) x3 (if (less? x1 x4) x1 x4))))))
      (if (less? x3 x2)
        (if (less? x1 x3)
          (if (less? x5 x3)
            (if (less? x5 x1) (if (less? x6 x1) x6 x1) x5)
            (if (less? x4 x3) x3 (if (less? x2 x4) x2 x4)))
          (if (less? x1 x2)
            (if (less? x5 x1)
              (if (less? x5 x3) (if (less? x6 x3) x6 x3) x5)
              (if (less? x4 x1) x1 (if (less? x2 x4) x2 x4)))
            (if (less? x5 x2)
              (if (less? x5 x3) (if (less? x6 x3) x6 x3) x5)
              (if (less? x4 x2) x2 (if (less? x1 x4) x1 x4)))))
        (if (less? x1 x2)
          (if (less? x5 x2)
            (if (less? x5 x1) (if (less? x6 x1) x6 x1) x5)
            (if (less? x4 x2) x2 (if (less? x3 x4) x3 x4)))
          (if (less? x1 x3)
            (if (less? x5 x1)
              (if (less? x5 x2) (if (less? x6 x2) x6 x2) x5)
              (if (less? x4 x1) x1 (if (less? x3 x4) x3 x4)))
            (if (less? x5 x3)
              (if (less? x5 x2) (if (less? x6 x2) x6 x2) x5)
              (if (less? x4 x3) x3 (if (less? x1 x4) x1 x4))))))))
  (if (less? x6 x5)
    (if (less? x3 x2)
      (if (less? x1 x3)
        (if (less? x5 x3)
          (if (less? x5 x1) (if (less? x4 x1) x4 x1) x5)
          (if (less? x6 x3) x3 (if (less? x2 x6) x2 x6)))
        (if (less? x1 x2)
          (if (less? x5 x1)
            (if (less? x5 x3) (if (less? x4 x3) x4 x3) x5)
            (if (less? x6 x1) x1 (if (less? x2 x6) x2 x6)))
          (if (less? x5 x2)
            (if (less? x5 x3) (if (less? x4 x3) x4 x3) x5)
            (if (less? x6 x2) x2 (if (less? x1 x6) x1 x6)))))
      (if (less? x1 x2)
        (if (less? x5 x2)
          (if (less? x5 x1) (if (less? x4 x1) x4 x1) x5)
          (if (less? x6 x2) x2 (if (less? x3 x6) x3 x6)))
        (if (less? x1 x3)
          (if (less? x5 x1)
            (if (less? x5 x2) (if (less? x4 x2) x4 x2) x5)
            (if (less? x6 x1) x1 (if (less? x3 x6) x3 x6)))
          (if (less? x5 x3)
            (if (less? x5 x2) (if (less? x4 x2) x4 x2) x5)
            (if (less? x6 x3) x3 (if (less? x1 x6) x1 x6))))))
    (if (less? x6 x4)
      (if (less? x3 x2)
        (if (less? x1 x3)
          (if (less? x6 x3)
            (if (less? x6 x1) (if (less? x4 x1) x4 x1) x6)
            (if (less? x5 x3) x3 (if (less? x2 x5) x2 x5)))
          (if (less? x1 x2)
            (if (less? x6 x1)
              (if (less? x6 x3) (if (less? x4 x3) x4 x3) x6)
              (if (less? x5 x1) x1 (if (less? x2 x5) x2 x5)))
            (if (less? x6 x2)
              (if (less? x6 x3) (if (less? x4 x3) x4 x3) x6)
              (if (less? x5 x2) x2 (if (less? x1 x5) x1 x5)))))
        (if (less? x1 x2)
          (if (less? x6 x2)
            (if (less? x6 x1) (if (less? x4 x1) x4 x1) x6)
            (if (less? x5 x2) x2 (if (less? x3 x5) x3 x5)))
          (if (less? x1 x3)
            (if (less? x6 x1)
              (if (less? x6 x2) (if (less? x4 x2) x4 x2) x6)
              (if (less? x5 x1) x1 (if (less? x3 x5) x3 x5)))
            (if (less? x6 x3)
              (if (less? x6 x2) (if (less? x4 x2) x4 x2) x6)
              (if (less? x5 x3) x3 (if (less? x1 x5) x1 x5))))))
      (if (less? x3 x2)
        (if (less? x1 x3)
          (if (less? x4 x3)
            (if (less? x4 x1) (if (less? x6 x1) x6 x1) x4)
            (if (less? x5 x3) x3 (if (less? x2 x5) x2 x5)))
          (if (less? x1 x2)
            (if (less? x4 x1)
              (if (less? x4 x3) (if (less? x6 x3) x6 x3) x4)
              (if (less? x5 x1) x1 (if (less? x2 x5) x2 x5)))
            (if (less? x4 x2)
              (if (less? x4 x3) (if (less? x6 x3) x6 x3) x4)
              (if (less? x5 x2) x2 (if (less? x1 x5) x1 x5)))))
        (if (less? x1 x2)
          (if (less? x4 x2)
            (if (less? x4 x1) (if (less? x6 x1) x6 x1) x4)
            (if (less? x5 x2) x2 (if (less? x3 x5) x3 x5)))
          (if (less? x1 x3)
            (if (less? x4 x1)
              (if (less? x4 x2) (if (less? x6 x2) x6 x2) x4)
              (if (less? x5 x1) x1 (if (less? x3 x5) x3 x5)))
            (if (less? x4 x3)
              (if (less? x4 x2) (if (less? x6 x2) x6 x2) x4)
              (if (less? x5 x3) x3 (if (less? x1 x5) x1 x5)))))))))
)

(let ()
  (cerr nl "Exhaustive check of median6" nl)
  (exhaustive-loop 6
    (lambda (l)
      (let ((found (apply median6 l))
	    (need  (apply unopt-median l)))
	(assert (= found need)))))
  (cerr nl "done" nl)
)
(assert (= 1 (median6 5 1 3 1 3 1)))
(av-path 6 median6)