(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))

(defun ptimesfft(r s)
  ;;compute the size Z of the answer. Z is a power of 2.
  ;; compute complex array r, increased to size S
  ;; compute complex array r, increased to size S
  ;; compute two FFTs
  ;; multiply pointwise
  ;; compute inverse FFT  * 1/n
  ;; convert back to array and return answer
  (let* ((lr (length r))
	 (ls (length s))
	 (lans (+ lr ls -1))
	 (z (ash 1 (ceiling (log lans 2)))) ; round up to power of 2
	 (rfft (four1 (v2dfa r z) z))
	 (sfft (four1 (v2dfa s z) z))
	 (prod (prodarray rfft sfft z ))) 
  ;;  (setf *r2* rfft)
    
    
    (dfa2v(four1 prod z :isign -1) lans)))

(defun v2dfa(a &optional (m (length a))	)
  ;;coerce a vector of numbers of length m to an array of length
  ;; 2m, since here complex numbers are stored in 2 adjacent
  ;; locations.
  (let* ((k (length a))
	 (ans (make-array (* 2 m) :element-type 'double-float))
	 (zz  0.0d0))
    (dotimes (i k)
      (setf (aref ans (* 2 i)) (coerce (aref a i) 'double-float)) ;; here we convert.
      (setf (aref ans (1+ (* 2 i))) zz))
    (loop for i from (* 2 k) to (1-(* 2 m)) do (setf  (aref ans i) zz))
    ans))

(defun dfa2v(a &optional (m (ash (length a)-1)))
  ;; Coerce real parts back to integers, more or less.  a is an an
  ;; array of even length.  If you know that there are trailing zeros,
  ;; set the actual length with the optional second parameter m.
  (let* ((k (ash (length a) -1))
	 (ans (make-array m)))
    (declare (type (simple-array double-float (*)) a))
    (dotimes (i m ans)
      (setf (aref ans i)(round  (aref a (* 2 i)) k)))))

;; a simple fft, not optimized much

(defun four1 (data nn &key (isign 1))
 (declare (type (simple-array double-float (*)) data)) 
 (declare (fixnum nn isign))

 (prog ((wr 0d0) (wi 0d0) (wpr 0d0) (wpi 0d0) (wtemp 0d0) 
        (theta 0d0) (tempr 0d0) (tempi 0d0) (j 0) (n 0) (m 0) 
        (mmax 0) (istep 0))
  (declare (double-float wr wi wpr wpi wtemp theta tempr tempi)) 
  (declare (fixnum j n m mmax istep))


  (setf n (* 2 nn)) 
  (setf j 1) 
  (do ((i 1 (+ i 2)))
      ((> i n) t)
      (declare (type fixnum i))
    (when (> j i) 
     (setf tempr (aref data (1- j)))
     (setf tempi (aref data j)) 
     (setf (aref data (1- j)) (aref data (1- i)))
     (setf (aref data j) (aref data i)) 
     (setf (aref data (1- i)) tempr)
     (setf (aref data i) tempi))
    (setf m (floor (/ n 2)))
    label1
    (when (and (>= m 2) (> j m))
     (setf j (- j m)) (setf m (floor (/ m 2)))
     (go label1))
    (setf j (+ j m))) 

  (setf mmax 2) 
  label2 
  (when (> n mmax)
   (setf istep (* 2 mmax))
   (setf theta (/ 6.28318530717959d0 (* isign mmax)))
   (setf wpr (* -2.0d0 (expt (sin (* 0.5d0 theta)) 2)))
   (setf wpi (sin theta)) (setf wr 1.0d0) (setf wi 0.0d0)
   (do ((m 1 (+ m 2)))
       ((> m mmax) t)
       (declare (type fixnum m))
     (do ((i m (+ i istep)))
         ((> i n) t)
         (declare (type fixnum i))
       (setf j (+ i mmax))
       (setf tempr (+ (* wr (aref data (1- j)))
                      (* (* -1d0 wi) (aref data j))))
       (setf tempi (+ (* wr (aref data j))
                      (* wi (aref data (1- j)))))
       (setf (aref data (1- j)) (+ (aref data (1- i)) (- tempr)))
       (setf (aref data j) (+ (aref data i) (* -1d0 tempi)))
       (setf (aref data (1- i)) (+ (aref data (1- i)) tempr))
       (setf (aref data i) (+ (aref data i) tempi)))
     (setf wtemp wr)
     (setf wr (+ (+ (* wr wpr) (* (- wi) wpi)) wr))
     (setf wi (+ (+ (* wi wpr) (* wtemp wpi)) wi)))
   (setf mmax istep)
   (go label2)) 
   (return data)))


(defun prodarray(r s len)
  ;; r and s are the same length arrays
  ;; compute, for i=0, 2, ..., len-2
  ;; ans[i]:=  r[i]*s[i]-r[i+1]*s[i+1] ;; real part
  ;; ans[i+1]:=r[i]*s[i+1]+s[i]*r[i+1] ;; imag part
  (let ((ans (make-array (* 2 len) :element-type 'double-float)))
   (declare (type (simple-array double-float (*)) r s))
  (dotimes (i len ans)
    (let* ((ind (* 2 i))
	   (ind1 (1+ ind))
	   (a (aref r ind))
	   (b (aref r ind1))
	   (c (aref s ind))
	   (d (aref s ind1)))
      (declare (fixnum ind ind1) (double-float a b c d))
      (setf (aref ans ind) (the double-float (- (* a c)(* b d))))
      (setf (aref ans ind1) (+ (* a d)(* b c)))))))