;; -*- Lisp -*-
(defprop fft (fft fasl dsk liblsp) autoload)
(defprop ift (fft fasl dsk liblsp) autoload)
(defprop complex-to-polar (fft fasl dsk liblsp) autoload)
(defprop polar-to-complex (fft fasl dsk liblsp) autoload)

(defun *make-array args
       (cond ((equal (errset (arraydims (arg 1.)) nil)
                     (listify (- 1 args)))
	      (cond ((eq (arg 2.) 'fixnum) (fillarray (arg 1.) '(0.)))
		    ((eq (arg 2.) 'flonum) (fillarray (arg 1.) '(0.0)))
		    (t (fillarray (arg 1.) '(nil)))))
	     ((apply (function array) (listify args))
	      (set (arg 1.) (get (arg 1.) 'array))))
       (eval (arg 1.)))

(defun $fft (rary iary)
       (fft-arg-check rary) (fft-arg-check iary)
       (fft (get rary 'array) (get iary 'array))
       (list '(mlist) rary iary))

(defun $ift (rary iary)
       (fft-arg-check rary) (fft-arg-check iary)
       (ift (get rary 'array) (get iary 'array))
       (list '(mlist) rary iary))

(defun $recttopolar (rary iary)
       (fft-arg-check rary) (fft-arg-check iary)
       (complex-to-polar1 (get rary 'array) (get iary 'array))
       (list '(mlist) rary iary))

(defun $polartorect (rary iary)
       (fft-arg-check rary) (fft-arg-check iary)
       (polar-to-complex (get rary 'array) (get iary 'array))
       (list '(mlist) rary iary))

(defun fft-arg-check (ary)
       (cond
	((not (and (get ary 'array) (eq (car (arraydims ary)) 'flonum)))
	 (displa ary)
	 (error
	  '|arg to fft//ift//recttopolar//polartorect must be floating point array|))))

(declare (flonum x y xa ya)
	 (fixnum i n m1 m2))

(defun complex-to-polar1 (rary iary)
       ((lambda (ll m1 m2 n)
		(cond ((or (> (length ll) 2)
			   (not (equal ll (cdr (arraydims iary)))))
		       (error '|bad array dimensions - complex-to-polar|)))
		(cond ((= (length ll) 1)
		       (setq n (car ll) ll nil))
		      (t (setq m1 (car ll) m2 (cadr ll)
			       n (* m1 m2) ll t)
			 (*rearray rary 'flonum n)
			 (*rearray iary 'flonum n)))
		(do ((i (1- n) (1- i)) (x 0.0) (y 0.0) (xa 0.0) (ya 0.0))
		    ((< i 0)
		     (cond (ll (*rearray rary 'flonum m1 m2)
			       (*rearray iary 'flonum m1 m2)))
		     (list rary iary))
		    (setq x (arraycall flonum rary i)
			  y (arraycall flonum iary i))
		    (setq xa (abs x) ya (abs y))
		    (and (> ya xa) (setq xa (prog2 nil ya (setq ya xa))))
		    (setq ya (//$ ya xa))
		    (and (or (> ya 1.0) (< ya 1.e-12)) (setq ya 0.0))
		    ; takes care of underflow in division and in squaring
		    (store (arraycall flonum rary i)
			   (setq xa (*$ xa (sqrt (1+$ (*$ ya ya))))))
		    (store (arraycall flonum iary i)
			   (cond ((= xa 0.0) 0.0)
				 (t (setq ya (atan (abs y) x))
				    (cond ((< y 0.0) (-$ ya))
					  (t ya)))))))
	(cdr (arraydims rary)) 0 0 0))