ergun-test.lisp

(in-package "NAPA-FFT.TESTS")

(defun make-forward-fun (size)
  (compile nil `(lambda (vec)
                  (declare (type complex-sample-array vec)
                           (optimize speed))
                  (let ((twiddle ,(make-twiddle size))
                        (start   0))
                    (declare (type complex-sample-array twiddle)
                             (type (eql 0) start))
                    twiddle
                    ,(gen-dif size))
                  vec)))

(defun check-eqv (a b &optional (name "Test"))
  (multiple-value-bind (ok diff index)
      (m= a b)
    (unless ok
      (error "~A failed with delta ~A (~A)~%" name diff index))))

(defun %forward-test-1 (size repeat function)
  (let ((a (make-vector size))
        (b (make-vector size))
        (sum (make-vector size)))
    (loop repeat repeat do
      (random-vector size a)
      (random-vector size b)
      (m+ a b sum)
      (funcall function a)
      (funcall function b)
      (funcall function sum)
      (m+ a b a)
      (check-eqv a sum "Forward-test-1"))
    t))

(defun %forward-test-2 (size repeat function)
  (let ((r    (make-vector size))
        (f[r] (make-vector size))
        (a    (make-vector size))
        (diff (make-vector size)))
    (setf (aref r 0) (complex 1d0 0d0))
    (fill f[r] (complex 1d0 0d0))
    (loop repeat repeat do
      (random-vector size a)
      (m- r a diff)
      (funcall function a)
      (funcall function diff)
      (m+ a diff a)
      (check-eqv a f[r] "Forward-test-2"))))

(defun rol (vec &optional (dst vec))
  (declare (type complex-sample-array vec dst))
  (let ((last (aref vec (1- (length vec)))))
    (replace dst vec :start1 1)
    (setf (aref dst 0) last)
    dst))

(defvar *bit-reversed* nil)

(defun %forward-test-3 (size outer-repeat inner-repeat function)
  (let ((a (make-vector size))
        (b (make-vector size))
        (diff (make-vector size))
        (y1 (make-vector size))
        (y2 (make-vector size)))
    (loop repeat outer-repeat do
      (random-vector size a)
      (setf (aref a (1- size)) (complex 0d0 0d0))
      (replace y1 a)
      (funcall function y1)
      (loop repeat inner-repeat do
        (random-vector size b)
        (m- a b diff)
        (funcall function b)
        (funcall function diff)
        (m+ b diff b)
        (check-eqv y1 b))
      (rol a)
      (replace y2 a)
      (funcall function y2)
      (loop repeat inner-repeat do
        (random-vector size b)
        ;; FIXME: Ergun subtracts (rol b) here. How does that make sense?
        (m- a b diff)
        (funcall function b)
        (funcall function diff)
        (m+ b diff b)
        (check-eqv y2 b))
      (let ((y1 (if *bit-reversed*
                    (slow-bit-reverse y1)
                    y1))
            (y2 (if *bit-reversed*
                    (slow-bit-reverse y2)
                    y2))
            (root (exp (* -2 pi (complex 0 1d0) (/ size))))
            (mul (complex 1d0 0d0)))
        (declare (type complex-sample root mul)
                 (type complex-sample-array y1 y2))
        (dotimes (i size)
          (setf (aref y1 i) (* mul (aref y1 i))
                mul         (* mul root)))
        (if *bit-reversed*
            (check-eqv (slow-bit-reverse y1)
                       (slow-bit-reverse y2))
            (check-eqv y1 y2))))))

(defun forward-test (size &key (prob 1d-5)
                            (maker 'make-forward-fun)
                            ((:bit-reversed *bit-reversed*) t))
  (let ((repeat (ceiling (log (/ 2d0 prob) 2d0)))
        (fun    (funcall maker size)))
    (assert (plusp repeat))
    (%forward-test-1 size repeat fun)
    (%forward-test-2 size repeat fun)
    (%forward-test-3 size repeat (+ repeat 2) fun)))

(defun run-forward-tests (max-size &optional fancy (*fancy-in-order* t))
  (loop for i upto max-size
        do (forward-test (ash 1 i)
                         :maker (if fancy
                                    (lambda (n)
                                      (get-fft n :in-order *fancy-in-order*))
                                    'make-forward-fun)
                         :bit-reversed (not (and fancy
                                                 *fancy-in-order*)))))