;;; -*- Mode: LISP; Syntax: Common-Lisp; Package: qd; Base: 10 -*-
;;;;;;;;;;;;;;;;;;;FFT!!!;;;;;;;;;;;;;;;;;

;; this is part of qd.lisp now.

(in-package :qd)
(eval-when (compile)
  (load "ga.fasl")
    (load "qd.fasl"))

(eval-when (compile load eval)
; Fourier Transform Spectral Methods
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;Routines translated with permission by Kevin A. Broughan from ;;;;;;;;;;;
;;Numerical Recipies in Fortran Copyright (c) Numerical Recipies 1986, 1989;;;;
;;;;;;;;;;;;;;;Modified by Ken Olum for Common Lisp, April 1996;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; notes 4/27/05, RJF.  
;; more notes, 2/28/06 RJF. Converting h:/lisp/fft.lisp to QD arithmetic.
;; see comments at fft.lisp.
;; We use four1 only
;; the input data will be over-written by the result.
;; the inverse fft is indicated by :isign -1, though it must be
;; multiplied by  1/nn to work.
; functions:
;	four1: fourier transform (FFT) in one dimension

(defun v2dfa(a &optional (m (length a))	)
  ;;coerce a vector of QD numbers of length m to a QD array of length
  ;; 2m, since here complex numbers are stored in 2 adjacent QD
  ;; locations.
  (let* ((k (length a))
	 (ans (make-array (cl::* 2 m) ))
	 (zz (qd::into 0)))
    (dotimes (i k)			;just copy the actual numbers, or convert?
      (setf (aref ans (cl::* 2 i))(qd::into (aref a i))) ;; here we convert.
      (setf (aref ans (cl::1+ (cl::* 2 i))) (copy zz)))
    (loop for i from (* 2 k) to (1-(* 2 m)) do (setf  (aref ans i) (copy zz)))
    ans))

;; (v2dfa #(30 40 50) 4)
;;  --> #(0.3Q2 0.Q0 0.4Q2 0.Q0 0.5Q2 0.Q0 0.Q0 0.Q0)

(defparameter *zz* (qd::into 0))
(defparameter *one* (qd::into 1))

(defun dfa2v(a &optional (m (/ (length a)2)))
  ;; 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 (/ (length a) 2))
	 (ans (make-array m)))
    (dotimes (i m ans)
      (setf (aref ans i)(round (ga::outof (aref a (* 2 i))) k)))))

;;(dfa2v  (v2dfa #(256 256 256) 4))
;; -->  #(64 64 64 0)  is correct.

;; this works!
(defun polymultfftqd(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 ))) 
    (dfa2v(four1 prod z :isign -1) lans)))

;; Utility to copy an array because this fft will clobber input.
;; We don't use it here, but here's a way to write it.
;; (defun copyarray (a)(make-array (length a):initial-contents a))

(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))))
  (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 ;(type aqd a b c d)
	       (fixnum ind ind1))
      (setf (aref ans ind)(copy(with-temps (- (* a c)(* b d)))))
      (setf (aref ans ind1)(copy(with-temps (+ (* a d)(* b c)))))))))

(defparameter qdtwopi
    (ga::/ (qd::into  165424160919322423196824703508232170249081435635340508251270944637 )
     (qd::into  26328072917139296674479506920917608079723773850137277813577744384))
  )
(defparameter onehalf (/ (qd::into 1) (qd::into 2)))


;;0.62831853071795864769252867665590057683943387987502116419498891846Q1

;;; this is a fairly generic FFT that works, but not optimized much.
#+ignore
(defun four1 (data nn &key (isign 1))

 ;;
  (declare (type fixnum nn isign))
 (prog ((wr (copy *zz*)) 
	(wi (copy *zz*)) 
	(wpr (copy *zz*))
	(wpi (copy *zz*))
	(wtemp (copy *zz*)) 
        (theta (copy *zz*)) 
	(tempr (copy *zz*)) 
	(tempi (copy *zz*))
	(j 0) (n 0) (m 0) (mmax 0) (istep 0))
   (declare
	;;(type aqd wr wi wpr wpi wtemp theta tempr tempi)
    (fixnum j n m mmax istep))
   
  (setf n (* 2 nn)) 
  (setf j 1) 
  (do ((i 1 (+ i 2)))
      ((> i n) t)
      (declare (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 (cl::* 2 mmax))
    (setf theta  (/ qdtwopi (* isign mmax)))
    (setf wpr  (identity(* -2 (expt (sin (* 1/2 theta)) 2))))
    (setf wpi (sin theta)) (setf wr (copy *one*)) (setf wi (copy *zz*))
    (do ((m 1 (+ m 2)))
	((cl::> m mmax) t)
      (declare (fixnum  m))
      (do ((i m (+ i istep)))
	  ((> i n) t)
	(declare (type fixnum i))
	(setf j (+ i mmax))
	(setf tempr (identity (- (* wr (aref data (1- j)))
				    (* wi (aref data j)))))
	(setf tempi (identity (+ (* wr (aref data j))
				    (* wi (aref data (1- j))))))
	(setf (aref data (1- j)) (identity (- (aref data (1- i)) tempr)))
	(setf (aref data j) (identity (- (aref data i) tempi)))
	(setf (aref data (1- i)) (identity (+ (aref data (1- i)) tempr)))
	(setf (aref data i) (identity (+ (aref data i) tempi))))
      (setf wtemp wr)
      (setf wr (identity (+ (+ (* wr wpr) (* (* -1 wi) wpi)) wr)))
      (setf wi (identity(+  (+ (* wi wpr) (* wtemp wpi)) wi))))
    (setf mmax istep)
    (go label2)) 
   (return data)))

;; hacking away at it to make it faster for QD....
(defun four1 (data nn &key (isign 1))

  (declare (type fixnum nn isign))
  (prog ((wr (copy *zz*)) 
	 (wi (copy *zz*)) 
	 (wpr (copy *zz*))
	 (wpi (copy *zz*))
	 (wtemp 0) 
	 (theta (copy *zz*)) 
	 (tempr 0) 
	 (tempi 0)
	 (j 0) (n 0) (m 0) (mmax 0) (istep 0))
	(declare  (fixnum j n m mmax istep))
	(setf n (cl::* 2 nn)) 
	(setf j 1) 
	(do ((i 1 (cl::+ i 2)))
	    ((cl::> i n) t)
	  (declare (fixnum i))
	  (when (cl::> j i) 
	    (setf tempr (aref data (cl::1- j)))
	    (setf tempi (aref data j)) 
	    (setf (aref data (cl::1- j)) (aref data (cl::1- i)))
	    (setf (aref data j) (aref data i)) 
	    (setf (aref data (cl::1- i)) tempr)
	    (setf (aref data i) tempi))
	  (setf m (cl::floor  n 2))
  label1
	  (when (and (cl::>= m 2) (cl::> j m))
	    (setf j (cl::- j m)) (setf m (cl::floor  m 2))
	    (go label1))
	  (setf j (cl::+ j m))) 
	(setf mmax 2) 
  label2 
	(when (> n mmax)
	  (setf istep (cl::* 2 mmax))
	  (dsetv theta (into (cl::* isign mmax)))
	  (setf theta (with-temps (/ qdtwopi theta)))
	  (dsetv wpr (with-temps  (sin (* 1/2 theta))))
	  (dsetv wpr  (with-temps (* -2 (* wpr wpr))))
	  (dsetv wpi (sin theta))
	  (setf wr (copy *one*))
	  (setf wi (copy *zz*))
	  (do ((m 1 (cl::+ m 2)))
	      ((cl::> m mmax) t)
	    (declare (fixnum m))
	    (do ((i m (cl::+ i istep)))
		((cl::> i n) t)
	      (declare (fixnum i))
	      (setf j (cl::+ i mmax))
	      ;; next line works, but not with dsetv
	      (setf tempr (with-temps (- (* wr (aref data (1- j)))
					 (* wi (aref data j)))))
	      (setf tempi (with-temps (+ (* wr (aref data j))
					 (* wi (aref data (1- j))))))

	      (setf (aref data (cl::1- j))  (- (aref data (cl::1- i)) tempr))
	      (setf (aref data j)   (- (aref data i) tempi))
	      (setf (aref data (cl::1- i))   (+ (aref data (cl::1- i)) tempr))
	      (setf (aref data i) (+  (aref data i) tempi)))
	      
	    (setf wtemp wr)
	    (setf wr (copy (with-temps (+ (+ (* wr wpr) (* (* -1 wi) wpi)) wr))))
	    (setf wi (copy (with-temps (+ (+ (* wi wpr) (* wtemp wpi)) wi)))))
	  (setf mmax istep)
	  (go label2)) 
	(return data)))
)

#|  what the answer should be ...
(defun t1()(polymultfftqd #(1  2 3 4 5 6) #(7 8 9)));; test
(t1)
#(7 22 46 70 94 118 93 54)
(four1 (v2dfa #(1 2 3 4 5 6) 16) 16) ;; except higher precision
#(21.0d0 0.0d0 4.203712543852026d0 17.12548253340269d0
  -9.656854249492387d0 2.999999999999995d0 1.4918055861931159d0
  -4.857754915068683d0 2.999999999999997d0 4.0d0 ...)

(defun randpol(n m)
  (let ((ans (make-array n))
	(lim (expt 2 m)))
    (dotimes (i n ans)(setf (aref ans i)(random lim)))))

|#