;; -*- Mode: Lisp; Syntax: Common-Lisp ;package: user; -*-

#| Lisp routines for finding polynomial zeros based on
W. Kahan's "Notes on Laguerre's Iteration"  4 Dec. 1992


A polynomial p of degree n-1 is represented by a lisp simple-vector of
length n, a[0..n-1]. The polynomial p(x) is then

      n-1     n-i
      sum a  x
      i=0  i 

The polynomial zero is a vector of length 0, #().
Constants look like vectors of length 1. For example, 2 = #(2.0d0).
The vector #(1 2) represents x+2.
You can create double-float polynomials conveniently by, for example,
(setf p (make-dfa 3 1/2 1/3 2.718 0.577215d0)) ==>
#(3.0d0 0.5d0 0.3333333333333333d0 2.7179999351501465d0 0.577215d0)

note conversions from integer, rational single and double.

|#

;; Assume IEEE 754 arithmetic for now.

(eval-when (compile eval load)
  (defvar double-eps (expt 2.0d0 -53)))
(declaim (double-float double-eps))

#| Implementation notes.. When we divide by 0.0d0 we
expect to get #.EXCL::*INFINITY-DOUBLE*
or maybe #.EXCL::*NEGATIVE-INFINITY-DOUBLE* .
We do not want a lisp error "divide by zero" which
is likely to happen if these programs are run with
error checking, or interpreted.
|#


;;;EVALUATION

(defun prreval (a x)
"evaluate a polynomial with real double precision coefficients
 at a real argument: return p, e, q, r"
;; p is approx p(x) = sum(a[k]*x^k)
;; e is a bound on the error in evaluating p:  e>= |p-p(x)|
;; q is approx p'(x)/p(x)
;; r is approx q^2-p''(x)/p(x)

  (declare (double-float x)
	   (type (simple-array double-float 1) a)
	   (optimize speed (safety 0) (debug 0)))

  (let* ((p (aref a 0)) (q 0.0d0) (r 0.0d0) (mu (abs x)) 
			(e (* 0.5d0 (abs p))) (n (1- (length a))))
    (declare	       ;;(:explain :calls :types)
		       (double-float p q r mu e)
		       (fixnum n))

      (do ((j 1 (1+ j)))
	  ((> j n) 
	   (values	;; return p e q r
	    p
	    (* #.double-eps (+ (- e (abs p)) e))
	    (/ q p)
	    (- (* q q) (* 2.0d0 (/ r p)))))
	(declare (fixnum j))
	(setf r (+ (* x r) q))
	(setf q (+ (* x q) p))
	(setf p (+ (* x p) (aref a j)))
	(setf e (+ (* mu e) (abs p)))
	)))



(defun prreval2 (a x)
"evaluate a polynomial with real double precision coefficients
 at a real argument: return p, error"

  (declare (double-float x)
	   (type (simple-array double-float 1) a)
	   (optimize speed (safety 0) (debug 0)))

  (let* ((p (aref a 0)) (mu (abs x)) 
			(e (* 0.5d0 (abs p))) (n (1- (length a))))
    (declare	       ;;(:explain :calls :types)
		       (double-float p mu e)
		       (fixnum n))

      (do ((j 1 (1+ j)))
	  ((> j n) 
	   (values	;; return p e
	    p
	    (* #.double-eps (+ (- e (abs p)) e))))
	(declare (fixnum j))
	(setf p (+ (* x p) (aref a j)))
	(setf e (+ (* mu e) (abs p)))
	)))

(defun prreval3 (a x)
"evaluate a polynomial with real double precision coefficients
 at a real argument: return p only"
;; this function is almost all declarations..
  (declare (double-float x)
	   (type (simple-array double-float 1) a)
	   (optimize speed (safety 0) (debug 0)))

  (let* ((p (aref a 0)) (n (1- (length a))))
    (declare (double-float p) (fixnum n))

      (do ((j 1 (1+ j)))
	  ((> j n) p)
	(declare (fixnum j))
	(setf p (+ (* x p) (aref a j))))))

;; from knuth vol 2 4.6.4 2nd ed p 468
;;  "A good programmer does not make indiscrimate use of the
;; built-in "complex arithmetic" features of high-level
;; programming languages."

(defun prceval3 (u z)
  "evaluate a polynomial with real double precision coefficients
 at a COMPLEX argument: return (list a b) only"
  (declare (type (simple-array double-float 1) u)
	   (type z (complex double-float)))

  (let* ((x  (realpart z))
	 (y  (imagpart z))
	 (r  (+ x x))
	 (s (+ (* x x)(* y y)))
	 )
    ;; this function is almost all declarations..
    (declare (double-float x y r s a b a0 b0) )
    (let*   
	((a0 (aref u 0))
	 (b0 (aref u 1))
	 (a 0.0d0)
	 (b 0.0d0)
	 (n (1- (length u))))
      (declare (double-float a0 b0 a b) (fixnum n))
      (if (= n 1)
	  (+ b0 (* z a0))
      (do ((j 2 (1+ j)))
	  ((> j n) (+ (* z a) b))		;return value
	(declare (fixnum j))
	(setf a (+ b0 (* r a0)))
	(setf b (- (aref u j) (* s a0)))
	(setf a0 a b0 b))))))


(defmacro framevars(z);;utility to generate vars etc
  (newsym "W" -1)
  (newsymvarsclear)
  `(prog1
       (lda-let 
	(let ((compute ,z))
	  `(let* 
	       ,(mapcar #'(lambda(x)(list (car x)(simp(cadr x))))
			(newsymvars))
	     ,(simp compute))))
     ;; report on usage of arithmetic in temp variables
     (format t ";; ~s adds, ~s mults, ~s divides" 
	     (newsymplus)(newsymtimes)(newsymdivides))))


;; symbolically do expansion of u(z), u is sequence of symbols, z is
;; complex number

(defun prc3-s(u z)(framevars (prceval3-s u z)))

(defun prceval3-s (u z)
  "evaluate a polynomial with real double precision coefficients
 at a COMPLEX argument: return (list a b) only"

 
  (let* ((x (newname (simp `(realpart ,z))))
	 (y (newname (simp `(imagpart ,z))))
	 (r  (newname(+s x x)))
	 (s (newname(+s (*s x x)(*s y y))))
	 )
    ;; this function is almost all declarations..
    (let*   
	((a0 (elt u 0)) ;was aref
	 (b0 (elt u 1))
	 (a 0.0d0)
	 (b 0.0d0)
	 (n (1- (length u))))
      (if (= n 1)
	  (+s b0 (*s z a0))
      (do ((j 2 (1+ j)))
	  ((> j n) (+s (*s z a) b))		;return value
	(declare (fixnum j))
	(setf a (newname(+s b0 (*s r a0))))
	(setf b (newname(-s (elt u j) (*s s a0)))) ;was aref
	(setf a0 a b0 b))))))

;;; symbolically

(defun pn-s(u z) (framevars (pneval-s u z)))

(defun pneval-s (a x)
  "SYMBOLICALLY evaluate a polynomial with any numeric coefficients at
  at any numeric argument: return p, e, q, r"
  (let* ((p (elt a 0)) (q 0) (r 0)
	 (mu (newname(abs-s x)) )
	 (e (newname (*s 0.7d0 (abs-s p))))
	 (n (1- (length a)))
	 (compute 
	  (do ((j 1 (1+ j)))
	      ((> j n) 
	       (list 'values
		     p
		     (*s 'eps (+s (*s 2.31d0 (-s e (abs-s p))) e))
		     (/s q p)
		     (-s (*s q q) (*s 2 (/s r p)))))
	    (setf r (newname(+s (*s x r) q)))
	    (setf q (newname(+s (*s x q) p)))
	    (setf p (newname(+s (*s x p) (elt a j))))
	    (setf e (newname(+s (*s mu e) (abs-s p)))))
	  ))
    compute))

(defun abs-s(x)(simp (list 'abs x)))
(defun -s(&rest x)(cons '- x))

#| evaluate a polynomial with complex double precision coefficients at
a complex argument. This is a lot slower in our lisp, since complex
arithmetic and abs is not compiled in-line, and extra data is
allocated for complex quantities. This could be hacked to use such
tricks. |#

;; the program below just uses complex arith.
;; one could evaluate p(z) as preal(z)+i*pimage(z) using
;; other programs here. (Ask WK why/not?)
(defun pcceval (a x)
"evaluate a polynomial with complex double precision coefficients
 at a complex argument: return p, e, q, r"

  (declare (type (complex double-float) x)
	   (type (simple-array (complex double-float) 1) a)
	   (optimize speed (safety 0) (debug 0)))

  (let* ((p (aref a 0)) (q 0.0d0) (r 0.0d0) (mu (abs x)) 
			(e (* 0.7d0 (abs p))) (n (1- (length a))))
    (declare	       ;;(:explain :calls :types)
     (type (complex double-float) p q r)
     (double-float mu e)
     (fixnum n))

      (do ((j 1 (1+ j)))
	  ((> j n) 
	   (values	;; return p e q r
	    p
	    (* #.double-eps (+ (* 2.31d0 (- e (abs p))) e))
	    (/ q p)
	    (- (* q q) (* 2.0d0 (/ r p)))))
	(declare (fixnum j))
	(setf r (+ (* x r) q))
	(setf q (+ (* x q) p))
	(setf p (+ (* x p) (aref a j)))
	(setf e (+ (* mu e) (abs p))))))


(defun pneval (a x)
"evaluate a polynomial with any numeric coefficients at
 at any numeric argument: return p, e, q, r"

  (declare 
	   (type (simple-array t 1) a)
	   (optimize speed (safety 0) (debug 0)))

  (let* ((p (aref a 0)) (q 0.0d0) (r 0.0d0) (mu (abs x)) 
			(e (* 0.7d0 (abs p))) (n (1- (length a))))
    (declare	       ;;(:explain :calls :types)
     ;(type (complex double-float) p q r)
     ;(double-float mu e)
     (fixnum n))

      (do ((j 1 (1+ j)))
	  ((> j n) 
	   (values	;; return p e q r
	    p
	    (* #.double-eps (+ (* 2.31d0 (- e (abs p))) e))
	    (/ q p)
	    (- (* q q) (* 2.0d0 (/ r p)))))
	(declare (fixnum j))
	(setf r (+ (* x r) q))
	(setf q (+ (* x q) p))
	(setf p (+ (* x p) (aref a j)))
	(setf e (+ (* mu e) (abs p))))))

;; generalize pneval to handle complex bigfloats as defined in
;; cbf.lisp in mma.src


;; The program below is rather wasteful, in that it
;; is repeatedly consing-up new data even when the information
;; is pretty much guaranteed to be temporary.
;; A more careful program would use a kind of register organization
;; where data could be accumulated by over-writing old values
;; rather than creating new ones.
;; it's probably a good idea
;; to make sure that 
;; bigfloat-init-bin has been run at least once
;; before doing anything

(defun pbigneval (olda oldx precision)
"evaluate a polynomial with bigfloat complex coefficients at
 at any numeric or bigfloat or complex bigfloat argument: return p, e, q, r"

  (declare 
	   (type (simple-array t 1) olda a)
	   (optimize speed (safety 0) (debug 0)))

  (let* ((bigfloat-bin-prec precision)
	 (bigfloatone (intofp 1));;initialize one
	 (bigfloatzero (bcons2 0 0));; initialize zero
	 (a (make-array t (length olda)))
	 (x (cbigfloat-convert oldx))
	 (p (aref a 0)) (q bigfloatzero) (r bigfloatzero) 
	 (mu (cbigfloat-abs x)) 
	 (e (bigfloat-* (bigfloat-convert 7/10)
			(cbigfloat-abs p)))
	 (n (1- (length a))))
    (declare;;(:explain :calls :types)
     (special bigfloat-bin-prec bigfloatone bigfloatzero)
     (fixnum n))
    (dotimes (i (length a))
      (setf (aref a i)(cbigfloat-convert (aref olda i))))

    (do ((j 1 (1+ j)))
	((> j n) 
	 (values;; return p e q r
	  p
	  (bigfloat-* (bigfloat-expt (intofp 2) 
				     (- bigfloat-bin-prec))
		      (bigfloat-+ (bigfloat-* (bigfloat-convert 231/100)
					      (bigfloat-- e (cbigfloat-abs p))) e))
	  (cbigfloat-/ q p)
	  (cbigfloat-- (cbigfloat-* q q) (cbigfloat-* (intofp 2) (cbigfloat-/ r p)))))
      (declare (fixnum j))
      (setf r (cbigfloat-+ (cbigfloat-* x r) q))
      (setf q (cbigfloat-+ (cbigfloat-* x q) p))
      (setf p (cbigfloat-+ (cbigfloat-* x p) (aref a j)))
      (setf e (cbigfloat-+ (cbigfloat-* mu e) (cbigfloat-abs p))))))

(defun pneval2 (a x)
"evaluate a polynomial with any numeric coefficients at
 at any numeric argument: return p, e"

  (declare 
	   (type (simple-array t 1) a)
	   (optimize speed (safety 0) (debug 0)))

  (let* ((p (aref a 0)) (q 0.0d0) (r 0.0d0) (mu (abs x)) 
			(e (* 0.7d0 (abs p))) (n (1- (length a))))
    (declare	       ;;(:explain :calls :types)
     ;(type (complex double-float) p q r)
     ;(double-float mu e)
     (fixnum n))

      (do ((j 1 (1+ j)))
	  ((> j n) 
	   (values	;; return p e q r
	    p
	    (* #.double-eps (+ (* 2.31d0 (- e (abs p))) e))
	    ;(/ q p)
	    ;(- (* q q) (* 2.0d0 (/ r p)))
	    ))
	(declare (fixnum j))
	;;(setf r (+ (* x r) q))
	;;(setf q (+ (* x q) p))
	(setf p (+ (* x p) (aref a j)))
	(setf e (+ (* mu e) (abs p))))))

(defun pnneval3 (a x)
"evaluate a polynomial of any types at any argument: return p only"
  (declare  (optimize speed (safety 0) (debug 0)))
  (let* ((p (aref a 0)) (n (1- (length a))))
    (declare(fixnum n))
      (do ((j 1 (1+ j)))
	  ((> j n) p)
	(declare (fixnum j))
	(setf p (+ (* x p) (aref a j))))))

(defun pnneval3-s (a x)
  "symbolically evaluate a polynomial at any argument: return p only"
  (newsym "W" -1)
  (newsymvarsclear)
  (let ((compute 
	 (let* ((p (elt a 0)) (n (1- (length a))))
	   (do ((j 1 (1+ j)))
	       ((> j n) p)
	     (setf p (newname(+s (*s x p) (elt a j))))))))
    `(let* ,(newsymvars) , compute)))

(defun evens(a)
  "every other element of a list: 0, 2 .."
  (if (null a) nil (cons (car a)(evens (cddr a)))))
       
(defun pnneval4-s (a x)
  "symbolically evaluate a polynomial at any argument: return p only,
   generate 2 non-interfering sequences."
  ;; example (pnneval4-s '(a b c d e f g) 'x)
  (framevars
   (let* ((evenp (evenp (length a)))
	  (ee (if evenp (evens (cdr a))(evens a)))
	  (oo (if evenp (evens a)(evens (cdr a))))
	  (x2 (newname (*s x x)))
	  (p (car ee))
	  (p2 (car oo)))
     (pop ee)(pop oo)
     (loop while (or ee oo) do
	   (if ee (setf p (newname(+s (*s x2 p) (pop ee)))))
	   (if oo (setf p2 (newname(+s (*s x2 p2) (pop oo))))));end loop
     (+s p (*s x p2)))))




(defun pevalgeneric (a x plus times)
  "evaluate a polynomial of any types at any argument: return p only.
   use generic arithmetic plus and times, which could be bigfloat,
   complex bigfloat, interval, etc."
  (declare  (optimize speed (safety 0) (debug 0)))
  (let* ((p (aref a 0)) (n (1- (length a))))
    (declare  (fixnum n))

      (do ((j 1 (1+ j)))
	  ((> j n) p)
	(declare (fixnum j))
	(setf p (funcall plus (funcall times  x p) (aref a j))))))

;;; DEFLATION

(defun deflaterr (a z)
"deflate a real polynomial a by the real root z"
  (declare (double-float z)
	   (type (simple-array double-float 1) a)
	   (optimize speed (safety 0) (debug 0)))

  (let* ((n (1- (length a)))
	 (b (make-array n :element-type 'double-float)))

    (declare (type (simple-array double-float 1) b)  (fixnum n))
    
    (setf (aref b 0) (aref a 0))
    (do ((j 1 (1+ j)))
	((= j n) b)
      
      (declare (fixnum j))
       
    (setf (aref b j) (+ (* z (aref b (1- j))) (aref a j)))
      )))

(defun deflatenn (a z)
"deflate a  polynomial a by the root z <any numeric types>"
  (declare ;(double-float z)
	  ; (type (simple-array double-float 1) a)
	   (optimize speed (safety 0) (debug 0)))

  (let* ((n (1- (length a)))
	 ;;(b (make-array n :element-type 'double-float))
	 (b (make-array n)))

    (declare ;;(type (simple-array double-float 1) b)
     (fixnum n))
    
    (setf (aref b 0) (aref a 0))
    (do ((j 1 (1+ j)))
	((= j n) b)
      
      (declare (fixnum j))
       
    (setf (aref b j) (+ (* z (aref b (1- j))) (aref a j)))
    )))
  
;;; INITIALIZATION

#| notes .. to create a simple-array of double-floats, it is not
enough to just type in (setf a #(1.0d0 2.0d0 3.0d0))
but it is necessary to do something like this
(coerce a '(simple-array double-float *))
or to define such an array initially  by, for example
(make-dfa  1 2 3) defined below.
 this creates the encoding for the polynomial x^2+2*x+3, by the way.

Note that Common Lisp distinguishes between an array of pointers to
floats and an array of floats.  The latter is more efficient, and is
required by the programs here. One is created by make-dfa.

|#

#+ignore
(defmacro make-dfa(&rest l)
  "make an array of double-floats"
  `(make-array ,(length l) :element-type 'double-float :initial-contents 
	       ',(mapcar #'(lambda(n)(coerce n 'double-float)) l)))

(defun make-dfa(&rest l)
  "make an array of double-floats"
  (make-array (length l) :element-type 'double-float :initial-contents 
	      (mapcar #'(lambda(n)(coerce n 'double-float)) l)))

(defun make-cdfa(&rest l)
  "make an array of double-floats"
  (make-array (length l) :element-type '(complex double-float) :initial-contents 
	       (mapcar #'(lambda(n)(coerce n '(complex double-float))) l)))

#+ignore
(defmacro make-cdfa(&rest l)
  `(make-array ,(length l) :element-type '(complex double-float) 
	       :initial-contents 
	       ',(mapcar #'(lambda(n)(coerce n '(complex double-float))) l)))


;; the next calc is sloppy-- should make sure divisor is non-zero..
;; find a small one..

(defun ahatbound (a)
  #| should compute
  min ( abs( choose(n,j)*a[n]/a[n-j])^(1/j)) for j>0 
  |#
  (let ((n (length a))
	(an (aref a (1- n))))
	(abs(* n (/ an (aref a (- n 2)))))))

(defun ahat(a)
  "return an estimate for smallest zero z
   ahat < |z| < ahat/(2^(1/n)-1)"
  (let*((n (length a))
	(an (aref a (1- n)))
	(capp (subseq a 0 (1- n)))

	(ahatbd (ahatbound a)))
	(expg (log ahatbd)))
    ;;capp= capital P in notes
    (dotimes (i (1- n))
      (setf (aref capp i)(/ (aref capp i) an)))
    ;;   (format t "~%capp =~s ahatbd= ~s" capp ahatbd)
    
  ;; Newton iteration of ln(h(exp(g)))=0 where g=ln(ahat)
 (dotimes (i 4) ;; just a guess that 4 iterations would be ok.
   (multiple-value-bind  (p e q r) 
       (prreval capp expg)			; q is p'/p
     ;;     (format t "~% expg = ~s p=~s q=~s " expg p q)
     (setf expg (- expg (/ (log (abs p)) q)))))

 (values g (/ g (- (expt 2.0 (/ 1 (1+ (length p)))  ) 1))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

#| Some other stuff. 
These routines provide a simple way of making up polynomials with
given roots by doing arithmetic on polynomials themselves.
These could be made more efficient by providing workspace for 
the polynomials rather than having lisp make up new
vectors each time.
|#

(defun addrr(r s) 
  "add 2 real double-float polynomials"  
  (let* ((lr (length r))
	 (ls (length s)))
    (declare (fixnum lr ls)
	     (type (simple-array double-float 1) r s))
    (if (> lr ls) (addrr s r)
      ;; now s is at least as long as r
      (let((res (copy-seq s))  (del (- ls lr)))
	(do ((j del (1+ j)))
	    ((= j ls) (normalize res))
	  
	  (declare (fixnum j del) (type (simple-array double-float 1) res))
	  (setf (aref res j)(+ (aref res j) (aref r (- j del)))))))))


(defun normalize (x)
  ;; remove leading zero coefficients
  (let ((pos (position-if-not #'zerop x)))
    (declare (fixnum pos))
    (if (zerop pos) x) x  (subseq x pos)))

(defun timesfr(f s &optional (offset 0))
  "multiply  a float f*x^(offset) by s, a real double-float polynomial"  
  
  (let ((res (make-array (+ (length s) offset)
			 :element-type 'double-float
			 :initial-element 0.0d0)))
    (declare (type (simple-array double-float 1) res s)
	     (double-float f))
    (if (zerop f) #()
      (dotimes (i (length s) (normalize res)) ; in case of underflow to 0.0
	(declare (fixnum i))
	(setf (aref res i)(* f (aref s i)))))))

(defun timesrr(r s)
  "multiply 2 real double-float polynomials"  
  (let ((res #())
	(n (1- (length r))))
  (dotimes (i (length r) res)
    (setf res (addrr res (timesfr (elt r i) s (- n i)))))))


(defun polywrz(&rest l)
  "double-float polynomial with prescribed zeros"
  (let((res (make-dfa 1)))
    (dolist (i l res)
      (setf res(timesrr 
		(coerce (vector 1.0d0 (* -1.0d0 i))
			'(simple-array double-float *))
		res)))))


(defun newt(h start count)
 " run newton iteration on polynomial h from start, count times."
  (dotimes (i count start)
    (format t "~%x[~a] = ~a" i start)
    (multiple-value-bind (p e q r)
	(prreval h start)
      ;; need to figure how to stop this iteration somehow
      (setq start (- start (/ 1 q))))))


(defun newt(h start count)
 " run newton iteration on polynomial h from start, count times."
  (dotimes (i count start)
    (format t "~%x[~a] = ~a" i start)
    (multiple-value-bind (p e q r)
	(prreval h start)
      (when (<= (abs p) e) 
	;;(format t "~%Can't do better than ~s" p) ;6/00/RJF
	    (return start))
	(setq start (- start (/ 1 q))))))


(defun weier(z p)
  "Weierstrass iteration: update z in place,
   a vector of complex zeros of poly p"
  ;; It works, any lisp numbers including in particular, complex 
  ;; 1/26/99 RJF
  (let ((deg (length z)))
    (declare (fixnum deg)(optimize speed (safety 0)))
    (dotimes(i deg z)			;return z
      (declare (fixnum i))
      (let ((zi (aref z i)))
	(setf (aref z i)
	      (- zi
		 (/ (pnneval3 p zi) ;evaluate any arg, any coeff
		    (let ((prod 1))
		      (dotimes (j deg prod)
			(declare (fixnum j))
			(if (= i j) nil 
			  (setf prod (* prod (- zi (aref z j))))))))))))))

;; here one could pass in #'+ #'* #'- #'/ for lisp numbers  or
;; bigfloat-only equivalents.  or
;; complex-bigfloat-only equivalents.
;; or generic (e.g. promote when necessary).
;; A good strategy might be to destructively upgrade the coefs in polynomial
;; p to higher precision bigfloats.

(defun weiergeneric(z p plus times difference quotient)
  "Weierstrass iteration: update z in place,
   a vector of complex zeros of poly p"
  (let ((deg (length z)))
    (declare (fixnum deg)(optimize speed (safety 0)))
    (dotimes(i deg z)			;return z
      (declare (fixnum i))
      (let ((zi (aref z i)))
	(setf (aref z i)
	      (funcall difference zi
		 (funcall quotient
			  (pevalgeneric p zi plus times) ;evaluate any arg, any coeff
		    (let ((prod 1))
		      (dotimes (j deg prod)
			(declare (fixnum j))
			(if (= i j) nil 
			  (setf prod (funcall times prod (funcall difference zi (aref z j))))))))))))))


#| roots of x^5+x-1 = #(1.0d0 0.0d0 0.0d0 0.0d0 1.0d0 -1.0d0)

#(
 #c(0.50000101661976d0 0.86602373046046d0)
 #c(0.50000101661976d0 -0.86602373046046d0)
#c(-0.8774382922239d0  0.74486261511982d0)
#c(-0.8774382922239d0  -0.74486261511982d0)
0.75487455120826d0)

better values after 3 iterations of weier

#(#c(0.5d0 0.8660254037844386d0) 
  #c(0.5d0 -0.8660254037844387d0)
  #c(-0.8774388331233464d0 0.7448617666197442d0)
  #c(-0.8774388331233464d0 -0.7448617666197442d0)
  #c(0.7548776662466927d0 -6.842277657836021d-49))

  |#

(defun acceptroots(z p)
  "if a root in sequence z satisfies p(z)<2*(error in evaluating p(z)) then
   accept it."	
  ;;Returns 3 values: list of good roots, list of bad roots, and
  ;;the deflated polynomial newp.
  (let ((goodroots nil)
	(badroots nil)
	(newp p)
	)
    (map nil
      #'(lambda(zero)(multiple-value-bind (val e)
		      (pneval2 p zero)
		    (cond((<= (abs val) (* 2 e) )
			  (push zero goodroots)
			  (setf newp (deflatenn newp zero)))
			 (t (push zero badroots)))))
      z)
    (values goodroots badroots newp)))



;;;; simple lisp hacks.  A numerical  polynomial is
;;;; represented by a list of its numerical coefficients,
;;;; highest to lowest, e.g. (a b c d) means a*x^3+b*x^2+c*x+d
;;;;

;;evaluate the polynomial


(defun ev(p x) (reduce #'(lambda( s r)(+ r (* x s))) p))

(defun ev-s(p x)
  (framevars(reduce #'(lambda( s r)(newname(+s r (*s x s)))) p)))

(defun ev-mt(p x)
  (framevars (reduce #'(lambda(s r)(mtrename(+m r (*m x s)))) p)))

(defun mtrename (x)
  (list 'mt (newname x) (mtlength x)))


(defun evabs(p x) (reduce #'(lambda(s r)(+ (abs r) (* x s))) p))

;;return the derivative of the polynomial

(defun df(p) (df1 p (1- (length p))))

(defun df1(p r)(cond ((null(cdr p))nil)
		     (t (cons (* r (car p))(df1 (cdr p) (1- r))))))
	       
;;;;;;;;;;chain of recurrences stuff  Bachmann et al ISSAC 94
;;; I didn't do too much of this; O. Bachmann should have
;;; this coded...

(defun accum(op init fun seq)
	     (if (null seq) init
		 (funcall op (funcall fun (car seq))
			  (accum op init fun (cdr seq)))))

;; given 3 arguments  for basic recurrence
;; produces a function.
;; e.g.
;;(br 1 #'* #'(lambda(r) (1+ r))) --> 


(defun br(phi op f1)
  `(lambda(i)(accumn ,op (quote ,phi) ,f1 0 (1- i))))

(defun accumn(op init fun low hi)
	     (if (> low hi) init
		 (funcall-s op (funcall-s fun low)
			  (accumn op init fun (1+ low) hi))))


(defun funcall-s(h &rest args)
  (if (or
       (functionp h)
       (and (listp h)(eq (car h) 'lambda)))
      (apply h args) ;; we can't apply. just simplify args
	(cons h (mapcar #'simp args))))


;; Bachmann et al ISSAC 94.

(defun creval (cr n)
  ;; res is an array 0..n-1
  ;; cr is an array of 2*k+1  of alternating ops and values
  (let* ((lim (-(length cr) 2))
	 (res (make-array n)))
  (loop for i to (1- n) do
	(setf (elt res i) (elt cr 0))
	(loop for j from 0 to lim by 2 do
	     (setf (elt cr j)(funcall (elt cr (1+ j)); an operator + or *
				      (elt cr j)
				      (elt cr (+ j 2))))))
  res))

(defun creval-s (cr n)
  ;; generate a program
;; res is an array 0..n-1
;; cr is an array of 2*k+1  of alternating ops and values
  (let* ((lim (-(length cr) 2))
	 (res (make-array n)))
    (framevars
     (let()
       (loop for i to (- n 1) do
	(setf (elt res i) (newname(elt cr 0)))
	(loop for j from 0 to (min lim  (* 2 (- n i 2))) by 2 do
	      (setf (elt cr j)
		(newname
		 (if (eq(elt cr (1+ j)) '+)
		       (+s  (elt cr j) (elt cr (+ j 2)))
		       (*s  (elt cr j) (elt cr (+ j 2))))))))

       (cons 'vector (coerce res 'list))))))

;;(creval '(0 + 1 + 6 + 6) 5) returns #(0 1 8 27 64)
;; because that CR is the encoding for x^3+0x^2+0*x+0.

  
  

;; exercise 7 section 4.6.4 in Knuth ;; need to look at 3rd ed of vol 2 for right ans

(defun tab1 (pol x0 h); pol is coefficients in polynomial
  (let*((deg (1- (length pol)))		;length = 1+degree of polynomial
	(y (apply #'vector 
		  (loop for i to deg collect (ev pol (+ x0 (* i h)))))))
    (loop for k from 1 to deg do
	  (loop for j from deg downto k  do
		(decf (aref y j) (aref y (1- j)))))
    y))

(defun march(y) ; a0 <- a0+a1, .... a[n-1]<- a[n-1]+a[n]
   (loop for k to (- (length y) 2) do
	 (incf (aref y k) (aref y (1+ k))))
   y)

;; (march (tab1 '(1 2 3) 0 0.5) 10)  ; look at the leftmost element for f(0)
;; (march *)  ;look at leftmost element again for f(0.5) ...
;; (march *)

;;; make a program that takes a polynomial with numeric coefficients
;;; and produces a program that generates many results.

(defun tab2-s (pol x0 h howmany)	; pol is coefficients in polynomial
  (if (not (numberp howmany))
      (tab2-sn pol x0 h howmany)
  (let ((deg (1- (length pol)));;(a b c) is length 3 degree 2
	;;(theprog (compile nil `(lambda (thevar),(ev-s pol 'thevar)))) ;faster?
	(evprog `(function(lambda (thevar),(ev-s pol 'thevar)))))
    (framevars
     (let*((theprog (newname evprog))
	   (y 
	    (cons (newname`(funcall ,theprog , x0 ))
	    (loop for i from 1 to deg collect
		  (prog2
		      (setf x0 (newname (+s x0 h)))
		  (newname`(funcall ,theprog , x0 )))
		    ))))
       (loop for k from 1 to deg do
	     (loop for j from deg downto k do
		   (setf (elt y j) (newname(-s (elt y j)(elt y (1- j)))))))
    
       (cons 'vector
	     (loop for count from 0 to (1- howmany) collect
		       
		   (let ((items (loop for k to (min (- howmany (+ 2 count))
						    (- (length y) 2)) collect
						    `(incf ,(elt y k),(elt y (1+ k))))))
		     (if  items  `(prog1 ,(elt y 0),@items)
		       (elt y 0)))
		   )))))))

(defun tab2-sn (pol x0 h howmany)	; pol is coefficients in polynomial
  (let ((deg (1- (length pol)));;(a b c) is length 3 degree 2
	;;(theprog (compile nil `(lambda (thevar),(ev-s pol 'thevar)))) ;faster?
	(evprog `(function(lambda (thevar),(ev-s pol 'thevar)))))
    (framevars
     (let*((theprog (newname evprog))
	   (m0 x0)
	   (y 
	    (loop for i to deg collect
		  (prog1
		      (newname`(funcall ,theprog , m0 ))
		    (setf m0 (newname (+s m0 h))))))
	   (marchlist `(prog1 ,(elt y 0),@(loop for k to (- (length y) 2) collect
				      `(incf ,(elt y k),(elt y (1+ k)))))))
       (loop for k from 1 to deg do
	     (loop for j from deg downto k do
		   (setf (elt y j) (newname(-s (elt y j)(elt y (1- j)))))))
       `(loop for count from 0 to (1- ,howmany) collect ,marchlist)))))



;;; We need some utility programs for the halving scheme for monic
;;; polys from Knuth, vol 2 3rd ed solution to problem 65.  invented
;;; by S. Winograd.


;; synthetic division of polynomials

(defun sd (a b);; a simple implementation of knuth Alg D 4.6.1
  (let* ((u (nreverse(apply #'vector a))); reverse to get indexes right.
	 (v (nreverse(apply #'vector b)))
	 (m (1-(length u)))
	 (n (1-(length v)))
	 (q (make-array (1+(- m n))))
	)
    (loop for k from (- m n) downto 0 do
	  (let ((qk (/ (elt u (+ n k))(elt v n))))
	    (setf (elt q k) qk)
	    (loop for j from (1- (+ n k)) downto k do
		  (decf (elt u j)(* qk (elt v (- j k)))))))
    (values
     (coerce (nreverse q)'list)			;quot
     (coerce (nreverse (subseq u 0  n))'list)	;rem
     )))

(defun part(pol n) ;; (a b c d e f) 4  --> (a b c d-1)  (1 e f)
  (let* ((left (subseq pol 0 n))
	 (right (subseq pol  n ))
	 (lastl (last left)))
    (rplaca lastl (1- (car lastl)))
    (values left (cons 1 right))))


;; u is a MONIC polynomial (lead coeff 1). Scheme due to S. Winograd

(defun win65(u x)
  (let* ((n (1- (length u)))
	 (l (floor (log n 2)))
	 (xpow (make-array (1+ l))))
    
    (if (not (equalp (car u) 1)) (error "not monic ~s" u))
    (if (not (> n 1)) (error "degree ~s is not interesting here" n))
    
    (framevars
     (let()
       ;; compute some powers of x in advance
	 (setf (aref xpow 0) (newname x));x^(2^0)=x
	 (loop for i from 1 to l do
	       ;;xpow[i]contains the temporary for
	       ;;x^(2^i) = x^2, x^4, ..
	       (setf (aref xpow i)
		 (newname (*s (aref xpow (1- i))
				(aref xpow (1- i))))))
	 ;; generate the rest of the program
	 (w65x u xpow)))))

;; interior program for win65

(defun w65x (u xpow)
  (if (not (equalp (car u) 1)) (error "not monic ~s" u))
  ;; u must be MONIC
  (let* ((n (1- (length u)))		;n is degree
	 (l (floor (log n 2)))) 
    (cond
     ((= n 0)(error "not monic")(car u)); should not happen
     ;; dispense with remaining base cases.
     ((= n 1)
      (newname
       (+s  (aref xpow 0) (cadr u))))	;linear (x+a))
     ((= n 2)
      (newname
       (+s				;quadratic (x^2+a*x+b)
	(aref xpow 1)			; x^2
	(*s (aref xpow 0) (cadr u))	; x*a
	(caddr u)			; b
	)))
     ((= n (expt 2 l))			;n is 2^l
      ;; use one step of horner's rule to reduce degree
      (let ((h (reverse u)))
	(newname (+s (car h)		; last term /constant
		       (*s (aref xpow 0)
			   (w65x (reverse (cdr h)) xpow))))))
   
     ((= n (1- (expt 2 (1+ l))))	;n is 2^(l+1)-1 = 2m+1
      ;; use synthetic division
      (let* ((m (/ (1- n) 2))
	     (alpha (1- (elt u (1+ m))))
	     (divisor (cons 1;; x^(m+1)+alpha
			    (append (make-list  m :initial-element 0)
				    (list alpha))))
	     (v nil)
	     (w nil))
	(multiple-value-setq
	    (v w) (sd u divisor))
	(newname(+s (*s (+s (aref xpow l) alpha)
			  (w65x v xpow ))
		      (w65x w  xpow)))))
     ;; 2^l < n < 2^(l+1)-1 is the remaining case
     (t  (multiple-value-bind (v w)
	     (part u (1+ (- n(expt 2 l))))
	   (newname (+s (*s (aref xpow l) (w65x v xpow)) (w65x w xpow))))))))

;; driver program for non monic case

(defun win65nm(u x)
  (let* ((n (1- (length u)))
	 (l (floor (log n 2)))
	 (xpow (make-array (1+ l)))
	 (mul nil))
    
    (if (equalp (car u) 1) (return (win65 u m)))
    (if (not (> n 1)) (error "degree ~s is not interesting here" n))
    (setf mul (/ 1 (car u)))
    (framevars
     (let()
       ;; compute some powers of x in advance
	 (setf (aref xpow 0) (newname x));x^(2^0)=x
	 (loop for i from 1 to l do
	       ;;xpow[i]contains the temporary for
	       ;;x^(2^i) = x^2, x^4, ..
	       (setf (aref xpow i)
		 (newname (*s (aref xpow (1- i))
				(aref xpow (1- i))))))
	 ;; generate the rest of the program
	 (*s mul (w65x (cons 1 (mapcar #'(lambda(r)(* r mul)) u)) xpow))))))

#+ignore ;debugging version
(defun w65x (u xpow)
  (if (not (equalp (car u) 1)) (error "not monic ~s" u))
  ;; u must be MONIC
  (let* ((n (1- (length u)))		;n is degree
	 (l (floor (log n 2)))) 

    (cond
     ((= n 0)(error "not monic")(car u)); should not happen
     ;; dispense with remaining base cases.
     ((= n 1)
      (format t "~% case n=1, u= ~s" u)
      (newname
	      (+s  (aref xpow 0) (cadr u))));linear (x+a))
     ((= n 2)
      (format t "~% case n=2, u= ~s" u)
      (newname
	      (+s			;quadratic (x^2+a*x+b)
	       (aref xpow 1)		; x^2
	       (*s (aref xpow 0) (cadr u)); x*a
	       (caddr u)		; b
	       )))
     ((= n (expt 2 l))			;n is 2^l
      (format t "~%case n=2^l")
      ;; use one step of horner's rule to reduce degree
      (let ((h (reverse u)))
	(newname (+s (car h)		; last term /constant
		       (*s (aref xpow 0)
			   (w65x (reverse (cdr h)) xpow))))))
   
     ((= n (1- (expt 2 (1+ l))))	;n is 2^(l+1)-1 = 2m+1
            (format t "~%case n=2^(l+1)-1")
      ;; use synthetic division
      (let* ((m (/ (1- n) 2))
	     (alpha (1- (elt u (1+ m))))
	     (divisor (cons 1;; x^(m+1)+alpha
			    (append (make-list  m :initial-element 0)
				    (list alpha))))
	     (v nil)
	     (w nil))
	(multiple-value-setq
	    (v w) (sd u divisor))
	(newname(+s (*s (+s (aref xpow l) alpha)
			  (w65x v xpow ))
		      (w65x w  xpow)))))
     ;; 2^l < n < 2^(l+1)-1 is the remaining case
     (t 
      (format t "~%case 2^l<n<2^(l+1)-1, l=~s" l)
      (multiple-value-bind (v w)
	  (part u (1+ (- n(expt 2 l))))
	      (format t "~%v= ~s w= ~s" v w)
	      (newname (+s (*s (aref xpow l) (w65x v xpow)) (w65x w xpow))))))))


(defun win65e(u x) ;; use the same idea but keep track of errors.
  
  (let* ((n (1- (length u)))
	 (l (floor (log n 2)))
	 (xpow (make-array (1+ l))))
    
    (if (not (equalp (car u) 1)) (error "not monic ~s" u))
    (if (not (> n 1)) (error "degree ~s is not interesting here" n))
    
    (framevars
     (let()
       ;; compute some powers of x in advance
	 (setf (aref xpow 0) (newname x));x^(2^0)=x
	 (loop for i from 1 to l do
	       ;;xpow[i]contains the temporary for
	       ;;x^(2^i) = x^2, x^4, ..
	       (setf (aref xpow i)
		 (newname (*s (aref xpow (1- i))
				(aref xpow (1- i))))))
	 ;; generate the rest of the program
	(cons 'values (multiple-value-list (w65xe u xpow)))))))
      

;; inside of w65 with errors returns 2 values.  The name of the
;; variable bound to the value of u(x) and the name of the variable
;; bound to the error in u(x)

(defun w65xe (u xpow)
  (if (not (equalp (car u) 1)) (error "not monic ~s" u))
  ;; u must be MONIC
  (let* ((n (1- (length u)))		;n is degree
	 (l (floor (log n 2)))) 
    (cond
     ((= n 0)(error "not monic")(car u)); should not happen
     ;; dispense with remaining base cases.
     ((= n 1)
      (values
	(newname
      	 (+s  (aref xpow 0) (cadr u))) 	;linear (x+a)
	(newname 
	 (+s `(abs,(aref xpow 0)) `(abs (cadr)))) ;error is (|x|+|a|)* eps
      ))
      
     ((= n 2)
      (let ((ax (newname (*s (aref xpow 0) (cadr u)))))
      (values
	(newname
       (+s				;quadratic (x^2+a*x+b)
	(aref xpow 1)			; x^2
	ax				; x*a
	(caddr u)			; b
	))
	(newname 
	 (+s (aref xpow 1) (abs (caddr u)) `(abs ,ax))))) ;error bound
     )
     ((= n (expt 2 l))			;n is 2^l
      ;; use one step of horner's rule to reduce degree
      (let* ((h (reverse u)))
	(multiple-value-bind
	    (p er)
	    (w65x (reverse (cdr h)) xpow)
	  (let ((ap (newname (*s (aref xpow 0) p))))
	(values
	  (newname (+s (car h)		; last term /constant
		       ap))
	  (newname (+s `(abs, p) (abs (car h)) er)))))))

     ;;;; edited up to here. Need to add more error cases.
     ((= n (1- (expt 2 (1+ l))))	;n is 2^(l+1)-1 = 2m+1
      ;; use synthetic division
      (let* ((m (/ (1- n) 2))
	     (alpha (1- (elt u (1+ m))))
	     (divisor (cons 1;; x^(m+1)+alpha
			    (append (make-list  m :initial-element 0)
				    (list alpha))))
	     (v nil)
	     (w nil))
	(multiple-value-setq
	    (v w) (sd u divisor))
	(newname(+s (*s (+s (aref xpow l) alpha)
			  (w65x v xpow ))
		      (w65x w  xpow)))))
     ;; 2^l < n < 2^(l+1)-1 is the remaining case
     (t  (multiple-value-bind (v w)
	     (part u (1+ (- n(expt 2 l))))
	   (newname (+s (*s (aref xpow l) (w65x v xpow)) (w65x w xpow))))))))