;;; FFT - Fast Fourier Transform, translated from "Numerical Recipes in C"
  
(import (rnrs base)
        (rnrs io simple)
        (rnrs arithmetic flonums))

;(define flsin sin)

(define (four1 data)
  (let ((n (vector-length data))
        (pi*2 6.28318530717959)) ; to compute the inverse, negate this value

    ; bit-reversal section

    (let loop1 ((i 0) (j 0))
      (if (< i n)
        (begin
          (if (< i j)
            (begin
              (let ((temp (vector-ref data i)))
                (vector-set! data i (vector-ref data j))
                (vector-set! data j temp))
              (let ((temp (vector-ref data (+ i 1))))
                (vector-set! data (+ i 1) (vector-ref data (+ j 1)))
                (vector-set! data (+ j 1) temp))))
          (let loop2 ((m (div n 2)) (j j))
            (if (and (>= m 2) (>= j m))
              (loop2 (div m 2) (- j m))
              (loop1 (+ i 2) (+ j m)))))))

    ; Danielson-Lanczos section

    (let loop3 ((mmax 2))
      (if (< mmax n)
        (let* ((theta
                (fl/ pi*2 (inexact mmax)))
               (wpr
                (let ((x (flsin (fl* 0.5 theta))))
                  (fl* -2.0 (fl* x x))))
               (wpi
                (flsin theta)))
          (let loop4 ((wr 1.0) (wi 0.0) (m 0))
            (if (< m mmax)
              (begin
                (let loop5 ((i m))
                  (if (< i n)
                    (let* ((j
                            (+ i mmax))
                           (tempr
                            (fl-
                              (fl* wr (vector-ref data j))
                              (fl* wi (vector-ref data (+ j 1)))))
                           (tempi
                            (fl+
                              (fl* wr (vector-ref data (+ j 1)))
                              (fl* wi (vector-ref data j)))))
                      (vector-set! data j
                        (fl- (vector-ref data i) tempr))
                      (vector-set! data (+ j 1)
                        (fl- (vector-ref data (+ i 1)) tempi))
                      (vector-set! data i
                        (fl+ (vector-ref data i) tempr))
                      (vector-set! data (+ i 1)
                        (fl+ (vector-ref data (+ i 1)) tempi))
                      (loop5 (+ j mmax)));***))
                (loop4 (fl+ (fl- (fl* wr wpr) (fl* wi wpi)) wr)
                       (fl+ (fl+ (fl* wi wpr) (fl* wr wpi)) wi)
                       (+ m 2)))))
));******
          (loop3 (* mmax 2)))))))

(define data
  (make-vector 1024 0.0))
 
(define (run data)
  (four1 data)
  (vector-ref data 0))

(define (main)
  (let* ((count (read))
         (input1 (read))
         (input2 (read))
         (output (read))
         (s2 (number->string count))
         (s1 (number->string input1))
         (name "fft"))
    (run-r6rs-benchmark
     (string-append name ":" s1 ":" s2)
     count
     (lambda ()
       (run (hide count (make-vector input1 input2))))
     (lambda (result) (equal? result output)))))