1(define times 1000)
2(define size 1024)
3
4(define (b_fft areal aimag)
5  (let ((ar 0)
6        (ai 0)
7        (i 1)
8        (j 1)
9        (k 0)
10        (m 0)
11        (n 0)
12        (le 1)
13        (le1 0) (le2 0)
14        (ip 0)
15        (nv2 0)
16        (nm1 0)
17        (ur 0.0)
18        (ui 0.0)
19        (wr 0.0)
20        (wi 0.0)
21        (tr 0.0)
22        (ti 0.0))
23    ;; initialize
24    (set! ar areal)
25    (set! ai aimag)
26    (set! n (length ar))
27    (set! n (- n 1))
28    (set! nv2 (quotient n 2))
29    (set! nm1 (- n 1))
30
31    (do ()
32	((>= i n))
33      (set! m (+ m 1))
34      (set! i (+ i i)))
35    (set! i 1)
36
37    (do ()
38	((>= i n))
39      (when (< i j)
40	(set! tr (float-vector-ref ar j))
41	(set! ti (float-vector-ref ai j))
42	(float-vector-set! ar j (float-vector-ref ar i))
43	(float-vector-set! ai j (float-vector-ref ai i))
44	(float-vector-set! ar i tr)
45	(float-vector-set! ai i ti))
46      (set! k nv2)
47      (do ()
48	  ((>= k j))
49	(set! j (- j k))
50	(set! k (quotient k 2)))
51
52      (set! j (+ j k))
53      (set! i (+ i 1)))
54
55    (do ((l 1 (+ l 1)))
56        ((> l m))
57      (set! le1 le)
58      (set! le2 (+ le 1))
59      (set! le (* le 2))
60      (set! ur 1.0)
61      (set! ui 0.)
62      (set! wr (cos (/ pi le1)))
63      (set! wi (sin (/ pi le1)))
64      (do ((j1 1 (+ j1 1)))
65	  ((= j1 le2))
66	(do ((i1 j1 (+ i1 le)))
67	    ((> i1 n))
68	  (set! ip (+ i1 le1))
69	  (set! tr (- (* (float-vector-ref ar ip) ur)
70		      (* (float-vector-ref ai ip) ui)))
71	  (set! ti (+ (* (float-vector-ref ar ip) ui)
72		      (* (float-vector-ref ai ip) ur)))
73	  (float-vector-set! ar ip (- (float-vector-ref ar i1) tr))
74	  (float-vector-set! ai ip (- (float-vector-ref ai i1) ti))
75	  (float-vector-set! ar i1 (+ (float-vector-ref ar i1) tr))
76	  (float-vector-set! ai i1 (+ (float-vector-ref ai i1) ti))))
77      (set! tr (- (* ur wr) (* ui wi)))
78      (set! ti (+ (* ur wi) (* ui wr)))
79      (set! ur tr)
80      (set! ui ti))
81    #t))
82
83(if (not (defined? 'complex))
84  (define complex make-rectangular))
85(if (not (defined? 'when))
86  (define-macro (when test . body) `(if ,test (begin ,@body))))
87
88(define* (cfft data n (dir 1)) ; complex data
89  (unless n (set! n (length data)))
90  (do ((i 0 (+ i 1))
91       (j 0))
92      ((= i n))
93    (if (> j i)
94	(let ((temp (data j)))
95	  (set! (data j) (data i))
96	  (set! (data i) temp)))
97    (do ((m (/ n 2) (/ m 2)))
98        ((or (< m 2)
99             (< j m))
100         (set! j (+ j m)))
101     (set! j (- j m))))
102  (do ((ipow (floor (log n 2)))
103       (prev 1)
104       (lg 0 (+ lg 1))
105       (mmax 2 (* mmax 2))
106       (pow (/ n 2) (/ pow 2))
107       (theta (complex 0.0 (* pi dir)) (* theta 0.5)))
108      ((= lg ipow))
109    (do ((wpc (exp theta))
110         (wc 1.0)
111         (ii 0 (+ ii 1)))
112	((= ii prev)
113	 (set! prev mmax))
114      (do ((jj 0 (+ jj 1))
115           (i ii (+ i mmax))
116           (j (+ ii prev) (+ j mmax)))
117          ((>= jj pow))
118        (let ((tc (* wc (data j))))
119          (set! (data j) (- (data i) tc))
120          (set! (data i) (+ (data i) tc))))
121      (set! wc (* wc wpc))))
122  data)
123
124(when (defined? 'equivalent?)
125  (unless (equivalent? (cfft (vector 0.0 1+i 0.0 0.0)) #(1+1i -1+1i -1-1i 1-1i))
126    (format *stderr* "cfft 1: ~S~%" (cfft (vector 0.0 1+i 0.0 0.0))))
127  (let-temporarily (((*s7* 'equivalent-float-epsilon) 1e-14))
128    (unless (equivalent? (cfft (vector 0 0 1+i 0 0 0 1-i 0)) #(2 -2 -2 2 2 -2 -2 2))
129      (format *stderr* "cfft 2: ~S~%" (cfft (vector 0 0 1+i 0 0 0 1-i 0))))))
130
131(define (fft-bench)
132  (let ((*re* (make-float-vector (+ size 1) 0.10))
133	(*im* (make-float-vector (+ size 1) 0.10)))
134    (do ((ntimes 0 (+ ntimes 1)))
135	((= ntimes times))
136      (b_fft *re* *im*)))
137
138  (let* ((n 256)
139	 (cdata (make-vector n 0.0)))
140    (do ((i 0 (+ i 1)))
141	((= i times))
142      (fill! cdata 0.0)
143      (vector-set! cdata 2 1+i)
144      (vector-set! cdata (- n 1) 1-i)
145      (cfft cdata))))
146
147(fft-bench)
148
149(when (> (*s7* 'profile) 0)
150  (show-profile 200))
151(exit)
152