;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     The data in this file contains enhancments.                    ;;;;;
;;;                                                                    ;;;;;
;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
;;;     All rights reserved                                            ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     (c) Copyright 1981 Massachusetts Institute of Technology         ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package :maxima)

(macsyma-module cpoly)

;;;This is a lisp version of algorithm 419 from the Communications of
;;;the ACM (p 97 vol 15 feb 1972) by Jenkins and Traub.
;;;That algorithm is followed very closely.
;;;Note the following modifications: arrays are indexed from 0 instead
;;;of 1.  This means that the variables n and nn are one less than the 
;;;acm verson.  The zeros are put into the arrays pr-sl and pi-sl,
;;;rather than into their own arrays.
;;;The algorithm seems to benefit be taking are 
;;;mre 0.01 times the published values.

(defvar my-type 0)

(declare-top (special logbas infin smalno are mre cr ci sr si tr ti zr zi n nn bool
		      conv pvr pvi $partswitch $keepfloat $demoivre $listconstvars
		      $algebraic $polyfactor polysc polysc1 $ratfac $programmode))


(declare-top (special *pr-sl* *pi-sl* *shr-sl* *shi-sl* *qpr-sl* *qpi-sl* *hr-sl*
		      *hi-sl* *qhr-sl* *qhi-sl*))
       
#+clisp (eval-when (:compile-toplevel :load-toplevel :execute)
	  (ext:without-package-lock ()
	    (setq sys::*inhibit-floating-point-underflow* t)))

(defmfun $allroots (expr)
  (prog (degree nn var res $partswitch $keepfloat $demoivre $listconstvars
	 $algebraic complex $ratfac den expr1)
     (setq $keepfloat t $listconstvars t $algebraic t)
     (setq expr1 (setq expr (meqhk expr)))
     (setq var (delete '$%i (cdr ($listofvars expr)) :test #'eq))
     (or var (setq var (list (gensym))))
     (cond ((not (= (length var) 1))
	    (merror "`allroots': polynomial not univariate: ~M" var))
	   ((setq var (car var))))
     (setq expr ($rat expr '$%i var)
	   res (reverse (car (cdddar expr))))
     (do ((i (- (length res) (length (caddar expr))) (1- i)))
	 ((= i 0))
       (setq res (cdr res)))
     (setq den (cddr expr) expr (cadr expr))
;;;check denominator is a complex number
     (cond ((numberp den) (setq den (list den 0)))
	   ((eq (car den) (cadr res))
	    (setq den (cddr den))
	    (cond ((numberp (car den))
		   (cond ((null (cddr den)) (setq den (list 0 (car den))))
			 ((numberp (caddr den))
			  (setq den (list (caddr den) (car den))))
			 (t (cpoly-err expr1))))
		  (t (cpoly-err expr1))))
	   (t (cpoly-err expr1)))
;;;if the name variable has disappeared, this is caught here
     (setq nn 0)
     (cond ((numberp expr) (setq expr (list expr 0)))
	   ((eq (car expr) (car res)) (setq nn 1))
	   ((eq (car expr) (cadr res))
	    (setq expr (cddr expr))
	    (cond ((numberp (car expr))
		   (cond ((null (cddr expr)) (setq expr (list 0 (car expr))))
			 ((numberp (caddr expr))
			  (setq expr (list (caddr expr) (car expr))))
			 (t (cpoly-err expr1))))
		  (t (cpoly-err expr1))))
	   (t (cpoly-err expr1)))
     (cond ((= nn 0)
	    (cond ($polyfactor
		   ((lambda (cr ci)
		      (cdivid-sl (float (car expr)) (float (cadr expr))
				 (float (car den)) (float (cadr den)))
		      (return (simplify (list '(mplus)
					      (simplify (list '(mtimes) '$%i ci))
					      cr))))
		    0d0 0d0))
		  (t (return (list '(mlist simp)))))))
     (setq degree (cadr expr) nn (1+ degree))
     (setq *pr-sl* (make-array nn :initial-element 0d0))
     (setq *pi-sl* (make-array nn :initial-element 0d0))
     (or (catch 'notpoly
	   (errset (do ((expr (cdr expr) (cddr expr)) (l) (%i (cadr res)))
		       ((null expr))
		     (setq l (- degree (car expr)) res (cadr expr))
		     (cond ((numberp res) (setf (aref *pr-sl* l) (float res)))
			   (t 
			    (or (eq (car res) %i)
				(throw 'notpoly nil))
			    (setq res (cddr res))
			    (setf (aref *pi-sl* l) (float (car res)))
			    (setq res (caddr res))
			    (and res (setf (aref *pr-sl* l) (float res)))
			    (setq complex t))))))
	 ;;this should catch expressions like sin(x)-x
	 (cpoly-err expr1))
     (setq *shr-sl* (make-array nn :initial-element 0d0))
     (setq *shi-sl* (make-array nn :initial-element 0d0))
     (setq *qpr-sl* (make-array nn :initial-element 0d0))
     (setq *hr-sl*  (make-array degree :initial-element 0d0))
     (setq *qhr-sl* (make-array degree :initial-element 0d0))
     (when complex
       (setq *qpi-sl* (make-array nn :initial-element 0d0))
       (setq *hi-sl*  (make-array degree :initial-element 0d0))
       (setq *qhi-sl* (make-array degree :initial-element 0d0)))
     (setq nn degree)
     (if complex
	 (setq res (errset (cpoly-sl degree)))
	 (setq res (errset (rpoly-sl degree))))
     (unless res
       (mtell "~%Unexpected error. Treat results with caution."))
     (when (= nn degree)
       (merror "~%No roots found"))
     (setq res nil)
     (cond ((not (zerop nn))
	    (mtell "~%Only ~S out of ~S roots found " (- degree nn) degree)
	    (setq expr 0d0)
	    (do ((i 0 (1+ i)))
		((> i nn))
	      (setq expr
		    (simplify
		     (list '(mplus) expr
			   (simplify (list '(mtimes)
					   (simplify (list '(mplus)
							   (simplify (list '(mtimes) '$%i (aref *pi-sl* i)))
							   (aref *pr-sl* i)))
					   (simplify (list '(mexpt) var (- nn i)))))))))
	    (setq res (cons expr res)))
	   ($polyfactor
	    (setq expr ((lambda (cr ci)
			  (cdivid-sl (aref *pr-sl* 0) (aref *pi-sl* 0)
				     (float (car den))
				     (float (cadr den)))
			  (simplify (list '(mplus) (simplify (list '(mtimes) '$%i ci)) cr)))
			0d0 0d0)
		  res (cons expr res))))
     (do ((i degree (1- i)))
	 ((= i nn))
       (setq expr (simplify (list '(mplus)
				  (simplify (list '(mtimes) '$%i (aref *pi-sl* i)))
				  (aref *pr-sl* i))))
       (setq res
	     (cond ($polyfactor (cons (cond ((or complex (zerop (aref *pi-sl* i)))
					     (simplify (list '(mplus) var (simplify (list '(mminus) expr)))))
					    (t (setq i (1- i))
					       (simplify (list '(mplus)
							       (simplify (list '(mexpt) var 2))
							       (simplify (list '(mtimes) var
									       (aref *pr-sl* i)))
							       (aref *pr-sl* (1+ i))))))
				      res))
		   ((cons ((lambda (expr) (if $programmode expr (displine expr)))
			   (simplify (list '(mequal) var expr)))
			  res)))))
     (return (simplify (if $polyfactor
			   (cons '(mtimes) res)
			   (cons '(mlist) (nreverse res))))))) 

(defun cpoly-err (expr)
  (merror "`allroots': not a polynomial:~%~M" expr))

(defun cpoly-sl (degree) 
  ((lambda (logbas infin smalno are mre xx yy cosr sinr cr ci sr si tr
	    ti zr zi bnd n polysc polysc1 conv) 
     (setq mre (* 2d0 (sqrt 2d0) are)
	   yy (- xx))
     (do ((i degree (1- i)))
	 ((not (and (zerop (aref *pr-sl* i)) (zerop (aref *pi-sl* i))))
	  (setq nn i
		n (1- i))))
     (setq degree nn)
     (do ((i 0 (1+ i)))
	 ((> i nn))
       (setf (aref *shr-sl* i) (cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
     (scale-sl)
     (do ()
	 ((> 2 nn)
	  (cdivid-sl (- (aref *pr-sl* 1)) (- (aref *pi-sl* 1))
		     (aref *pr-sl* 0) (aref *pi-sl* 0))
	  (setf (aref *pr-sl* 1) cr)
	  (setf (aref *pi-sl* 1) ci)
	  (setq nn 0))
       (do ((i 0 (1+ i)))
	   ((> i nn))
	 (setf (aref *shr-sl* i) (cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
       (setq bnd (cauchy-sl))
       (catch 'newroot
	 (do ((cnt1 1 (1+ cnt1)))
	     ((> cnt1 2))
	   (noshft-sl 5)
	   (do ((cnt2 1 (1+ cnt2)))
	       ((> cnt2 9))
	     (setq xx (prog1
			  (- (* cosr xx) (* sinr yy))
			(setq yy (+ (* sinr xx) (* cosr yy)))) 
		   sr (* bnd xx) 
		   si (* bnd yy))
	     (fxshft-sl (* 10 cnt2))
	     (cond (conv (setf (aref *pr-sl* nn) zr)
			 (setf (aref *pi-sl* nn) zi)
			 (setq nn n n (1- n))
			 (do ((i 0 (1+ i)))
			     ((> i nn))
			   (setf (aref *pr-sl* i) (aref *qpr-sl* i))
			   (setf (aref *pi-sl* i) (aref *qpi-sl* i)))
			 (throw 'newroot t))))))
       (or conv (return t)))
     (do ((i (1+ nn) (1+ i)))
	 ((> i degree))
       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) polysc1))
       (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) polysc1)))
     (do ((i 0 (1+ i)) (j (- polysc (* polysc1 degree)) (+ j polysc1)))
	 ((> i nn))
       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))
       (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) j)))
     nn)
   (log 2d0)
   most-positive-double-float
   least-positive-double-float
   double-float-epsilon
   0d0 0.70710677d0 0d0 -0.069756474d0 0.99756405d0
   0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0 0 0 nil))

(defun noshft-sl (l1) 
  (do ((i 0 (1+ i))
       (xni (float nn) (1- xni))
       (t1 (/ (float nn))))
      ((> i n))
    (setf (aref *hr-sl* i) (* (aref *pr-sl* i) xni t1))
    (setf (aref *hi-sl* i) (* (aref *pi-sl* i) xni t1)))
  (do ((jj 1 (1+ jj)))
      ((> jj l1))
    (cond ((> (cmod-sl (aref *hr-sl* n) (aref *hi-sl* n))
	      (* 10d0 are (cmod-sl (aref *pr-sl* n) (aref *pi-sl* n))))
	   (cdivid-sl (- (aref *pr-sl* nn)) (- (aref *pi-sl* nn)) (aref *hr-sl* n) (aref *hi-sl* n))
	   (setq tr cr ti ci)
	   (do ((j n (1- j)) (t1) (t2))
	       ((> 1 j))
	     (setq t1 (aref *hr-sl* (1- j)) t2 (aref *hi-sl* (1- j)))
	     (setf (aref *hr-sl* j) (- (+ (aref *pr-sl* j) (* t1 tr)) (* t2 ti)))
	     (setf (aref *hi-sl* j) (+ (aref *pi-sl* j) (* t1 ti) (* t2 tr))))
	   (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
	   (setf (aref *hi-sl* 0) (aref *pi-sl* 0)))
	  (t (do ((j n (1- j)))
		 ((> 1 j))
	       (setf (aref *hr-sl* j) (aref *hr-sl* (1- j)))
	       (setf (aref *hi-sl* j) (aref *hi-sl* (1- j))))
	     (setf (aref *hr-sl* 0) 0d0)
	     (setf (aref *hi-sl* 0) 0d0))))) 

(defun fxshft-sl (l2) 
  ((lambda (test pasd otr oti svsr svsi bool pvr pvi) 
     (polyev-sl)
     (setq conv nil)
     (calct-sl)
     (do ((j 1 (1+ j)))
	 ((> j l2))
       (setq otr tr oti ti)
       (nexth-sl)
       (calct-sl)
       (setq zr (+ sr tr) zi (+ si ti))
       (cond ((and (not bool) test (not (= j l2)))
	      (cond ((> (* 0.5 (cmod-sl zr zi))
			(cmod-sl (- tr otr) (- ti oti)))
		     (cond (pasd (do ((i 0 (1+ i)))
				     ((> i n))
				   (setf (aref *shr-sl* i) (aref *hr-sl* i))
				   (setf (aref *shi-sl* i) (aref *hi-sl* i)))
				 (setq svsr sr svsi si)
				 (vrshft-sl 10.)
				 (when conv (return nil))
				 (setq test nil)
				 (do ((i 0 (1+ i)))
				     ((> i n))
				   (setf (aref *hr-sl* i) (aref *shr-sl* i))
				   (setf (aref *hi-sl* i) (aref *shi-sl* i)))
				 (setq sr svsr si svsi)
				 (polyev-sl)
				 (calct-sl))
			   ((setq pasd t))))
		    ((setq pasd nil))))))
     (or conv (vrshft-sl 10))
     nil)
   t nil 0d0 0d0 0d0 0d0 nil 0d0 0d0)) 

(defun vrshft-sl (l3) 
  (setq conv nil sr zr si zi)
  (do ((i 1 (1+ i))
       (bool1 nil)
       (mp) (ms) (omp) (relstp) (tp) (r1))
      ((> i l3))
    (polyev-sl)
    (setq mp (cmod-sl pvr pvi) ms (cmod-sl sr si))
    (cond ((> (* 20d0 (errev-sl ms mp)) mp)
	   (setq conv t zr sr zi si)
	   (return t)))
    (cond ((= i 1) (setq omp mp))
	  ((or bool1 (> omp mp) (not (< relstp 0.05)))
	   (if (> (* 0.1d0 mp) omp)
	       (return t)
	       (setq omp mp)))
	  (t (setq tp relstp bool1 t)
	     (when (> are relstp) (setq tp are))
	     (setq r1 (sqrt tp) 
		   sr (prog1
			  (- (* (1+ r1) sr) (* r1 si))
			(setq si (+ (* (1+ r1) si) (* r1 sr)))))
	     (polyev-sl)
	     (do ((j 1 (1+ j))) ((> j 5)) (calct-sl) (nexth-sl))
	     (setq omp infin)))
    (calct-sl)
    (nexth-sl)
    (calct-sl)
    (or bool
	(setq relstp (/ (cmod-sl tr ti) (cmod-sl sr si)) 
	      sr (+ sr tr) 
	      si (+ si ti))))) 

(defun calct-sl nil 
  (do ((i 1 (1+ i))
       ($t)
       (hvr (setf (aref *qhr-sl* 0) (aref *hr-sl* 0)))
       (hvi (setf (aref *qhi-sl* 0) (aref *hi-sl* 0))))
      ((> i n)
       (setq bool (not (> (cmod-sl hvr hvi) (* 10d0 are (cmod-sl (aref *hr-sl* n) (aref *hi-sl* n))))))
       (cond ((not bool) (cdivid-sl (- pvr) (- pvi) hvr hvi) (setq tr cr ti ci))
	     (t (setq tr 0d0 ti 0d0)))
       nil)
    (setq $t (- (+ (aref *hr-sl* i) (* hvr sr)) (* hvi si)))
    (setf (aref *qhi-sl* i) (setq hvi (+ (aref *hi-sl* i) (* hvr si) (* hvi sr))))
    (setf (aref *qhr-sl* i) (setq hvr $t)))) 

(defun nexth-sl ()
  (cond (bool (do ((j 1 (1+ j)))
		  ((> j n))
		(setf (aref *hr-sl* j) (aref *qhr-sl* (1- j)))
		(setf (aref *hi-sl* j) (aref *qhi-sl* (1- j))))
	      (setf (aref *hr-sl* 0) 0d0)
	      (setf (aref *hi-sl* 0) 0d0))
	(t (do ((j 1. (1+ j)) (t1) (t2))
	       ((> j n))
	     (setq t1 (aref *qhr-sl* (1- j)) t2 (aref *qhi-sl* (1- j)))
	     (setf (aref *hr-sl* j) (- (+ (aref *qpr-sl* j) (* t1 tr)) (* t2 ti)))
	     (setf (aref *hi-sl* j) (+ (aref *qpi-sl* j) (* t1 ti) (* t2 tr))))
	   (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
	   (setf (aref *hi-sl* 0) (aref *qpi-sl* 0))))
  nil) 

(defun polyev-sl ()
  (setq pvr (setf (aref *qpr-sl* 0) (aref *pr-sl* 0))
	pvi (setf (aref *qpi-sl* 0) (aref *pi-sl* 0)))
  (do ((i 1 (1+ i)) ($t))
      ((> i nn))
    (setq $t (- (+ (aref *pr-sl* i) (* pvr sr)) (* pvi si)))
    (setf (aref *qpi-sl* i) (setq pvi (+ (aref *pi-sl* i) (* pvr si) (* pvi sr))))
    (setf (aref *qpr-sl* i) (setq pvr $t)))) 

(defun errev-sl (ms mp) 
  (- (* (do ((j 0 (1+ j))
	       (e (/ (* (cmod-sl (aref *qpr-sl* 0) (aref *qpi-sl* 0)) mre) (+ are mre))))
	      ((> j nn) e)
	    (setq e (+ (cmod-sl (aref *qpr-sl* j) (aref *qpi-sl* j)) (* e ms))))
	  (+ are mre))
      (* mp mre))) 

(defun cauchy-sl ()
  ((lambda (x xm) 
     (setf (aref *shr-sl* nn) (- (aref *shr-sl* nn)))
     (cond ((not (zerop (aref *shr-sl* n)))
	    (setq xm (- (/ (aref *shr-sl* nn) (aref *shr-sl* n))))
	    (cond ((> x xm) (setq x xm)))))
     (do ((f))
	 (nil)
       (setq xm (* 0.1d0 x) f (aref *shr-sl* 0))
       (do ((i 1 (1+ i))) ((> i nn)) (setq f (+ (aref *shr-sl* i) (* f xm))))
       (cond ((not (< 0d0 f)) (return t)))
       (setq x xm))
     (do ((dx x) (df) (f))
	 ((> 5d-3 (abs (/ dx x))) x)
       (setq f (aref *shr-sl* 0) df f)
       (do ((i 1 (1+ i)))
	   ((> i n))
	 (setq f (+ (* f x) (aref *shr-sl* i)) df (+ (* df x) f)))
       (setq f (+ (* f x) (aref *shr-sl* nn)) dx (/ f df) x (- x dx))))
   (exp (/ (- (log (aref *shr-sl* nn)) (log (aref *shr-sl* 0))) (float nn)))
   0d0)) 

(defun scale-sl ()
  (do ((i 0 (1+ i)) (j 0) (x 0d0) (dx 0d0))
      ((> i nn)
       (setq x (/ x (float (- (1+ nn) j)))
	     dx (/ (- (log (aref *shr-sl* nn)) (log (aref *shr-sl* 0))) (float nn))
	     polysc1 (floor (+ 0.5d0 (/ dx logbas)))
	     x (+ x (* (float (* polysc1 nn)) logbas 0.5d0))
	     polysc (floor (+ 0.5d0 (/ x logbas)))))
    (cond ((zerop (aref *shr-sl* i)) (setq j (1+ j)))
	  (t (setq x (+ x (log (aref *shr-sl* i)))))))
  (do ((i nn (1- i)) (j (- polysc) (+ j polysc1)))
      ((< i 0))
    (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))
    (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) j))))

(defun cdivid-sl (ar ai br bi) 
  ((lambda (r1)
     (cond ((and (zerop br) (zerop bi)) (setq cr (setq ci infin)))
	   ((> (abs bi) (abs br))
	    (setq r1 (/ br bi) 
		  bi (+ bi (* br r1)) 
		  br (+ ai (* ar r1)) 
		  cr (/ br bi) 
		  br (- (* ai r1) ar) 
		  ci (/ br bi)))
	   ((setq r1 (/ bi br) 
		  bi (+ br (* bi r1)) 
		  br (+ ar (* ai r1)) 
		  cr (/ br bi) 
		  br (- ai (* ar r1)) 
		  ci (/ br bi)))))
   0d0)
  nil) 

(defun cmod-sl (ar ai) 
  (setq ar (abs ar)
	ai (abs ai))
  (cond ((> ai ar) (setq ar (/ ar ai)) (* ai (sqrt (1+ (* ar ar)))))
	((> ar ai) (setq ai (/ ai ar)) (* ar (sqrt (1+ (* ai ai)))))
	((* 1.41421357d0 ar)))) 

;;;this is the algorithm for doing real polynomials.
;;;it is algorithm 493 from acm toms vol 1 p 178 (1975) by jenkins.
;;;note that array indexing starts from 0.
;;;the names of the arrays have been changed to be the same as for cpoly.
;;;the correspondence is:  p - pr-sl, qp - qpr-sl, k - hr-sl, qk - qhr-sl,
;;;svk - shr-sl, temp - shi-sl.  the roots are put in pr-sl and pi-sl.
;;;the variable si appears not to be used here

(declare-top (special sr u v a b c d a1 a3 a7 e f g h szr szi lzr lzi are
		      mre n nn nz ui vi s $polyfactor)) 

(defun rpoly-sl (degree) 
  ((lambda (logbas infin smalno are mre xx yy cosr sinr aa cc bb bnd sr u v t1 szr
	    szi lzr lzi nz n polysc polysc1 zerok conv1) 
     (setq mre are yy (- xx))
     (do ((i degree (1- i))) ((not (zerop (aref *pr-sl* i))) (setq nn i n (1- i))))
     (setq degree nn)
     (do ((i 0 (1+ i))) ((> i nn)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i))))
     (scale-sl)
     (do nil
	 ((< nn 3)
	  (cond ((= nn 2)
		 (quad-sl (aref *pr-sl* 0.) (aref *pr-sl* 1) (aref *pr-sl* 2))
		 (cond ((and $polyfactor (not (zerop szi)))
			(setf (aref *pr-sl* 2) (/ (aref *pr-sl* 2) (aref *pr-sl* 0)))
			(setf (aref *pr-sl* 1) (/ (aref *pr-sl* 1) (aref *pr-sl* 0)))
			(setf (aref *pi-sl* 2) 1d0))
		       (t (setf (aref *pr-sl* 2) szr)
			  (setf (aref *pi-sl* 2) szi)
			  (setf (aref *pr-sl* 1) lzr)
			  (setf (aref *pi-sl* 1) lzi))))
		(t (setf (aref *pr-sl* 1) (- (/ (aref *pr-sl* 1) (aref *pr-sl* 0))))))
	  (setq nn 0))
       (do ((i 0 (1+ i))) ((> i nn)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i))))
       (setq bnd (cauchy-sl))
       (do ((i 1 (1+ i)))
	   ((> i n))
	 (setf (aref *hr-sl* i) (/ (* (float (- n i)) (aref *pr-sl* i)) (float n))))
       (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
       (setq aa (aref *pr-sl* nn) bb (aref *pr-sl* n) zerok (zerop (aref *hr-sl* n)))
       (do ((jj 1 (1+ jj)))
	   ((> jj 5.))
	 (setq cc (aref *hr-sl* n))
	 (cond (zerok (do ((j n (1- j)))
			  ((< j 1))
			(setf (aref *hr-sl* j) (aref *hr-sl* (1- j))))
		      (setf (aref *hr-sl* 0) 0d0)
		      (setq zerok (zerop (aref *hr-sl* n))))
	       (t (setq t1 (- (/ aa cc)))
		  (do ((j n (1- j)))
		      ((< j 1))
		    (setf (aref *hr-sl* j) (+ (* t1 (aref *hr-sl* (1- j))) (aref *pr-sl* j))))
		  (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
		  (setq zerok (not (> (abs (aref *hr-sl* n))
				      (* (abs bb) are 10d0)))))))
       (do ((i 0 (1+ i))) ((> i n)) (setf (aref *shi-sl* i) (aref *hr-sl* i)))
       (do ((cnt 1 (1+ cnt)))
	   ((> cnt 20.) (setq conv1 nil))
	 (setq xx (prog1
		      (- (* cosr xx) (* sinr yy))
		    (setq yy (+ (* sinr xx) (* cosr yy)))) 
	       sr (* bnd xx) 
	       u (* -2d0 sr) 
	       v bnd)
	 (fxshfr-sl (* 20 cnt))
	 (cond ((> nz 0)
		(setf (aref *pr-sl* nn) szr)
		(setf (aref *pi-sl* nn) szi)
		(cond ((= nz 2)
		       (setf (aref *pr-sl* n) lzr)
		       (setf (aref *pi-sl* n) lzi)
		       (cond ((and $polyfactor (not (zerop szi)))
			      (setf (aref *pr-sl* nn) v)
			      (setf (aref *pr-sl* n) u)
			      (setf (aref *pi-sl* nn) 1d0)))))
		(setq nn (- nn nz) n (1- nn))
		(do ((i 0 (1+ i))) ((> i nn)) (setf (aref *pr-sl* i) (aref *qpr-sl* i)))
		(return nil)))
	 (do ((i 0 (1+ i))) ((> i n)) (setf (aref *hr-sl* i) (aref *shi-sl* i))))
       (or conv1 (return nil)))
     (cond ($polyfactor
	    (do ((i degree (1- i)))
		((= i nn))
	      (cond ((zerop (aref *pi-sl* i))
		     (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) polysc1)))
		    (t (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) (* 2 polysc1)))
		       (setq i (1- i))
		       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) polysc1))))))
	   (t (do ((i (1+ nn) (1+ i)))
		  ((> i degree))
		(setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) polysc1))
		(setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) polysc1)))))
     (do ((i 0 (1+ i)) (j (- polysc (* polysc1 degree)) (+ j polysc1)))
	 ((> i nn))
       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))))
   (log 2d0)
   most-positive-double-float
   least-positive-double-float
   double-float-epsilon
   0d0 0.70710677d0 0d0 -0.069756474d0 0.99756405d0
   0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0 0 0 0 0 t))  

(defun fxshfr-sl (l2) 
  ((lambda (my-type a b c d e f g h a1 a3 a7)
     (declare (special my-type))
     (setq nz 0)
     (quadsd-sl)
     (calcsc-sl)
     (do ((j 1 (1+ j))
	  (betav 0.25d0)
	  (betas 0.25d0)
	  (oss sr)
	  (ovv v)
	  (tvv) (tss) (ss) (vv) (tv) (ts) (ots) (otv)
	  (ui) (vi) (s) (svv) (svu) (iflag) (vpass) (spass) (vtry) (stry))
	 ((> j l2))
       (nextk-sl)
       (calcsc-sl)
       (newest-sl)
       (setq vv vi
	     ss 0d0)
       (or (zerop (aref *hr-sl* n))
	   (setq ss (- (/ (aref *pr-sl* nn) (aref *hr-sl* n)))))
       (setq tv 1d0 ts 1d0)
       (cond ((not (or (= j 1) (= my-type 3)))
	      (or (zerop vv) (setq tv (abs (/ (- vv ovv) vv))))
	      (or (zerop ss) (setq ts (abs (/ (- ss oss) ss))))
	      (setq tvv 1d0)
	      (and (< tv otv) (setq tvv (* tv otv)))
	      (setq tss 1d0)
	      (and (< ts ots) (setq tss (* ts ots)))
	      (setq vpass (< tvv betav) spass (< tss betas))
	      (cond ((or spass vpass)
		     (setq svu u svv v)
		     (do ((i 0 (1+ i)))
			 ((> i n)) (setf (aref *shr-sl* i)
					  (aref *hr-sl* i)))
		     (setq s ss vtry nil stry nil)
		     (and (do ((bool (not (and spass (or (not vpass) (< tss tvv)))) t)
			       (l50 nil nil))
			      (nil)
			    (cond (bool (quadit-sl)
					(and (> nz 0) (return t))
					(setq vtry t
					      betav (* 0.25d0 betav))
					(cond ((or stry (not spass))
					       (setq l50 t))
					      (t (do ((i 0 (1+ i)))
						     ((> i n))
						   (setf (aref *hr-sl* i)
							  (aref *shr-sl* i)))))))
			    (cond ((not l50)
				   (setq iflag (realit-sl))
				   (and (> nz 0) (return t))
				   (setq stry t betas (* 0.25d0 betas))
				   (cond ((zerop iflag) (setq l50 t))
					 (t (setq ui (- (+ s s)) 
						  vi (* s s))))))
			    (cond (l50 (setq u svu v svv)
				       (do ((i 0 (1+ i)))
					   ((> i n))
					 (setf (aref *hr-sl* i)
						(aref *shr-sl* i)))
				       (and (or (not vpass) vtry)
					    (return nil)))))
			  (return nil))
		     (quadsd-sl)
		     (calcsc-sl)))))
       (setq ovv vv
	     oss ss
	     otv tv
	     ots ts)))
   0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0))

(defun quadit-sl nil 
  (setq nz 0 u ui v vi)
  (do ((tried) (j 0) (ee) (zm) (t1) (mp) (relstp) (omp))
      (nil)
    (quad-sl 1d0 u v)
    (and (> (abs (- (abs szr) (abs lzr))) (* 1d-2 (abs lzr)))
	 (return nil))
    (quadsd-sl)
    (setq mp (+ (abs (- a (* szr b))) (abs (* szi b))) 
	  zm (sqrt (abs v)) 
	  ee (* 2d0 (abs (aref *qpr-sl* 0)))
	  t1 (- (* szr b)))
    (do ((i 1 (1+ n)))
	((> i n)) (setq ee (+ (* ee zm) (abs (aref *qpr-sl* i)))))
    (setq ee (+ (* ee zm) (abs (+ a t1))) 
	  ee (- (* (+ (* 5d0 mre) (* 4d0 are)) ee)
		 (* (+ (* 5d0 mre) (* 2d0 are))
		     (+ (abs (+ a t1)) (* (abs b) zm)))
		 (* -2d0 are (abs t1))))
    (cond ((not (> mp (* 2d1 ee))) (setq nz 2)
	   (return nil)))
    (setq j (1+ j))
    (and (> j 20) (return nil))
    (cond ((not (or (< j 2) (> relstp 1d-2) (< mp omp) tried))
	   (and (< relstp are) (setq relstp are))
	   (setq relstp (sqrt relstp) 
		 u (- u (* u relstp)) 
		 v (+ v (* v relstp)))
	   (quadsd-sl)
	   (do ((i 1 (1+ i)))
	       ((> i 5)) (calcsc-sl) (nextk-sl))
	   (setq tried t j 0)))
    (setq omp mp)
    (calcsc-sl)
    (nextk-sl)
    (calcsc-sl)
    (newest-sl)
    (and (zerop vi) (return nil))
    (setq relstp (abs (/ (- vi v) vi)) u ui v vi)))

(defun realit-sl ()
  (setq nz 0)
  (do ((j 0) (pv) (ee) (ms) (mp) (kv) (t1) (omp))
      (nil)
    (setq pv (aref *pr-sl* 0))
    (setf (aref *qpr-sl* 0) pv)
    (do ((i 1 (1+ i)))
	((> i nn))
      (setq pv (+ (* pv s) (aref *pr-sl* i)))
      (setf (aref *qpr-sl* i) pv))
    (setq mp (abs pv)
	  ms (abs s)
	  ee (* (/ mre (+ are mre)) (abs (aref *qpr-sl* 0))))
    (do ((i 1 (1+ i)))
	((> i nn)) (setq ee (+ (* ee ms) (abs (aref *qpr-sl* i)))))
    (cond ((not (> mp (* 2d1 (- (* (+ are mre) ee) (* mre mp)))))
	   (setq nz 1 szr s szi 0d0)
	   (return 0)))
    (setq j (1+ j))
    (and (> j 10) (return 0))
    (cond ((not (or (< j 2)
		    (> (abs t1) (* 1d-3 (abs (- s t1))))
		    (not (> mp omp))))
	   (return 1)))
    (setq omp mp kv (aref *hr-sl* 0))
    (setf (aref *qhr-sl* 0) kv)
    (do ((i 1 (1+ i)))
	((> i n))
      (setq kv (+ (* kv s) (aref *hr-sl* i)))
      (setf (aref *qhr-sl* i) kv))
    (cond ((> (abs kv) (* (abs (aref *hr-sl* n)) 1d1 are))
	   (setq t1 (- (/ pv kv)))
	   (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
	   (do ((i 1 (1+ i)))
	       ((> i n))
	     (setf (aref *hr-sl* i)
		    (+ (* t1 (aref *qhr-sl* (1- i))) (aref *qpr-sl* i)))))
	  (t (setf (aref *hr-sl* 0) 0d0)
	     (do ((i 1 (1+ i)))
		 ((> i n)) (setf (aref *hr-sl* i) (aref *qhr-sl* (1- i))))))
    (setq kv (aref *hr-sl* 0))
    (do ((i 1 (1+ i)))
	((> i n)) (setq kv (+ (* kv s) (aref *hr-sl* i))))
    (setq t1 0d0)
    (and (> (abs kv) (* (abs (aref *hr-sl* n)) 10d0 are))
	 (setq t1 (- (/ pv kv))))
    (setq s (+ s t1))))

(defun calcsc-sl ()
  (declare (special my-type))
  (setq d (aref *hr-sl* 0))
  (setf (aref *qhr-sl* 0) d)
  (setq c (- (aref *hr-sl* 1) (* u d)))
  (setf (aref *qhr-sl* 1) c)
  (do ((i 2 (1+ i))
       (c0))
      ((> i n))
    (setq c0 (- (aref *hr-sl* i) (* u c) (* v d)))
    (setf (aref *qhr-sl* i) c0)
    (setq d c c c0))
  (cond ((not (or (> (abs c) (* (abs (aref *hr-sl* n)) 1d2 are))
		  (> (abs d) (* (abs (aref *hr-sl* (1- n))) 1d2 are))))
	 (setq my-type 3))
	((not (< (abs d) (abs c)))
	 (setq my-type 2 
	       e (/ a d) 
	       f (/ c d) 
	       g (* u b) 
	       h (* v b) 
	       a3 (+ (* (+ a g) e) (* h (/ b d))) 
	       a1 (- (* b f) a) 
	       a7 (+ (* (+ f u) a) h)))
	(t (setq my-type 1
		 e (/ a c) 
		 f (/ d c) 
		 g (* u e) 
		 h (* v b) 
		 a3 (+ (* a e) (* (+ (/ h c) g) b)) 
		 a1 (- b (* a (/ d c))) 
		 a7 (+ a (* g d) (* h f)))))
  nil)

(defun nextk-sl ()
  (declare (special my-type))
  (cond ((= my-type 3)
	 (setf (aref *hr-sl* 0) 0d0)
	 (setf (aref *hr-sl* 1) 0d0)
	 (do ((i 2 (1+ i)))
	     ((> i n)) (setf (aref *hr-sl* i) (aref *qhr-sl* (- i 2)))))
	((> (abs a1) (* (abs (if (= my-type 1) b a)) 1d1 are))
	 (setq a7 (/ a7 a1) a3 (/ a3 a1))
	 (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
	 (setf (aref *hr-sl* 1) (- (aref *qpr-sl* 1) (* a7 (aref *qpr-sl* 0))))
	 (do ((i 2 (1+ i)))
	     ((> i n))
	   (setf (aref *hr-sl* i)
		  (+ (* a3 (aref *qhr-sl* (- i 2)))
		      (- (* a7 (aref *qpr-sl* (1- i))))
		      (aref *qpr-sl* i)))))
	(t (setf (aref *hr-sl* 0) 0d0)
	   (setf (aref *hr-sl* 1) (- (* a7 (aref *qpr-sl* 0))))
	   (do ((i 2 (1+ i)))
	       ((> i n))
	     (setf (aref *hr-sl* i)
		    (- (* a3 (aref *qhr-sl* (- i 2)))
			(* a7 (aref *qpr-sl* (1- i))))))))
  nil)

(defun newest-sl ()
  (declare (special my-type))
  ((lambda (a4 a5 b1 b2 c1 c2 c3 c4) 
     (cond ((= my-type 3) (setq ui 0d0 vi 0d0))
	   (t (if (= my-type 2)
		  (setq a4 (+ (* (+ a g) f) h) 
			a5 (+ (* (+ f u) c) (* v d)))
		  (setq a4 (+ a (* u b) (* h f)) 
			a5 (+ c (* (+ u (* v f)) d))))
	      (setq b1 (- (/ (aref *hr-sl* n) (aref *pr-sl* nn))) 
		    b2 (- (/ (+ (aref *hr-sl* (1- n)) (* b1 (aref *pr-sl* n))) (aref *pr-sl* nn))) 
		    c1 (* v b2 a1) 
		    c2 (* b1 a7) 
		    c3 (* b1 b1 a3) 
		    c4 (- c1 c2 c3) 
		    c1 (+ a5 (* b1 a4) (- c4)))
	      (if (zerop c1)
		  (setq ui 0d0 vi 0d0)
		  (setq ui (- u (/ (+ (* u (+ c3 c2))
				      (* v (+ (* b1 a1) (* b2 a7))))
				   c1))
			vi (* v (1+ (/ c4 c1)))))))
     nil)
   0d0 0d0 0d0 0d0 0d0 0d0 0d0 0d0)) 

(defun quadsd-sl ()
  (setq b (aref *pr-sl* 0))
  (setf (aref *qpr-sl* 0) b)
  (setq a (- (aref *pr-sl* 1) (* u b)))
  (setf (aref *qpr-sl* 1) a)
  (do ((i 2 (1+ i))
       (c0))
      ((> i nn))
    (setq c0 (- (aref *pr-sl* i) (* u a) (* v b)))
    (setf (aref *qpr-sl* i) c0)
    (setq b a
	  a c0)))

(defun quad-sl (a0 b1 c0) 
  (setq szr 0d0 szi 0d0 lzr 0d0 lzi 0d0)
  ((lambda (b0 d0 e) 
     (cond ((zerop a0) (or (zerop b1) (setq szr (- (/ c0 b1)))))
	   ((zerop c0) (setq lzr (- (/ b1 a0))))
	   (t (setq b0 (/ b1 2d0))
	      (cond ((< (abs b0) (abs c0))
		     (setq e a0)
		     (and (< c0 0.0d0) (setq e (- a0)))
		     (setq e (- (* b0 (/ b0 (abs c0))) e) 
			   d0 (* (sqrt (abs e)) (sqrt (abs c0)))))
		    (t (setq e (- 1d0 (* (/ a0 b0) (/ c0 b0))) 
			     d0 (* (sqrt (abs e)) (abs b0)))))
	      (cond ((< e 0.0d0)
		     (setq szr (- (/ b0 a0)) 
			   lzr szr 
			   szi (abs (/ d0 a0)) 
			   lzi (- szi)))
		    (t (or (< b0 0d0) (setq d0 (- d0)))
		       (setq lzr (/ (- d0 b0) a0))
		       (or (zerop lzr) (setq szr (/ c0 lzr a0)))))))
     nil)
   0d0 0d0 0d0)) 

(declare-top (unspecial logbas infin smalno are mre cr ci sr si tr ti zr zi
			n nn bool conv pvr pvi polysc polysc1 sr u v
			a b c d a1 a3 a7 e f g h szr szi lzr lzi are mre
			n nn nz	ui vi s))