;;; attempts by RJF to make OCT work fast in Allegro, based on 
;;; translating RTOY's code to macros, and converting code from RJF's QD to OCT.
;;; see comments in file octi.lisp too
;;;  certainly until I include all of RTOY's stuff for sin/cos/.
;;;  Now has fft, polynomial eval, quadrature

;;; 10/18/07 RJF

(defpackage :oct
  (:use :ga  :octi :cl)
    (:shadowing-import-from  
   :ga    
   "+" "-" "/" "*" "expt"		;... n-ary arith
   "=" "/=" ">" "<" "<=" ">="		;... n-ary comparisons
   "1-" "1+" "abs" "incf" "decf"
   "min" "max"
   "ftruncate" "ffloor" "fround" "fceiling" 
   sin cos tan 
   asin acos 
   atan log   ;; these are trickier, 1 or 2 args.
   exp
   sinh cosh tanh 
   asinh acosh atanh 
   sqrt 
   tocl 
   numerator denominator
   expt 
   rational  ;; convert to rational, exactly
   )
    (:export into)
    )

(in-package :oct)   ;;; NOT OCTI


(eval-when (:execute :load-toplevel :compile-toplevel :execute)
  (declaim (optimize (speed 3) (safety 1) (space 0) (compilation-speed 0)))
  (load "octi")
  
(defun into(x &optional (targ 0 targp)) ;if a target supplied, store into it
    (if targp (progn (octi::into x (oct-real targ))targ)
      (make-oct :real (octi::into x))))



(defstruct oct (real  #+allegro
		      (make-array 4 :element-type 'double-float 
			       :initial-element 0.0d0 
			       ;; next line is Allegro-specific
			       :allocation :lispstatic-reclaimable
			       )
		      #-allegro
		      (make-array 4 :element-type 'double-float 
			       :initial-element 0.0d0)

		      :type (simple-array double-float (4)))))

(defmethod make-load-form ((a oct)&optional environment)
  (make-load-form-saving-slots a :environment environment)

;;(defmethod make-oct-d ((a array))  ;; make an object from the array
;;  (make-instance 'oct :real a))
)
(defmacro copy_mac(a);;make a fresh copy of oct a.
  `(make-oct :real
   (make-array 4 :element-type 
			    'double-float 
			    :initial-contents (oct-real ,a)
			    :allocation :lispstatic-reclaimable)))

(defun copy_into(op1 op2) ;fill op2's memory with op1's value.
  (octi::copy-octi-into-octi (oct-real op1)(oct-real op2)) op2)

(defmacro defarithmetic (op pgm)
    (let ((two-arg
	   (intern (concatenate 'string "two-arg-" (symbol-name op))
		   :ga))
	  (oo-entry ;; two oct args
	   (intern (concatenate 'string (symbol-name pgm) "-oct-t")
		   :octi))
	  ;;(od-entry ;; oct and double args
	  ;; (intern (concatenate 'string (symbol-name pgm) "-oct-d-t") :octi))
	  )
	
 ;;           (format t "~% defining ~s" two-arg)
      `(progn
	 
	 ;; new defmethods for oct. ..
	 (defmethod ,two-arg ((arg1 oct) (arg2 oct))
	   (let* ((r (oct::make-oct))
		  (in (oct::oct-real r))
		  (a1 (oct::oct-real arg1))
		  (a2 (oct::oct-real arg2)))
	     (declare (optimize speed)
		      (type(simple-array double-float (4)) in a1 a2))
	     (,oo-entry a1 a2 in) r))
       
	 (defmethod ,two-arg((arg1 real) (arg2 oct))
	  #+ignore (,two-arg (into arg1) arg2)
	   (let* ((r (oct::make-oct :real (octi::into arg1)))
		  (in (oct::oct-real r))
		  (a2 (oct::oct-real arg2)))
	     (declare (optimize speed)
		      (type(simple-array double-float (4)) in a1 a2))
	     (,oo-entry in a2 in) r)
	   )
       
	 (defmethod ,two-arg ((arg1 oct) (arg2 real))
	  #+ignore    (,two-arg arg1 (into arg2))

		   (let* ((r (oct::make-oct :real (octi::into arg2)))
		  (in (oct::oct-real r))
		  (a1 (oct::oct-real arg1) ))
	     (declare (optimize speed)
		      (type(simple-array double-float (4)) in a1 a2))
	   
	     (,oo-entry a1 in in) r))
	 (setf (get ',op 'argnum) 2) ;used by with-temps, dsetv
	 (setf (get ',op 'oct-program) ',oo-entry) ;used by with-temps, dsetv
	 (setf (get ',two-arg 'oct-program) ',oo-entry) ;used after macroexpand-all
	 (setf (get ',two-arg 'argnum) 2)
	 
	 )))


(defarithmetic + add)
(defarithmetic * mul)
(defarithmetic - sub)
(defarithmetic / div)

;; need to do other elementary functions etc as well.
;; see lisp/generic/qd and  lisp/oct/qd-fun..

;; this special case doesn't fit into defarithmetic macro.
(defmethod ga::two-arg-expt ((base oct)(n integer))
  (let ((ans (make-oct)))
    (octi::pow-oct-i-t (oct-real base) n (oct-real ans))
    ans))


       
(defmethod print-object ((a oct) stream)(format stream "~a" 
						(octi::oct2string (oct::oct-real a))))

(defmethod octi::lisp2oct((s string)(ans array))(string2oct s))

(defmethod octi::lisp2oct((a oct) (ans array)) (octi::copy-octi-into-octi (oct-real a) ans ))

(defun oct-reader(stream subchar arg)
  (declare (ignore subchar arg))
  (let ((s (read stream)))

    (string2oct (format nil "~s" s))))

(defun string2oct(s);; read integerQinteger like 10Q2   = 100 = 0.1Q3  or  -3.1q-1
  ;; examples of acceptable numbers [put " " around them]
  ;; .1q0 ; 1Q0; 1.2q3;  -1.23q-4;  - .23 q -4;  2.3q+4; +4.5q6; q0; --3Q1 same as 3q1
  ;; 3 ; 3.4; 123.45q+100; -.q0 ; .q0 ;
  ;; examples of not acceptable numbers:
  ;;    q ; .q;   123.45q+100  [=NaN(qd)]
  (let ((p0 0)(sign 1))
    (setf s (delete #\  s));; remove leading or other spaces from string.
    (if (char= #\- (aref s 0))(setf p0 1 sign -1));; set sign if negative
    (multiple-value-bind (frac pos);; read fraction to left of .
	(parse-integer s :start p0 :radix 10 :junk-allowed t)
  ;;    (format t "~% frac= ~s pos=~s" frac pos)
      (if (null frac)(setf frac 0));; empty fraction is zero
      (if (= pos (length s)) (* sign (into frac))
	(case (aref s pos);; look at next char
	  ((#\Q #\q);; 10Q2
	   (multiple-value-bind (expon pos2)
	       (parse-integer s :start (1+ pos);skip the Q
			      :radix 10)
	;;     (format t "~% sign=~s frac=~s expon=~s" sign frac expon)
	     (* sign (into frac) (expt 10 expon))))
	  (#\.;; 
	   (multiple-value-bind (frac2 pos2)
	       (parse-integer s :start (1+ pos);skip the "."
			      :radix 10 :junk-allowed t )
	     (if (null frac2)(setf frac2 0))
	     (setf frac
		   (+ (into frac)
		      (* (if (< frac 0) -1 1) frac2 (cl::expt 10 (cl::1+(cl::- pos pos2))))))
	      ;; (format t "~% string pos ~s frac= ~s" pos2 frac)
	     (if (cl::= pos2 (length s))(* sign (into frac))
	       (case (aref s pos2)
		 ((#\Q #\q);; 10Q2
		  (multiple-value-bind (expon pos3)
		      (parse-integer s :start (1+ pos2);skip the Q
				     :radix 10 )
		    (* sign  frac (cl::expt 10 expon))
		    ))
		 (otherwise 
		  (format t "next char is ~a -- Not an oct spec: ~s"  
			  (aref s pos2) s)
		  (into (or frac 0)))))))
	  )))))



;; reading OCT numbers can be done this way...
(set-dispatch-macro-character #\# #\q #'oct::oct-reader) 

(defun time-mul (n)
  (declare (fixnum n)
	   (optimize (speed 3) (safety 1) (debug 0)))
  (let* ((h (into 1/3)))
    (print h)
    (time (dotimes (k n)
	    (declare (fixnum k))
	    (* h h)))))

(defun time-poly (n)
  (declare (fixnum n)
	   (optimize (speed 3) (safety 1) (debug 0)))
  (let* ((h (list (into 5)(into 2)(into 1)))
	 (arg (into 100)))
    (print h)
    (time (dotimes (k n)
	    (declare (fixnum k))
	    (polyeval h arg)))))
(defun time-polyd (n)
  (declare (fixnum n)
	   (optimize (speed 3) (safety 1) (debug 0)))
  (let* ((h (list 5d0 2d0 1d0))
	 (arg 100d0))
    (print h)
    (time (dotimes (k n)
	    (declare (fixnum k))
	    (polyevald h arg)))))



;;Compiler for  oct data type
;; RJF
(in-package :oct)

;;; We want to make better use of the state-based programs like
;;; mul-oct-t Assuming octs.. for a, b, and c: (dsetv a (+ b c))
;;; destroys the value in a.  Compare this to (setf a (+ b c)) which
;;; creates a new value and points a to it.

;; dsetv,  data driven
(defmacro dsetv (targ ex)
  ;; try  (dsetv a (+ b c)) 
  ;; should be faster than (setf a (+ b c)). maybe 2X.
  ;; All the logic below is done during macro-expansion,
  ;; which means it is usually done at compile time. Run time
  ;; is therefore not penalized.  If you use dsetv from an interpreted
  ;; program it will be slow, however, because it will do the macro
  ;; expansion followed by the execution, each time it is used.
  (setf ex (macroexpand ex))  
  (cond 
   ((atom ex) `(into ,ex ,targ))
   ((eq (car ex) 'into) `(into ,@(cdr ex)  ,targ))
   ((eq (car ex) 'setq) 
    (let ((gg (gensym))) ;; need to protect against capturing z in (setq z ..))
    `(let ((,gg  ,(with-temps (caddr ex))))
       (copy_into  (oct-real ,gg) (oct-real ,(cadr ex)))
       ,gg)))
   (t 
    (let* ((op (car ex))
	   (args (cdr ex))
	   (the-op (get op 'oct-program))
	   (argnum (get op 'argnum)))
      (cond 	       
       ((not the-op);; not a previously listed op
	`
	   (let* ((lval ,targ)
		  (a1 (oct-real  (,op ,@ args)))
		  (tt (oct-real lval)))
	     (declare (optimize speed)
		      (type (simple-array double-float (4)) a1 tt))
	     (copy_into a1 tt)
	     lval))
       ((not (eql argnum (length args))) 
	(error "dsetv was given operator ~s which expects ~s args, but was given ~s --  ~s" 
	       op argnum (length args) args))
       (t
	(case argnum
	  (1;; one argument.
	   `(let ((a1 (oct-real ,(macroexpand `(with-temps ,(car args)))))
		    (tt (oct-real ,targ)))
		(declare (optimize speed)(type (simple-array double-float (4)) a1 tt))
		;; could also check other args for being type qd
		;; could also allow for args to be si, ui, dd, etc.
		;; could also check number of args to be appropriate for operation
		(,the-op a1 tt)
	      ,targ))
	  (2
	   `(let ((a1 (oct-real ,(macroexpand `(with-temps ,(car args)))))
		  (a2 (oct-real ,(macroexpand `(with-temps ,(cadr args)))))
		    (tt (oct-real ,targ)))
		(declare (optimize speed)(type (simple-array double-float (4)) a1 a2 tt))
		(,the-op a1 a2 tt)
		,targ
		))
	  (otherwise (error "argnum is wrong for op ~s " op))
	  )))))))

;;;;;;;;;;;;;;;;;;more efficiency hackery follows.;;;;;;;;;;;;;;;;


;; We just allocate a few private "registers" say, for a
;; function, or an inner loop, and re-use them, if we are in a
;; loop. No need to tell anyone else about a few temp locations,
;; especially if they are GC'd when truly inaccessible.  That's what
;; is below.

(defmacro with-temps(expr)
  (let ((*names* nil)
	(*howmany* 0))
    (labels ((genlist(n)(loop for i from 1 to n collect (into i))) ;make a list of fresh qd items
	     (ct1 (r) ;; count temporaries needed
	       (cond ((numberp r) (incf *howmany*))
		     ((not (consp r)) r)
		     (t (incf *howmany*)
			(mapc #'ct1 (cdr r)))))
		
	   (maketemps(r) ;change r=(+ a (* b c)) to  temp storage .
		     (cond ((numberp r) (into r))
			   ((atom r) r)
			   ((get (car r) 'argnum); known operator
			    `(dsetv ,(pop *names*)
				    ,(cons (car r)(mapcar #'maketemps (cdr r)))))
			   ;; just a symbol name? maybe aref? better be the right type, aqd.
			   (t  r))))
      (setf expr (macroexpand expr))
       (ct1 expr)
     ;; (ct1 expr); count the temporaries
    (setf *names* (genlist *howmany*))
    (maketemps expr))))


;;  try (pprint (macroexpand '(with-temps (+ x (* 3 z)))))
;; or  (defun hypot(x y)(copy_mac (with-temps (sqrt (+ (* x x)(* y y))))))
;; need the call to copy to make a  copy of the result before calling hypot again. 

;; set up the environment for with-temps and dsetv here

(eval-when (compile load eval)
  (mapc #'(lambda(h) (setf (get h 'argnum) 2)) '(atan2 log2 setq))
;; for now assume they are given only one arg and if they are given 2 signal an error with dsetv.
  (mapc #'(lambda(h) (setf (get h 'argnum) 1)) '(atan log))
;;  (defun qd_setq (a b)(qd_copy_into b a) b)
 )

;; time trial

(defun time-mul-d (n)
  (declare (fixnum n)
	   (optimize (speed 3) (safety 1) (debug 0)))
  (let* ((h (into 1/3))
	 (ans (into 0)))
    (print h)
    (time (dotimes (k n)
	    (declare (fixnum k))
	    (dsetv ans (* h h))))))


;;;;;;;;;;;::::There's a lot more in generic/qd.lisp
;;;;;;;;;not yet converted to OCT.

(defmacro defcomparison (op)
    (let ((two-arg-ga (intern (concatenate 'string "two-arg-" 
					   (symbol-name op)) :ga))
	  (two-arg-octi (intern (concatenate 'string "two-arg-" 
					(symbol-name op)) :octi))
	  )
    ;;  (format t "~% defining ~s" two-arg)
      `(progn
	 ;; only extra methods not in ga are defined here.
	 ;; qd_comp returns -1 0 1   for < = >
	 (defmethod ,two-arg-ga ((arg1 oct) (arg2 oct))   
	   (,two-arg-octi (oct-real arg1)(oct-real arg2)))

	 (defmethod ,two-arg-ga ((arg1 real) (arg2 oct))
	   (let ((arg1 (into arg1)))
	     (,two-arg-octi (oct-real arg1)(oct-real arg2))))
	 
	 (defmethod ,two-arg-ga ((arg1 oct) (arg2 real))
	   (let ((arg2 (into arg2)))
	     (,two-arg-octi (oct-real arg1)(oct-real arg2))))
	 ',op)))

(defcomparison >)
(defcomparison <)
(defcomparison =)
(defcomparison /=)
(defcomparison >=)
(defcomparison <=)

;; allows (< (into 1)(into 2)(into 3)) 



;;How about the single-arg functions like sin cos log...

(defmacro r (op);;
  (let ((fun-name  (intern op :ga ))
	(octi-name (intern (concatenate 'string op "-oct-t") :octi)))

    `(progn
       (defmethod ,fun-name ((arg oct))
	 (let* ((h (make-oct)) (in (oct-real h)))
	   (declare (optimize speed)
		    (type (simple-array double-float (4)) in))
	   (,octi-name (oct-real arg) in)
	   h))
       
       (setf (get ',fun-name 'argnum) 1)
       (setf (get ',fun-name 'oct-program) ',octi-name)
       ',fun-name
       )))


;; (r "sqrt")   take care of complex version. see later.

(r "abs")
(r "sin")(r "cos")

;;(r "log")(r "exp") ;; etc
(r "1-")
(r "1+")
(r "exp")


(defmethod ga::one-arg-atan((arg oct))
           (let* ((h (make-oct)) (in (oct-real h)))
           (declare (optimize speed)
            (type (simple-array double-float (4)) in))
           (octi::atan-oct-t (oct-real arg) in)
           h
	   ))

(defun ftruncate ( x &optional (y 1))
  (octi::ftruncate-oct (oct-real x) y))
(defun fceiling (x &optional (y 1)) ;; fceiling for oct only
  (octi::fceiling-oct (oct-real x) y))
(defun fround (x &optional (y 1))
  (octi::fround-oct (oct-real x) y))

;; complex version
;; we do not propose to implement complex numbers at the octi level.


(defstruct octz (real (into 0) oct)
	   (imag (into 0)oct))

;; try these out.

(defmethod ga::two-arg-+ ((a octz)(b octz))
  (let ((ans (make-octz)))
    (setf (octz-real ans) (+ (octz-real a)(octz-real b))
	  (octz-imag ans) (+ (octz-imag a)(octz-imag b)))
    ans))

(defmethod ga::two-arg-* ((a octz)(b octz))
  (let ((ans (make-octz)))
    (setf (octz-real ans) (* (octz-real a)(octz-real b))
	  (octz-imag ans) (* (octz-imag a)(octz-imag b)))
    ans))

(defun lisp2octz(a)
  (let ((ans (make-octz)))
    (setf (octz-real ans) (into (realpart a))
	  (octz-imag ans) (into (imagpart a)))
    ans))

;; we could set up a separate package for complex oct..
;; but then the sqrt(neg) issue becomes complicated.

(defmethod ga::exp ((z octz))
  ;; exp(a+bi)= exp(a)*(cos b + i sin b)
  
  (let* ((ans (make-octz))
	 (a (octz-real z))		;real
	 (b (octz-imag z))		;real
	 (ea (exp a))			;real
	 )
    (setf (octz-real ans)   (* ea(cos b))
	  (octz-imag ans)   (* ea (sin b)))
    ans))

;; re-implement ga::sqrt((z oct))  as well as ga::sqrt((z octz)).


(defmethod sqrt((x oct))
  (if (< (octi::oct-0 (oct-real x)) 0d0) ;; negative answer.
      (make-octz :imag (sqrt (- x)))
     (make-oct :real (octi::sqrt-oct (oct-real x)))))
    
;;; need to allow for octz+oct, oct+octz, octz+real real+octz  etc.
;;; whew, macro writing can do all this too
  
;;;;;; some fun.
;; need abs, < >
#|

;; A BUNCH OF PROGRAMS PLAYING OCT GAMES and ROOTFINDING.

;; Refining a polynomial root via Newton's method.
;; First, deriv computes the derivative of a polynomial, in a list.
;; recall that we can create a polynomial as a list:
;; (setf testpoly (list 5 3 1)) ;; 5*x^2+3*x+1
;; or as an OCT polynomial
;; (setf testpoly (list (into 5) (into 3) (into 1))) ;; 5*x^2+3*x+1
|#

(defun deriv(coefs)
  ;;given coefs of a polynomial. return coefs of derivative.
  ;; constant coeff is last.  Any kind of numbers are OK.
  (let ((ans nil) (i 0))
    (dolist (c (cdr (reverse coefs)) ans) 
      (push (* (incf i) c) ans))  ans))

(defun polyeval (coeflist x) ;coeflist of octs, x is an oct
      (let ((sum (into 0)))
      (dolist (i coeflist sum)
	(setf sum (with-temps (+ i (* x sum)))))))

;; example  (polyeval (list (into 1d0)(into 2d0)(into 3d0))(into 2d0))

;;  a fairly general Newton iteration requiring three parameters, and some optional ones.
;;  f  is a function to evaluate f(x),
;;  df is a function to evaluate f'(x).
;;  x, a starting point,
;; optional: threshold for absolute value of f at which to stop,
;; optional: iters, a maximum number of iterations.

;;  oneroot uses whatever arithmetic is appropriate, based on what f and df use.

(defun oneroot (f df x threshold iters &aux fval)
  (dotimes  (i iters (error "rootfinder failed to converge. Residual is ~s after ~s iterations." 
			    fval i))
    (setf fval (funcall f x))
    (if (< (abs fval) threshold) 
	(return (values x fval i))		;return x, residual and iteration count
      (decf x (/ fval (funcall df x))))))


;; if we replace (decf x (/ fval (funcall df x)))
;; by            (dsetv x (with-temps (- x (/ fval (funcall df x))))))))
;; this would save 2 allocated temporaries per iteration.


;; A carefully declared polynomial eval using
;; Lisp double-floats. I don't see how to make this more optimized.
(defun polyevald (dlist x) 
  (let ((sum 0.0d0))
    (declare (double-float sum x) (optimize (speed 3)(debug 0)))
    (dolist (i dlist sum)
      (declare (double-float i))
      (setf sum (the double-float 
		  (cl::+ i
			 (the double-float 
			   (cl::* x sum)))))))) ;; this is about 130 X faster than the oct version.
;; example  (polyevald '(1d0 2d0 3d0)2d0) is 11d0

;; Next, set up a call to oneroot given a list of coefficients representing
;; a polynomial.  polyrootd makes everything computing in (real) double-float.
;; defaults are provided, so all you need is the list and a starting point.
;; the inputs coefs and x may be given as exact integers, floats, rationals.
;; example  (polyrootd '(1 0 -2) 1); returns 1.4142135623730951d0, i.e. sqrt(2)

(defun polyrootd (coefs x &optional (threshold 1.0d-14) (iters 20))
  (setf coefs (mapcar #'(lambda(z)(coerce z 'double-float))coefs))
    (let ((dp (deriv coefs)))
      (oneroot  #'(lambda(x)(polyevald coefs x))
		#'(lambda(x)(polyevald dp x))
		(* 1.0d0 x)
		(* 1.0d0 threshold)
		iters)))

;; this is the same functionality, but using higher (QD) precision.
;; the inputs coefs and x may be given as exact integers, floats, rationals.
;; (oct::polyroot '(1 0 -2) 1); returns sqrt(2), but if you want more precision:
;; (oct::polyroot '(1 0 -2) 1 1.0d-60)  returns
;;0.141421356237309504880168872420969807856967187537694807317667973796Q1

;;this is called oct:polyroot. You may prefer other ways to do it, 
;;e.g. mpfr::polyroot or some interval method, or some method
;; computing derivatives.

(defun polyroot (coefs x &optional (threshold 1.0d-50) (iters 20))
  (setf coefs (mapcar #'into coefs))
  (let ((dp (deriv coefs)))
    (oneroot  #'(lambda(x)(polyeval coefs x))
	      #'(lambda(x)(polyeval dp x))
	      (into x)
	      (into threshold)
	      iters)))

;; try (polyroot '(1 0 -4) 1 1d-60 20)

;; ONEROOT is considerably more versatile: it can find root of other
;; functions, not just polynomials.
#|
(defun tt (x)(- (cos x) x))  ;; easy near 0.73
(defun ttd(x) (- (+ (sin x) 1))  )
(oneroot #'tt #'ttd (into 0) (into 1.0d-64) 20)  ;; uses oct
(oneroot #'tt #'ttd 0.0d0 1.0d-38 20) ;;uses double-float.

|#
	

;;;;;;;;;;;;;;;;;;;FFT!!!;;;;;;;;;;;;;;;;;

(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 OCT arithmetic.
;; slightly rearranged for OCT, 10/29/07  
  
;; see comments at lisp/fft.lisp.  We use four1 only.  The input data
;; will be overwritten 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
  
(defparameter *zz* (into 0))
(defparameter *one* (into 1))

(defun v2dfa(a &optional (m (length a))	)
  ;;coerce a vector of oct numbers (or lisp numbers) of length m to a
  ;; oct array of length 2m, since here complex numbers are stored in 2
  ;; adjacent oct locations. Caution: same objects are being used. Do not
  ;; mutate them.
  (let* ((k (length a))
	 (ans (make-array (cl::* 2 m) :allocation :lispstatic-reclaimable ))
	 (zz (copy_mac *zz*));zero
	 (h nil))
    (declare (fixnum k m)(optimize (speed 3)(safety 1)))
    (dotimes (i k)	;just point to the actual numbers, or convert if needed
      (declare (fixnum i))
      (setf (aref ans (cl:* 2 i))
	(if (oct-p (setf h (aref a i))) h (oct:into h))) ;; here we convert.
      ;;(setf (aref ans (cl:1+ (cl:* 2 i))) (copy_mac zz))
      (setf (aref ans (cl:1+ (cl:* 2 i))) zz )); just use zero
    (loop for i fixnum from (* 2 k) to (1-(* 2 m)) do 
	  (setf  (aref ans i) zz)) 
    ans))

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


(defmethod ga::outof ((x oct))(oct2lisp x))

(defmethod oct2lisp((x oct))
  (if (excl::nan-p (aref (oct-real x) 0)) (aref (oct-real x) 0)
					; if qd contains a NaN, use it.
      (apply #'cl:+  (map 'list #'rational (oct-real x)))))

(defun dfa2v(a &optional (m (/ (length a)2)))
  ;; Coerce real parts back
  ;; to integers, more or less.  a is 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)))
    (declare (fixnum k m)(optimize (speed 3)(safety 1)))
    (dotimes (i m ans)
      (declare (fixnum i))
      (setf (aref ans i)(round (ga::outof (aref a (cl::* 2 i))) k)))))

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

;; this works!
(defun polymultfft(r s)
  (declare (optimize (speed 3)))
  ;;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 (cl:+ lr ls -1))
	 (z (ash 1 (ceiling (cl:log lans 2)))) ; round up to power of 2
	 (rfft (four1 (v2dfa r z) z))
	 (sfft (four1 (v2dfa s z) z))
	 (ans (make-array (* 2 z))))
    (declare (fixnum z))
    (dotimes (i (* 2 z))(setf (aref ans i) (copy_mac *zz*)))
    (prodarray rfft sfft z ans)
    (dfa2v(four1 ans z :isign -1) lans)))

(defun polysquare(r) 
  ;;just square the polynomial r. saves one fft. so you can compare times.
  (let* ((lr (length r))
	 (lans (+ lr lr -1))
	 (z (ash 1 (ceiling (log lans 2)))) ; round up to power of 2
	 (rfft (four1 (v2dfa r z) z))
	 (ans (make-array (* 2 z))))
    (prodarray rfft rfft z ans)
    (dfa2v(four1 ans z :isign -1) lans)))

#+ignore
(defun polytime(r s)  ;; this shows that much of the time is in v2dfa etc
  (let* ((lr (length r))
	 (ls (length s))
	 (lans (+ lr ls -1))
	 (z (ash 1 (ceiling (log lans 2)))) ; round up to power of 2
	 (r1 (v2dfa r z))
	 (s1 (v2dfa s z))
	 (sfft 0)
	 (rfft  (four1 r1 z)))
    (start-profiler)
  (time (progn (setf sfft (four1 s1 z))
    (setf sfft (four1 s1 z))
    (setf sfft (four1 s1 z))
    (setf sfft (four1 s1 z))
    (setf sfft (four1 s1 z))))
    (show-flat-profile)
	 ;;(prod (prodarray rfft sfft z ))
	 ;;(ans (time(four1 prod z :isign -1)))
	 )
    ;;(dfa2v ans lans)
    )

#+ignore
(defun polytime2(r)  ;; this shows that much of the time is in v2dfa etc
  (let* ((lr (length r))
	 (lans (+ lr lr -1))
	 (z (ash 1 (ceiling (log lans 2)))) ; round up to power of 2
	 (r1 (v2dfa r z))
	 (start-profiler)
	 (rfft (time (four1 r1 z)))
	 (prod (prodarray rfft rfft z  rfft))
	 (ans (four1 prod z :isign -1)))
    (show-flat-profile)
  (dfa2v ans 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)			; of octs
;;  (map 'array #'(lambda(r)(copy_mac r)) a))

(defun prodarray(r s len ans)
  ;; 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
  (declare (fixnum len)(optimize (speed 3)(safety 1)))
  (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 i 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 oct2pi (make-oct :real octi::+oct-2pi+))
(defparameter octpi/2 (make-oct :real octi::+oct-pi/2+))

(defparameter onehalf (into 1/2))

;;0.62831853071795864769252867665590057683943387987502116419498891846Q1

;;; this is a fairly generic FFT that works, but not optimized much.
;;; we leave it here just in case you want to copy it for other 
;;; generic arithmetic packages
;;#+ignore  
;; this works, but is not optimized 
(defun copy(x)(if (oct-p x)(copy_mac x)(copy-tree x))) ;; whatever..
#+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  (/ oct2pi (* isign mmax)))
    (setf wpr  (* -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  (- (* wr (aref data (1- j)))
				    (* 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) 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) (* (* -1 wi) wpi) wr))
      (setf wi (+  (* wi wpr) (* wtemp wpi) wi)))
    (setf mmax istep)
    (go label2)) 
   (return data)))

;; hacking four1 for speed, keeping space consumption down

(defun four1 (data nn &key (isign 1))
  (declare (type fixnum nn isign) (optimize (speed 3)(safety 1)))
  (prog ((wr (copy_mac *zz*)) 
	 (wi (copy_mac *zz*)) 
	 (wpr (copy_mac *zz*))
	 (wpi (copy_mac *zz*))
	 (wtemp (copy_mac *zz*)) 
	 (theta (copy_mac *zz*)) 
	 (halftheta (copy_mac *zz*)) 
	 (cost (copy_mac *zz*))
	 (tempr (copy_mac *zz*)) 
	 (tempi (copy_mac *zz*))
	 (one (copy_mac *one*))
	 (zero (copy_mac *zz*))
	 (temprx 0) (tempix 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))
	  (declare (fixnum i))
	  (when (cl:> j i) 
	    (setf temprx (aref data (cl:1- j)))
	    (setf tempix (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)) temprx)
	    (setf (aref data i) tempix))
	  (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 (cl:> n mmax)
	  (setf istep (cl:* 2 mmax))
	  (dsetv theta (into (cl:* isign mmax)))
	  (dsetv theta (/ oct2pi theta))
	  (dsetv halftheta (* 1/2 theta))
	  (octi::sincos-oct-t (oct-real halftheta)(oct-real wpr)(oct-real cost))
	 (dsetv wpi  (* 2(* wpr cost))) ;; 2*sin(t/2)*cos(t/2)
	 (dsetv wpr  (* -2 (* wpr wpr)))
	 (dsetv wr one)
	 (dsetv wi zero)
	  (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))

	      (dsetv tempr  (- (* wr (aref data (cl:1- j)))
			       (* wi (aref data j))))
	      
	      (dsetv tempi (+ (* wr (aref data j))
			      (* wi (aref data (cl:1- j)))))

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

#|  what the answer should be ...
(defun t1()(polymultfft #(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 ...)  ;; works with unoptimized..

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

(setf r (randpol 512 332))
(time (progn (polymultfft r r) nil))
(time (progn (polymultfft qr qr) nil))
(setf s (make-array 512 :initial-element 1))
(time (polymultfft s s) ))
;; This is similar in speed to wxmaxima in GCL:
;; s:rat(sum(x^i,i,0, 511),x)$
;; s*s$  ;; time is .6 sec, vs .8 sec for fft.

|#


;;; quadrature
;;; Gaussian Quadrature, quad-double precision / OCT
;;; See quad-ga.lisp for basic ideas and comments
;;; quad-oqd.lisp for older versions.
;;; Jan 10, 2007.
;;; 10/29/2007  
;;;  by Richard Fateman

;; Legendre_pd returns legendre_p(k,x) and its derivative.
;; Note that d/dx P[n](x) =  (1/(x^2-1))*n*(x*P[n](x)-P[n-1](x)) except for x=+-1
;; This program is used with Newton iteration to refine roots of P[n](x).

;;(in-package :oct)



(defun legendre_pd (k x)
  (declare(double-float x) (fixnum k) (optimize (speed 3)(safety 0)))
  (assert (and (typep k 'fixnum)(<= 0 k)))
  (case k 
    (0 (values 1.0d0 0.0d0))
    (1 (values x 1.0d0))
    (otherwise
     (let ((t0 1.0d0)
	   (t1 x)
	   (ans 0.0d0))
       (declare(double-float t0 t1 ans x))
       (declare (fixnum k))
       (loop for i from 2 to k  do
	     (setf ans     (cl:/ (cl:- (cl:* (cl:1- (* 2.0d0 i)) x t1) (cl:* (cl:1- i) t0)) i)
		   t0 t1
		   t1 ans))
       (values        t1
		      ;; (1/(x^2-1))*k*(x*P[k](x)-P[k-1](x))
		      ;; except if abs(x)=1 then use
		      (if (= 1 (abs x))
			  ;; 1/2*k*(k+1)*x^(k+1)
			  (/ (* k (1+ k) (if (oddp k) 1 x)) 2)
			(cl:/ (cl:* k (cl:- (cl:* x t1)t0))  (cl:1-  (cl:* x x)))))))))

;;octfloat version
(defun legendre_pdbf (k x)
  (assert (and (typep k 'fixnum)(<= 0 k)))
  (case k 
    (0 (values (into 1)(into 0)))
    (1 (values (into x) (into 1)))
    (otherwise
     (let ((t0 (into 1))
	   (t1 (into x))
	   (ans (into 0))
	   (i (into 1)))
       (declare (fixnum k))
       (loop for ii from 2 to k  do
	     (dsetv i  (+ 1d0 i))
	     (setf ans  (with-temps(/(- (* (- (* 2 i)1) x t1) (* (- i 1) t0)) i)))
	     (dsetv t0 t1)
	     (dsetv t1 ans))
       (values        t1
		      ;; (1/(x^2-1))*k*(x*P[k](x)-P[k-1](x))
		      ;; except if abs(x)=1 then use
		      (if (= 1 (abs x))
			  ;; 1/2*k*(k+1)*x^(k+1)
			   (/ (* k (1+ k) (if (oddp k) 1 x)) 2)
			(with-temps (/ (* (into k) (- (* x t1)t0))  (1-  (* x x))))
			))))))

;; Here is a Newton iteration to converge to a root of legendre_p, the float version.
(defun improve_lg_root(n guess)
  (multiple-value-bind (val deriv)(legendre_pd n guess)
    (if (= deriv 0) (if (= val 0) 
			val 
		      (error "newton iteration failed at ~s"guess))
      ;;new_guess= guess-f(x)/f'(x)
      (- guess  (/ val deriv)))))

;; Here is a Newton iteration to converge to a root of legendre_p, bigfloat version.
(defun improve_lg_rootbf(n guess)
	   (multiple-value-bind (val deriv)(legendre_pdbf n guess)
	   ;;new_guess= guess-f(x)/f'(x)
	   ;; dsetv alters guess value in place
	   (dsetv guess  (- guess (/ val deriv)))) )


;; where x= a zero of P[k], m_wt computes the corresponding Legendre weight.
;; Lots of different forms for this; perhaps this is more stable than many others.

(defun m_wt (k x &aux (kf (into k)))		; compute -2 / (  k+1) * P'[k](x)*P[k+1](x))
  (assert (and (typep k 'fixnum)(< 1 k))) ; k>=2
  (setf x (into x)) 
	(let 
	    ((t0 (into 1))
	     (t1 x)
	     (ans 0)
	     (i (into 0))
	     (diffpk nil)(pkp1 nil))
	  (declare (fixnum k))
	  (loop for ii from 2 to k  do
		(setf i (into ii))
		(setf ans   
		  ;; make a new number here
		    (/  (with-temps(- (* (1- (* 2 i)) x t1) (* (1- i) t0))) i)
		      t0 t1
		      t1 ans))
	  ;; t0  is now legendre_p[k-1](x)
	  ;; ans is now legendre_p[k](x).
	  ;;deriv of p[k]  is k/(1-x^2) * p[k-1](x)-x*p[k](x)
	  (setf diffpk (with-temps (/ (*  kf (- t0 (* x t1))) (- 1 (* x x)))))
	  ;; compute  p[k+1](x)
	  (setf kf(+ kf 1))
	;;  (format t "~%kf=~s" kf)
	  ;; (setf pkp1  (/ (- (* (- (* 2 kf) 1) x t1) (* (- kf 1) t0)) kf)); 
	   (setf pkp1 (with-temps (/ (- (* (1- (* 2 kf)) x t1) (* (1- kf) t0)) kf)) )
;;	  (format t "~% diffpkm1=~s  pk+1=~s kf=~s" diffpkm1 pkp1 kf)
	 ;; (list diffpkm1 pkp1 kf (/ -2 (*(* kf diffpk) pkp1)))
	  (/ -2 (* kf diffpk pkp1))))

#|
:pa :oct
;; running the integration formula repeatedly is very cheap after
;; abscissae and weights are computed.  This  version computes
;; these items once as needed, then remembers them. Timing is much
;; faster after the first run.
(time (int_gl #'exp 20))
should give almost the same answer as
(defun right ()(let ((k (exp (into 1)))) (- k (/ 1 k)))) ; (right) is answer to above integral
|#

;; version 4
;;  to save allocation time and space we 
;;   0. initialize s to 0. For i= 0 to n/2 do lines 1-5:  ;; assume n is even.
;;   1. compute one abscissa x=x[i]
;;   2. compute one weight   w=w[i]
;;   3. compute two values    v=f(x)+f(-x))
;;   4. set s:= s+w*v.
;;   5. reusing x,w,v storage, repeat until n.
;; depending on whether n was even or odd, there may be one more point, at f(0)

;; Another improvement.  For the middle weight of an odd
;; number of terms, realize that the sum of all the weights will be
;; 2. So if we compute S, the sum weights up to n/2, the remaining
;; weight will be 2*(1 - S).

;; Also we can perhaps skip the full weight computation if n is even,
;; since one can be deduced by subtraction.

;; here's where we store precomputed wts and abscissae
(defparameter *legendabs* (make-hash-table))
(defparameter *legendwts* (make-hash-table))

;; integrate using gauss-legendre from -1 to 1
(defun int_gl(fun n)		; memoize the abscissae and weights
  (let ((abscissae (gethash n *legendabs*))
	(weights nil))

    (cond (abscissae (setf weights (gethash n *legendwts* )))
	  (t;; not found, so compute and remember them now.
		   (multiple-value-setq (abscissae weights)(ab_and_wts n))
	   (setf (gethash n *legendabs*) abscissae)
	   (setf (gethash n *legendwts*) weights)))
    ;; compute the sum.
    (let ((sum (into 0))	 
	  (halfn (ash n -1)))
      (loop for i from 0 to (1- halfn) do 
	    (dsetv sum 
		   (+ sum (*(+ (funcall fun (aref abscissae i))
			       (funcall fun (with-temps (-(aref abscissae i)))))
			    (aref weights i))))	 )
      (if (oddp n) 
	  (dsetv sum (+ sum (* (aref weights halfn)(funcall fun (into 0))))) )
      sum)))

;; pre-compute abscissae and weights
(defun ab_and_wts (n) ;precision is somewhat less than 215, qd. Stores only half the values (maybe +1)
  (let ((a 0)
	(v 0.0d0)
	(nph (/ 1.0d0(+ n 0.5d0)))
	(halfn (ash n -1))
	(abscissae (make-array (ceiling n 2)))
	(weights  (make-array (ceiling n 2))))
    (loop for i from 0 to (1- halfn) do
	  (setf v (cos(* pi (- (cl:1+ i) 0.25) nph)))
	  ;;	  (format t "~% compute cos(~s )=~s" (* pi (- (cl:1+ i) 0.25) nph) v)
	  (loop for i from 1 to 4  do	; 3 is almost enough.
		(multiple-value-bind (val deriv)(legendre_pd n v)
		  (declare (double-float val deriv))
		   ;;new_guess= guess-f(x)/f'(x)
		  (setf v (cl::- v (cl::/ val deriv)))))
	  ;; 	  (format t "~% --- improved v=~s" v)
	 (setf  a (into v))
        ;; compute more accurate abscissa.
	  (dotimes (k 4)
	  (multiple-value-bind (val deriv)(legendre_pdbf n a)
	    ;;new_guess= guess-f(x)/f'(x)
	    ;;  (format t "~% --- improved a=~s" a)
	    ;; dsetv alters a in place
	     (dsetv a  (- a (/ val deriv))) ))
	  
	  (setf (aref abscissae i) a)
	  (setf (aref weights i) (m_wt n a)))
    (cond ((oddp n)
	   (setf (aref abscissae  halfn) (into 0))
	   (setf (aref weights  halfn) 
      ;; fill in the middle element of an odd list by subtracting the others off
	     (let ((sum (into 0)))
		 (loop for i from 0 to (1- halfn)  do (dsetv sum (+ sum (aref weights i ))))
		 (dsetv sum (- 2 (+ sum sum)))))	   ))
    (values abscissae weights)))

(defun clear-leg-hash()			; call to deallocate storage from hash tables for legendre polys
  (clrhash *legendabs*)
  (clrhash *legendwts*))


(defun gaussab (fun lo hi n) ;; integrate from lo to hi
  (let ((a (* 1/2 (- hi lo)))
	(b (* 1/2 (+ hi lo)))
	)
    (* a (int_gl #'(lambda(r)(funcall fun (+ b (* a r))))
		 n))))

(defun gauss0inf (fun n) ;; integrate from 0 to inf, not a great way to sample.
  (int_gl #'(lambda(x)
		 (let ((d (- 1 x)))
		   (/ (funcall fun (/ (+ x 1) d)) (* d d))))
	  n))

;; Here is the  tanh/sinh method. This also integrates from -1 to 1.

(defun quadts(fun n)
  (declare (special octpi/2))
  (let*(
	(h (into(/ 4 n)))
	(sum (* octpi/2 (funcall fun 0)))
	(he (exp h))
	(relerr(into #.(expt 2.0d0 (- 50))))
	(t2 1)t3 t4 ab cor (oldsum (into 0d0)))
    (loop for j from 1 to (* 2 n) do
	  (setf t2 (* t2 he))
	  (setf t3 (exp (* octpi/2 (/ (- (* t2 t2) 1)
				    (* t2 2)))))
	  (setf ab (* t3 t3))
	  (setf t4 (/ (+ ab 1)(* 2 t3)))
	  (setf ab (/ (- ab 1)(* 2 t3 t4)))
	  (setf cor 
		(*(/ (* octpi/2 (+ (* t2 t2) 1))
		     (* 2 t2 t4 t4))
		  (+ (funcall fun ab)(funcall fun (- ab)))))
	  (setf oldsum sum)
	  (setf sum (+ sum cor)) 
	  (cond ((>= j n)
		 (if (= oldsum 0) (if (= sum 0) (return 0)))
		 (if (<(abs(/ (- sum oldsum) 
			      (if (= 0 oldsum)sum oldsum )))
		       relerr)
		     (return (* sum h))))))
    ;; fell out of loop
    (* sum h)))

(defun quadtsab (fun lo hi n) ;; integrate from lo to hi
  (let ((a (* 1/2 (- hi lo)))
	(b (* 1/2 (+ hi lo))))
    (* a (quadts #'(lambda(r)(funcall fun (+ b (* a r))))
		 n))))

#| some tests

(defun test2 ()
  ;;    ;;pi via Machin's atan formula~
  ;; pi/4 = 4 * arctan(1/5) - arctan(1/239)
  (* (- (* (atan (into 1/5)) 4) (atan(into 1/239))) 4))

(defun test3 ()
  (declare (optimize (speed 3)))
  ;; Salamin-Brent Quadratic formula for pi
  (let* ((a #q1)
	 (b (sqrt #q.5))
	 (s #q.5)
	 (m 1d0)
	 (p (/ (* (* a a)
		  2d0)
	       s)))
    (declare (double-float m))
    (dotimes (k 9)
      (setf m (* 2 m))
      (let* ((a-new (* (+ a b) .5d0))
	     (b-new (sqrt (* a b)))
	     (s-new (- s
		       (* (- (* a-new a-new)
			     (* b-new b-new))
			  m))))
	(setf a a-new)
	(setf b b-new)
	(setf s s-new)
	(setf p (/ (* (* a a) 2d0)
			s))))
    (format t "~2&Salamin-Brent Quadratic formula for pi~%")
  ;  (print-result p +pi+)
    p))

|#