;; Trefethen and Bau's linear algebra textbook has a nice plot
;; on page 93 demonstrating that polynomial root-finding is
;; an ill-conditioned problem. This little program shows how 
;; one would create that plot in Lush (not sure the point sizes
;; are appropriate for all terminals, it looks great with the
;; wxt terminal).

(libload "libnum/polynomials")
(libload "gnuplot/plot")

;; http://en.wikipedia.org/wiki/Wilkinson%27s_polynomial
(defparameter *p-wilkinson* [
   2432902008176640000   -8752948036761600000
  13803759753640704000  -12870931245150988800
   8037811822645051776   -3599979517947607200
   1206647803780373360    -311333643161390640
     63030812099294896     -10142299865511450
      1307535010540395       -135585182899530
        11310276995381          -756111184500
           40171771630            -1672280820
	      53327946               -1256850
	         20615                   -210
		     1])

(defparameter *roots-wilkinson*	(arange 1 20))

(defun wilkinson-demo (&optional (n-pertubations 100))
  (let* ((p *p-wilkinson*)
	 (n ($n *p-wilkinson*))
	 (s ($> [@ 1e-10] n-pertubations n))
	 (coeffs (* ($< p n-pertubations) (+ 1 (rand s))))
	 (plotter (new Gnuplot 'interactive ())) )
    (plot
      (title "Roots of Wilkinson's polynomial")
      (points (++> *roots-wilkinson* [0]) (lc 'black) (pt 7) (ps 1))
      (do ((c coeffs))
        (points (poly-solve c) (lc 'red) (pt 7) (ps 0.2)) ))
    plotter))