;;; attempt to make bisect work with maxima BIGFLOATS, based on bisect-mpfr based on f2cl version
;;; broken.
;;; maybe allow mixed mode arithmetic
;;; ("f2cl1.l,v 1.222 2010/10/08 03:05:30 rtoy Exp $"

(in-package :maxima)
(eval-when (compile)(proclaim '(special fpprec diag offdiag)))

(defun tomp(x)(if (equal 0 x) (list 0 0)(cdr($bfloat x)))) ; fast if x is fixnum, double,rational..

;; fdo mirrors the fortran DO. The variable remains defined after
;; loop exit. Only one variable is stepped.
;; copied from f2cl macro file.
(defmacro fdo (do_vble_clause predicate_clause &rest body)
  (let ((step (gensym (symbol-name '#:step-)))
	(iteration_count (gensym (symbol-name '#:cnt-)))
	(loop-var (first do_vble_clause)))
    `(prog* ((,step ,(third (third do_vble_clause)))
	     (,iteration_count 
	      (max 0 (the fixnum
		       (truncate (the fixnum
				   (+ (the fixnum (- ,(third (first predicate_clause))
						     ,(second do_vble_clause)))
				      ,step))
				 ,step))
		   )))
	    (declare (type fixnum ,step ,iteration_count))
       ;; initialise loop variable
	    (setq ,loop-var ,(second do_vble_clause))
	    loop
	    (return
	     (cond			; all iterations done
	      ((zerop ,iteration_count) nil)
	      ;; execute loop, in/de-crement loop vble and decrement cntr
	      ,(cons 't 
		     (append 
		      (append body
			      `((setq ,loop-var (the fixnum ,(third do_vble_clause))
				      ,iteration_count (the fixnum (1- ,iteration_count)))))
		      '((go loop)))))))))

;; these could be macros
(defun fp> (x y)(fpgreaterp x y))
(defun fp<= (x y)(not (fpgreaterp x y)))
(defun fpzerop (x)(and (consp x)(equal (car x) 0)))
(defun fpminusp(x)(and (consp x)(< (car x) 0)))
(defun fpposp(x)(and (consp x)(> (car x) 0)))
(defun fpnonnegp(x)(and (consp x)(>= (car x) 0)))

;;set up the eigenvalue programs

(let ((*last-bigfrprec-bisect* 0)(*bigfr-machep*))
  
  ;; the bigfloat arrays are NOT full bigfloats. 
  ;; in lisp they look like 
  ;; ((BIGFLOAT SIMP 56) 36028797018963968 2).
  ;; we just use the CDR of these guys.
  
  (defun bisect (n eps1 d e e2 lb ub mm m w ind rv4 rv5)
    (declare (special fpprec)
	     (type (array t (*)) rv5 rv4 w e2 e d)
	     (type fixnum ierr m mm n)
	  ;   (type mpfr ub lb eps1)
	     (type (array fixnum (*)) ind))
    (prog 
     ((u (tomp 0)) (v (tomp 0)) (t1 lb) 
      (t2 ub) (xu (tomp 0)) (x0 (tomp 0))
      (x1 ub) (s1 (tomp 0)) (s2 (tomp 0)) 
      (half (tomp 0.5)) (two (tomp 2))
      (i 0) (j 0) (k 0) (l 0) (p 1)
      (q n) (r 0) (s 0) (ii 0) (m1 0) (m2 0) 
      (tag 0) (ierr 0)(isturm 0) (prec fpprec) (machep nil))

     (declare
      (type fixnum isturm tag m2 m1 ii s r q p l k
	    j i prec))
      ;; compute machine epsilon as 
      ;; (expt 2 (1- (mpfr::get-prec)))
      ;; unless we already did that.
     (cond((= prec *last-bigfrprec-bisect*)) ;already set?
	  (t (setf *last-bigfrprec-bisect* prec)
	     (setf *bigfr-machep* (fpexpt two (- 1 prec)))))
      (setf machep *bigfr-machep*)
      
     (fdo  (i 1 (+ i 1))
	   ((> i n) nil)
	   (tagbody
	    (if (= i 1) (go label20))
	    (setf s1
		  (fpplus (fpabs (aref d (- i 1)))
			  (fpabs (aref d (- i 2)))))
	    (setf s2
		  (fpplus s1 (fpabs (aref e (1- i)))))
	    (if (fp> s2 s1) (go label40))
	    label20
	     (setf (aref e2 (- i 1)) '(0 0))
	    label40))
          
     (setf isturm 1)
     (go label320)
     label60 (setf m s)
     (setf x1 lb)
     (setf isturm 2)
     (go label320)
     label80 
     (setf m (- m s))
     (if (> m mm) (go label980))
     (setf q 0)
     (setf r 0)
     label100 (if (= r m) (go label1001))
     (setf tag (+ tag 1))
     (setf p (+ q 1))
     (setf xu (aref d (1- p)))
     (setf x0 xu)
     (setf u (tomp 0))

     (fdo  (q p (+ q 1))
	   ((> q n) nil)
	   (tagbody
	    (setf x1 u)
	    (setf u (tomp 0))
	    (setf v (tomp 0))
	    (if (= q n) (go label110))
	    (setf u  (fpabs (aref e q)))
	    (setf v(aref e2 q))
	    label110
	    (setf xu
		  (fpmin (fpdifference (aref d (1- q))
				      (fpplus x1 u))
			     xu))
	    (setf x0
		  (fpmax (fpplus (aref d (1- q))
				      (fpplus x1 u))
			     x0))
	    (if (fpzerop v) (go label140))))
	       
     label140 
     (setf x1 (fptimes* (fpmax (fpabs xu) (fpabs x0)) machep))
     (if (fpminusp eps1) (setf eps1 (fpminus x1)))
     (if (not(equal p q)) (go label180))
     (if (or (fp> t1 (aref d (1- p)))
	     (fp>= (aref d (1- p)) t2))
	 (go label940))
     (setf m1 p)
     (setf m2 p)
     (setf (aref rv5 (1- p))(aref d (1- p)))
     (go label900)
     label180 
     (setf x1
	   (fptimes* x1  (tomp(+ (-  q p) 1))))
     (setf lb (fpmax t1 (fpdifference xu x1)))
     (setf ub (fpmin t2 (fpplus x0 x1)))
     (setf x1 lb)
     (setf isturm 3)
     (go label320)
     label200 (setf m1 (+ s 1))
     (setf x1 ub)
     (setf isturm 4)
     (go label320)
     label220 
;      (format t "~% label 220 m1=~s, m2=~s, s=~s" m1 m2 s)
      (setf m2 s)
     (if (> m1 m2) (go label940))
     (setf x0 ub)
     (setf isturm 5)
     (fdo  (i m1 (+ i 1))
	   ((> i m2) nil)
	   (tagbody
	    (setf (aref rv5 (1- i))  ub)
	    (setf (aref rv4 (1- i))  lb)))
					; label240
     (setf k m2)
     label250 
;      (print 250)
      (setf xu lb)
     (fdo  (ii m1 (+ ii 1))
	   ((> ii k) nil)
	   (tagbody
	    (setf i
		  (- (+ m1 k) ii))
	    (if (not (fpgreaterp xu (aref rv4 (1- i))))
		(go label260))
	    (setf xu
		  (aref rv4 (1- i)))
	    (go label280)
	    label260))
     label280 
;      (print 280)
      (if (fp> x0
			   (aref rv5 (1- k)))
		  (setf x0 (aref rv5 (1- k))))
     label300 
;      (print 300)
;      (format t "~% xu=~x x0=~s half=~s two=~s eps1=~s s1=~s s2=~s" xu x0 half two eps1 s1 s2)
      (setf x1 (fptimes* (fpplus xu x0) half))
      
     (setf s1 (fptimes* two (fpplus (fpabs xu) (fpplus(fpabs x0) (fpabs eps1)))))
     (setf s2 (fpplus s1 (fpabs (fpdifference x0 xu))))
     (if (equal s2 s1) (go label420))
     label320 
;      (print 320)
      (setf s (- p 1))
     (setf u (tomp 1))
     (fdo  (i p (+ i 1))
	   ((> i q) nil)
	   (tagbody
	    (if (not(equal (car u) 0)) (go label325))
	    (setf v
		  (fpquotient (fpabs (aref e (1- i)))
			   machep))
	    (if (equal (car(aref e2 (1- i))) 0)
		(setf v (tomp 0)))
	    (go label330)
	    label325

	    (setf v
		  (fpquotient (aref e2 (1- i))
			   u))
	    label330
	    (setf u
		  (fpdifference (aref d (1- i))
			   (fpplus x1 v)))
	    (if (minusp (car u)) (setf s (+ s 1))) ))
					;label340
;      (format t "~%isturm =~s" isturm)
      (case isturm 
	(1 (go label60))(2 (go label80))
	(3 (go label200))(4 (go label220))(5 (go label360)))
     label360 
;      (print 360)
      (if (>= s k)  (go label400))
     (setf xu x1)
     (if (>= s m1) (go label380))
     (setf (aref rv4 (1- m1))	   x1)
     (go label300)

     label380 
;      (print 380)
      (setf (aref rv4 s)   x1)
     (if (fpgreaterp (aref rv5 (1- s))  x1)
	 (setf (aref rv5 (1- s))    x1))
     (go label300)
     label400 
;      (print 400)
      (setf x0 x1)
     (go label300)
     label420 
;      (print 420)
      (setf (aref rv5 (1- k))  x1)
     (decf k)
     (if (>= k m1) (go label250))
     label900 
;      (print 900)
      (setf s r)
     (setf r	 (+ (- (+ r m2) m1) 1))
     (setf j 1)
     (setf k m1)
     (fdo  (l 1 (+ l 1))
	   ((> l r) nil)
	   (tagbody
	    (if (> j s) (go label910))
	    (if (> k m2) (go label940))

	    (if (not(fplessp(aref rv5 (1- k))(aref w (1- l))))
		(go label915))
	    (fdo  (ii j (+ ii 1))
		  ((> ii s) nil)
		  (tagbody
		   (setf i (-  (+ l s) ii))
		   (setf (aref w i)(aref w (1- i)))
		   (setf (aref ind i)(aref ind (1- i)))))
					;label905
	    label910
;	     (print 910)
	    (setf (aref w (1- l))(aref rv5 (1- k)))
	    (setf (aref ind (1- l)) tag)
	    (incf k)
	    (go label920)
	    label915
;	     (print 915)
	    (incf j)
	    label920))
     label940 
;      (print 940)
      (if (< q n) (go label100))
     (go label1001)
     label980 (setf ierr (+ (* 3 n) 1))
     label1001 (setf lb t1)
      (setf ub t2)
      (if (= 0 ierr) nil (format t "~%error: eigenvalues by bisect, ierr=~s" ierr))
     (return (values w ierr)) ;; or check ierr if you want error info
     )))


;;
#| here's a 5X5  matrix, tridiagonal, symmetric.
{{0.0625, 0.0161374, 0., 0., 0.}, 
{0.0161374, 0.0208333, 0.00704295, 0., 0.},
{0., 0.00704295, 0.0104167, 0.00393713, 0.},
{0., 0.,   0.00393713, 0.00625, 0.00251259}, 
{0., 0., 0., 0.00251259,   0.00416667}}
whose eigenvalues are
{0.0681107, 0.0203175, 0.00960646, 0.00478623, 0.00134577}
|#

;; set up calling program for bigfloat eigenvalue computation
;; of symmetric tridiagonal matrix with diagonal d, off-diagonal of,
;; look only for number (default 1)
;; eigenvalues between lb lower bound and ub upper bound
;; to accuracy eps.

;; 
;; d is a maxima list, of is a maxima list
(defun $bfbisect(d of lb ub &optional (eps 1.0e-8) (number 1))
  (let* ((n (length (cdr d)))
	 (diag (make-array n   :initial-contents
			  (map 'list #'(lambda(x)(tomp x ))
			       (cdr d)))) ;; take off (mlist)

	(off (make-array n 
			  :initial-contents
			  (map 'list #'(lambda(x)(tomp x ))
			       (cdr of))))
	 (sqs (make-array n
			  :initial-contents
			  (map 'list #'(lambda(x)(fptimes* x x))off)))
	 (w (make-array n :initial-element (tomp 0)))
	 (rv4 (make-array n :initial-element (tomp 0) ))
	 (rv5 (make-array n :initial-element (tomp 0)))
	 (ind (make-array n :element-type 'fixnum
			  :initial-element 0)))


	 #+ignore (loop for i from 0 to (1- n) do  ;; no need to initialize unique elements
	       (setf (aref   w i) (tomp 0))
	       (setf (aref rv4 i) (tomp 0))
	       (setf (aref rv5 i) (tomp 0)))
	       
	 (bisect  n (tomp eps) diag off sqs 
		 (tomp lb )
		 (tomp ub )
		 number
		 0 w ind rv4 rv5)
	(cons '(mlist) (mapcar #'bcons (coerce w 'list)) )))

(defun test2()
 (setf $diag  '((mlist)0.0625 0.0208333 0.0104167 0.00625 0.00416667))
 (setf $offdiag '((mlist) 0.0 0.0161374 0.00704295 0.00393713 0.00251259))
  ($bfbisect $diag $offdiag 0 1 1.0d-17 5))

;; after executing this ... w contains something like
;((mlist) 0.0013457789 0.0047862385 0.009606467 0.02031746 0.06811072)