;;; Written by R. Matthew Emerson <rme@acm.org> in August 1999,
;;; and placed in the public domain.

(eval-when (compile)
  (declaim (optimize (speed 3) (safety 0) (space 0) (debug 0))))

(defun test-dct (block)
  (declare (type (simple-array single-float (8 8)) block))
  (let ((b (make-array '(8 8) :element-type 'single-float)))
    (declare (type (simple-array single-float (8 8)) b))
    (dotimes (runs 100000)
      (declare (fixnum runs))
      (dotimes (i 8)			;level shift
	(declare (fixnum i))
	(dotimes (j 8)
	  (declare (fixnum j))
	  (setf (aref b i j) (- (aref block i j) 128.0))))
      (forward-dct! b))))

(defparameter *image*
    (make-array '(8 8)
		:element-type 'single-float
		:initial-contents '((139.0 144.0 149.0 153.0 155.0 155.0 155.0 155.0)
				    (144.0 151.0 153.0 156.0 159.0 156.0 156.0 156.0)
				    (150.0 155.0 160.0 163.0 158.0 156.0 156.0 156.0)
				    (159.0 161.0 162.0 160.0 160.0 159.0 159.0 159.0)
				    (159.0 160.0 161.0 162.0 162.0 155.0 155.0 155.0)
				    (161.0 161.0 161.0 161.0 160.0 157.0 157.0 157.0)
				    (162.0 162.0 161.0 163.0 162.0 157.0 157.0 157.0)
				    (162.0 162.0 161.0 161.0 163.0 158.0 158.0 158.0))))

(defun forward-dct! (b)
  "Perform in-place forward scaled DCT on 8x8 array b"
  (declare (type (simple-array single-float(8 8)) b))
  (dotimes (row 8)
    ;; step one
    (let ((p0 (+ (aref b row 0) (aref b row 7)))
          (p1 (+ (aref b row 1) (aref b row 6)))
          (p2 (+ (aref b row 2) (aref b row 5)))
          (p3 (+ (aref b row 3) (aref b row 4)))
          (p4 (- (aref b row 3) (aref b row 4)))
          (p5 (- (aref b row 2) (aref b row 5)))
          (p6 (- (aref b row 1) (aref b row 6)))
          (p7 (- (aref b row 0) (aref b row 7))))
      ;; even half
      ;; step two
      (let ((q0 (+ p0 p3))
            (q1 (+ p1 p2))
            (q2 (- p1 p2))
            (q3 (- p0 p3)))
	;; step three
        (setf (aref b row 0) (+ q0 q1))
        (setf (aref b row 4) (- q0 q1))
        (let ((z (* (+ q2 q3) 0.7071067811865475)))    ; a1
	  ;; step five
          (setf (aref b row 2) (+ q3 z))
          (setf (aref b row 6) (- q3 z))))
      ;; odd half
      ;; step two
      (let* ((q4 (+ p4 p5))			       ; q4 sign reversed
             (q5 (+ p5 p6))
             (q6 (+ p6 p7))
	     ;; step three
             (z (* (- q6 q4) 0.38268343236508984)))    ; a5, q4 sign reversed
	;; step four
        (let ((s4 (- (* q4 0.5411961001461969) z))     ; a2, q4 sign reversed
              (s5 (* q5 0.7071067811865475))	       ; a3
              (s6 (- (* q6 1.3065629648763766) z)))    ; a4
	  ;; step five
          (let ((t5 (+ p7 s5))
                (t7 (- p7 s5)))
	    ;; step six
            (setf (aref b row 5) (+ t7 s4))
            (setf (aref b row 1) (+ t5 s6))
            (setf (aref b row 7) (- t5 s6))
            (setf (aref b row 3) (- t7 s4)))))))
  (dotimes (col 8)
    (let ((p0 (+ (aref b 0 col) (aref b 7 col)))
          (p1 (+ (aref b 1 col) (aref b 6 col)))
          (p2 (+ (aref b 2 col) (aref b 5 col)))
          (p3 (+ (aref b 3 col) (aref b 4 col)))
          (p4 (- (aref b 3 col) (aref b 4 col)))
          (p5 (- (aref b 2 col) (aref b 5 col)))
          (p6 (- (aref b 1 col) (aref b 6 col)))
          (p7 (- (aref b 0 col) (aref b 7 col))))
      ;; even half
      ;; step two
      (let ((q0 (+ p0 p3))
            (q1 (+ p1 p2))
            (q2 (- p1 p2))
            (q3 (- p0 p3)))
	;; step three
        (setf (aref b 0 col) (+ q0 q1))
        (setf (aref b 4 col) (- q0 q1))
        (let ((z (* (+ q2 q3) 0.7071067811865475)))    ; a1
	  ;; step five
          (setf (aref b 2 col) (+ q3 z))
          (setf (aref b 6 col) (- q3 z))))
      ;; odd half
      ;; step two
      (let* ((q4 (+ p4 p5))			       ; q4 sign reversed
             (q5 (+ p5 p6))
             (q6 (+ p6 p7))
	     ;; step three
             (z (* (- q6 q4) 0.38268343236508984)))    ; a5, q4 sign reversed
	;; step four
        (let ((s4 (- (* q4 0.5411961001461969) z))     ; a2, q4 sign reversed
              (s5 (* q5 0.7071067811865475))	       ; a3
              (s6 (- (* q6 1.3065629648763766) z)))    ; a4
	  ;; step five
          (let ((t5 (+ p7 s5))
                (t7 (- p7 s5)))
	    ;; step six
            (setf (aref b 5 col) (+ t7 s4))
            (setf (aref b 1 col) (+ t5 s6))
            (setf (aref b 7 col) (- t5 s6))
            (setf (aref b 3 col) (- t7 s4))))))))