dft   xtlang

;; real in -> complex out
(bind-func dft
  (lambda (in:double* out:Complexd* len)
    (let ((n:i64 0) (k:i64 0) (tmp 0.0) (cpx:Complexd* null) (N (i64tod len)))
      (memset (cast out i8*) 0 (convert (* len 16)))
      (dotimes (k len)
        (set! cpx (pref-ptr out k))
        (dotimes (n len)
          (set! tmp (* (i64tod (* k n)) (/ TWOPI N)))
          (tset! cpx 0
                 (+ (tref cpx 0)
                    (* (pref in n) (cos tmp))))
          (tset! cpx 1
                 (- (tref cpx 1)
                    (* (pref in n) (sin tmp))))))
      void)))