;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; File:         fft.sc
;;; Description:  FFT benchmark from the Gabriel tests.
;;; Author:       Harry Barrow
;;; Created:      8-Apr-85
;;; Modified:     6-May-85 09:29:22 (Bob Shaw)
;;;               11-Aug-87 (Will Clinger)
;;;               16-Nov-94 (Qobi)
;;;               31-Mar-98 (Qobi)
;;;               26-Mar-00 (flw)
;;; Language:     Scheme
;;; Status:       Public Domain
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define pi (atan2 0 -1))

(define count 0)			;(felix)

;;; FFT -- This is an FFT benchmark written by Harry Barrow.
;;; It tests a variety of floating point operations,
;;; including array references.

(define *re* (make-f64vector 1025 0.0))

(define *im* (make-f64vector 1025 0.0))

(define (fft areal aimag)
 (let ((ar areal)			;Qobi
       (ai aimag)			;Qobi
       (i 0)
       (j 0)
       (k 0)
       (m 0)
       (n 0)
       (le 0)
       (le1 0)
       (ip 0)
       (nv2 0)
       (nm1 0)
       (ur 0.0)				;Qobi
       (ui 0.0)				;Qobi
       (wr 0.0)				;Qobi
       (wi 0.0)				;Qobi
       (tr 0.0)				;Qobi
       (ti 0.0))			;Qobi
  ;; initialize
  (set! ar areal)
  (set! ai aimag)
  (set! n 1025)
  (set! n (- n 1))
  (set! nv2 (quotient n 2))		;(felix) was "quotient"
  (set! nm1 (- n 1))
  (set! m 0)				;compute m = log(n)
  (set! i 1)
  (let loop ()
   (if (< i n)
       (begin (set! m (+ m 1))
	      (set! i (+ i i))
	      (loop))))
  (cond ((not (= n (let ((px 0))
		     (let loop ((i m) (p 1)) ;Qobi
		       (if (zero? i) 
			   (set! px p)	;(felix) omit loop result
			   (loop (- i 1) (* 2 p))))
		     px) ) )
	 (display "array size not a power of two.")
	 (newline)))
  ;; interchange elements in bit-reversed order
  (set! j 1)
  (set! i 1)
  (let l3 ()
   (cond ((< i j)
	  (set! tr (f64vector-ref ar j))
	  (set! ti (f64vector-ref ai j))
	  (f64vector-set! ar j (f64vector-ref ar i))
	  (f64vector-set! ai j (f64vector-ref ai i))
	  (f64vector-set! ar i tr)
	  (f64vector-set! ai i ti)))
   (set! k nv2)
   (let l6 ()
    (cond ((< k j)
	   (set! j (- j k))
	   (set! k (quotient k 2))	;Qobi: was / but this violates R4RS
	   (l6))))
   (set! j (+ j k))
   (set! i (+ i 1))
   (cond ((< i n) (l3))))
  ;; loop thru stages (syntax converted from old MACLISP style \bs)
  (do ((l 1 (+ l 1))) ((> l m))
   (let loop ((i l) (p 1))	;Qobi
     (if (zero? i) 
	 (set! le p)			;(felix) omit loop result
	 (loop (- i 1) (* 2 p))))
   (set! le1 (quotient le 2))
   (set! ur 1.0)
   (set! ui 0.0)
   (set! wr (cos (/ pi le1)))
   (set! wi (sin (/ pi le1)))
   ;; loop thru butterflies
   (do ((j 1 (+ j 1))) ((> j le1))
    ;; do a butterfly
    (do ((i j (+ i le))) ((> i n))
      (set! count (+ count 1))		;(felix)
     (set! ip (+ i le1))
     (set! tr (- (* (f64vector-ref ar ip) ur) (* (f64vector-ref ai ip) ui)))
     (set! ti (+ (* (f64vector-ref ar ip) ui) (* (f64vector-ref ai ip) ur)))
     (f64vector-set! ar ip (- (f64vector-ref ar i) tr))
     (f64vector-set! ai ip (- (f64vector-ref ai i) ti))
     (f64vector-set! ar i (+ (f64vector-ref ar i) tr))
     (f64vector-set! ai i (+ (f64vector-ref ai i) ti))))
   (set! tr (- (* ur wr) (* ui wi)))
   (set! ti (+ (* ur wi) (* ui wr)))
   (set! ur tr)
   (set! ui ti))
  #t))

;;; the timer which does 10 calls on fft

(define (fft-bench)
 (do ((ntimes 0 (+ ntimes 1))) ((= ntimes 10))
  (fft *re* *im*)))

(fft-bench)
(display count)				;(felix)
(newline)