; hamming codes

; http://www2.rad.com/networks/1994/err_con/hamming.htm

(define (make-matrix rows columns . value)
  (do ((m (make-vector rows)) (i 0 (+ i 1)))
      ((= i rows) m)
    (if (null? value)
        (vector-set! m i (make-vector columns))
        (vector-set! m i (make-vector columns (car value))))))

(define (matrix-rows x) (vector-length x))

(define (matrix-cols x) (vector-length (vector-ref x 0)))

(define (matrix-ref m i j) (vector-ref (vector-ref m i) j))

(define (matrix-set! m i j x) (vector-set! (vector-ref m i) j x))

(define-syntax for
  (syntax-rules ()
    ((for (var first past step) body ...)
      (let ((ge? (if (< first past) >= <=)))
        (do ((var first (+ var step)))
            ((ge? var past))
          body ...)))
    ((for (var first past) body ...)
      (let* ((f first) (p past) (s (if (< first past) 1 -1)))
        (for (var f p s) body ...)))
    ((for (var past) body ...)
      (let* ((p past)) (for (var 0 p) body ...)))))

(define (column v k)
  (if (<= (vector-length (vector-ref v 0)) k)
      (error 'column "out of bounds")
      (let ((z (make-matrix (matrix-rows v) 1)))
        (do ((r 0 (+ r 1))) ((= r (matrix-rows v)) z)
          (matrix-set! z r 0 (matrix-ref v r k))))))

(define (vector-map f v)
  (for (i (vector-length v))
    (vector-set! v i (f (vector-ref v i))))
  v)

(define (subvector v first past)
  (let* ((len (- past first))
         (z (make-vector len)))
    (do ((from first (+ from 1)) (to 0 (+ to 1)))
        ((= from past) z)
      (vector-set! z to (vector-ref v from)))))

(define (mmul2 x y) ; matrix multiplication mod 2
  (let* ((x-rows (matrix-rows x)) (x-cols (matrix-cols x))
         (y-rows (matrix-rows y)) (y-cols (matrix-cols y))
         (z-rows x-rows) (z-cols y-cols))
    (if (not (= x-cols y-rows))
        (error 'mul2 "incompatible matrices")
        (let ((z (make-matrix z-rows z-cols 0)))
          (for (i x-rows)
            (for (j y-cols)
              (for (k x-cols)
                (matrix-set! z i j
                  (modulo
                    (+ (matrix-ref z i j)
                       (* (matrix-ref x i k)
                          (matrix-ref y k j)))
                    2)))))
          z))))

(define g #( #(1 0 0 0 1 1 1)
             #(0 1 0 0 0 1 1)
             #(0 0 1 0 1 0 1)
             #(0 0 0 1 1 1 0)))

(define h #( #(1 0 1 1 1 0 0)
             #(1 1 0 1 0 1 0)
             #(1 1 1 0 0 0 1)))

(define (encode c)
  (vector-ref (mmul2 (vector c) g) 0))

(define (decode r)
  (let* ((s (mmul2 h (vector-map vector r)))
         (len (matrix-rows g))
         (z (vector-map (lambda (v) (vector-ref v 0)) (subvector r 0 len))))
    (let loop ((c 0))
      (cond ((= c len) z)
            ((equal? s (column h c))
              (vector-set! z c (- 1 (vector-ref z c))) ; flip bit
              z)
            (else (loop (+ c 1)))))))

(display (encode #(1 0 0 1))) (newline)
(display (decode #(1 0 0 1 0 0 1))) (newline)
(display (decode #(0 0 0 1 0 0 1))) (newline)
(display (decode #(1 1 0 1 0 0 1))) (newline)
(display (decode #(1 0 1 1 0 0 1))) (newline)
(display (decode #(1 0 0 0 0 0 1))) (newline)
(display (decode #(1 0 0 1 1 0 1))) (newline)
(display (decode #(1 0 0 1 0 1 1))) (newline)
(display (decode #(1 0 0 1 0 0 0))) (newline)

#(1 0 0 1 0 0 1)
#(1 0 0 1)
#(1 0 0 1)
#(1 0 0 1)
#(1 0 0 1)
#(1 0 0 1)
#(1 0 0 1)
#(1 0 0 1)
#(1 0 0 1)