maxima-src 5.32.1-1 / usr / share / maxima / 5.32.1 / src / cpoly.lisp

;;; -*-  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.

(declare-top (special $partswitch $keepfloat $demoivre $listconstvars
		      $algebraic $ratfac $programmode))

(declare-top (special *logbas* *infin* *are* *mre* *cr* *ci* *sr* *si*
		      *tr* *ti* *zr* *zi* *n* *nn* *bool*
		      *conv* *pvr* *pvi* *polysc* *polysc1*))

(declare-top (special *pr-sl* *pi-sl* *shr-sl* *shi-sl* *qpr-sl* *qpi-sl* *hr-sl*
		      *hi-sl* *qhr-sl* *qhi-sl*))

(declare-top (special *u* *v* *a* *b* *c* *d* *a1* *a3* *a7* *e* *f* *g* *h*
		      *szr* *szi* *lzr* *lzi* *nz* *ui* *vi* *s*))

(defmvar $polyfactor nil
  "When T factor the polynomial over the real or complex numbers.")

(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 (intl:gettext "allroots: expected a polynomial in one variable; found variables ~M") `((mlist) ,@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
		   (let ((*cr* 0.0) (*ci* 0.0))
		     (cdivid-sl (float (car expr)) (float (cadr expr))
				(float (car den)) (float (cadr den)))
		     (return (simplify (list '(mplus)
					     (simplify (list '(mtimes) '$%i *ci*))
					     *cr*)))))
		  (t (return (list '(mlist simp)))))))
     (setq degree (cadr expr) *nn* (1+ degree))
     (setq *pr-sl* (make-array *nn* :initial-element 0.0))
     (setq *pi-sl* (make-array *nn* :initial-element 0.0))
     (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 0.0))
     (setq *shi-sl* (make-array *nn* :initial-element 0.0))
     (setq *qpr-sl* (make-array *nn* :initial-element 0.0))
     (setq *hr-sl*  (make-array degree :initial-element 0.0))
     (setq *qhr-sl* (make-array degree :initial-element 0.0))
     (setq *qpi-sl* (make-array *nn* :initial-element 0.0))
     (when complex
       (setq *hi-sl*  (make-array degree :initial-element 0.0))
       (setq *qhi-sl* (make-array degree :initial-element 0.0)))
     (setq *nn* degree)
     (if complex
	 (setq res (errset (cpoly-sl degree)))
	 (setq res (errset (rpoly-sl degree))))
     (unless res
       (mtell (intl:gettext "allroots: unexpected error; treat results with caution.~%")))
     (when (= *nn* degree)
       (merror (intl:gettext "allroots: no roots found.")))
     (setq res nil)
     (cond ((not (zerop *nn*))
	    (mtell (intl:gettext "allroots: only ~S out of ~S roots found.~%") (- degree *nn*) degree)
	    (setq expr 0.0)
	    (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 (let ((*cr* 0.0) (*ci* 0.0))
			 (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*)))
		  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 (let ((expr (simplify (list '(mequal) var expr))))
			    (if $programmode expr (displine expr)))
			  res)))))
     (return (simplify (if $polyfactor
			   (cons '(mtimes) res)
			   (cons '(mlist) (nreverse res)))))))

(defun cpoly-err (expr)
  (merror (intl:gettext "allroots: expected a polynomial; found ~M") expr))

(defun cpoly-sl (degree)
  (let ((*logbas* (log 2.0))
	(*infin* most-positive-flonum)
	(*are* flonum-epsilon)
	(*mre* 0.0)
	(xx (sqrt 0.5))
	(yy 0.0)
	(cosr (cos (float (* 94/180 pi))))
	(sinr (sin (float (* 94/180 pi))))
	(*cr* 0.0) (*ci* 0.0)
	(*sr* 0.0) (*si* 0.0)
	(*tr* 0.0) (*ti* 0.0)
	(*zr* 0.0) (*zi* 0.0)
	(bnd 0.0)
	(*n* 0)
	(*polysc* 0)
	(*polysc1* 0)
	(*conv* nil))
    (setq *mre* (* 2.0 (sqrt 2.0) *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))))
    (if (> *nn* 0) (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*))

(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*))
	      (* 10.0 *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) 0.0)
	     (setf (aref *hi-sl* 0) 0.0)))))

(defun fxshft-sl (l2)
  (let ((test t)
	(pasd nil)
	(otr 0.0) (oti 0.0)
	(svsr 0.0) (svsi 0.0)
	(*bool* nil)
	(*pvr* 0.0) (*pvi* 0.0))
    (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))

(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 ((> (* 20.0 (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.1 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) (* 10.0 *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* 0.0 *ti* 0.0)))
       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) 0.0)
	      (setf (aref *hi-sl* 0) 0.0))
	(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*)))

;; Compute a lower bound on the roots of the polynomial.  Let the
;; polynomial be sum(a[k]*x^k,k,0,n).  Then the lower bound is the
;; smallest real root of sum(|a[k]|*x^k,k,0,n-1)-a[n].  For our
;; purposes, this lower bound is computed to an accuracy of .005 or
;; so.
(defun cauchy-sl ()
  (let ((x (exp (/ (- (log (aref *shr-sl* *nn*)) (log (aref *shr-sl* 0))) (float *nn*))))
	(xm 0.0))
    (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.1 x) *f* (aref *shr-sl* 0))
      (do ((i 1 (1+ i))) ((> i *nn*)) (setq *f* (+ (aref *shr-sl* i) (* *f* xm))))
      (cond ((not (< 0.0 *f*)) (return t)))
      (setq x xm))
    (do ((dx x) (df) (*f*))
	((> 5e-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)))))

;; Scale the coefficients of the polynomial to reduce the chance of
;; overflow or underflow.  This is different from the algorithm given
;; in TOMS 419 and 493 which just scales the coefficients by some
;; fixed factor.  Instead, this algorithm computes a scale factor for
;; the roots and adjusts the coefficients of the polynomial
;; appropriately.
;;
;; The scale factor for the roots is in *polysc1*.  The roots are
;; scaled by 2^(*polysc1*).  There is an additional scaling *polysc*
;; such that the coefficients of the polynomial are all scaled by
;; 2^(*polysc*).
(defun scale-sl ()
  (do ((i 0 (1+ i)) (j 0) (x 0.0) (dx 0.0))
      ((> 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.5 (/ dx *logbas*)))
	     x (+ x (* (float (* *polysc1* *nn*)) *logbas* 0.5))
	     *polysc* (floor (+ 0.5 (/ 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)
  (let ((r1 0.0))
    (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)))))
  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)))))
	((* (sqrt 2.0)
	  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

(defun rpoly-sl (degree)
  (let ((*logbas* (log 2.0))
	(*infin* most-positive-flonum)
	(*are* flonum-epsilon)
	(*mre* 0.0)
	(xx (sqrt 0.5)) ;; sqrt(0.5)
	(yy 0.0)
	(cosr (cos (float (* 94/180 pi))))
	(sinr (sin (float (* 94/180 pi))))
	(aa 0.0)
	(cc 0.0)
	(bb 0.0)
	(bnd 0.0)
	(*sr* 0.0)
	(*u* 0.0)
	(*v* 0.0)
	(t1 0.0)
	(*szr* 0.0) (*szi* 0.0)
	(*lzr* 0.0) (*lzi* 0.0)
	(*nz* 0)
	(*n* 0)
	(*polysc* 0)
	(*polysc1* 0)
	(zerok 0)
	(conv1 t))
    (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))))
    (if (> *nn* 0) (scale-sl))
    ;; Loop to find all roots
    (do nil
	((< *nn* 3)
	 (cond ((= *nn* 2)
		;; Solve the final quadratic polynomial
		(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) 1.0))
		      (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
		;; Solve the final linear polynomial 
		(setf (aref *pr-sl* 1) (- (/ (aref *pr-sl* 1) (aref *pr-sl* 0))))))
	 (setq *nn* 0))
      ;; Calculate a lower bound on the modulus of the zeros.
      (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i))))
      (setq bnd (cauchy-sl))

      ;; Compute derivative of polynomial.  Save result in *hr-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 5 steps with no shift.
      (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) 0.0)
		     (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* 10.0)))))))
      (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* (* -2.0 *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*) 1.0)))))
	       (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)))))

(defun fxshfr-sl (l2)
  (let ((*my-type* 0)
	(*a* 0.0) (*b* 0.0) (*c* 0.0) (*d* 0.0) (*e* 0.0) (*f* 0.0) (*g* 0.0) (*h* 0.0)
	(*a1* 0.0) (*a3* 0.0) (*a7* 0.0))
    (declare (special *my-type*))
    (setq *nz* 0)
    (quadsd-sl)
    (calcsc-sl)
    (do ((j 1 (1+ j))
	 (betav 0.25)
	 (betas 0.25)
	 (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 0.0)
      (or (zerop (aref *hr-sl* *n*))
	  (setq ss (- (/ (aref *pr-sl* *nn*) (aref *hr-sl* *n*)))))
      (setq tv 1.0 ts 1.0)
      (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 1.0)
	     (and (< tv otv) (setq tvv (* tv otv)))
	     (setq tss 1.0)
	     (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.25 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.25 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))))

(defun quadit-sl nil
  (setq *nz* 0 *u* *ui* *v* *vi*)
  (do ((tried) (j 0) (ee) (mp) (relstp) (omp) (ms))
      (nil)
    (quad-sl 1.0 *u* *v*)
    (and (> (abs (- (abs *szr*) (abs *lzr*))) (* 1e-2 (abs *lzr*)))
	 (return nil))
    (quadsd-sl)
    (setf mp (+ (abs (- *a* (* *szr* *b*))) (abs (* *szi* *b*))))
    (setf ms (cmod-sl *szr* *szi*))
    (setf ee (errev-sl ms mp))
    (cond ((not (> mp (* 2e1 ee))) (setq *nz* 2)
	   (return nil)))
    (setq j (1+ j))
    (and (> j 20) (return nil))
    (cond ((not (or (< j 2) (> relstp 1e-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 (* 2e1 (- (* (+ *are* *mre*) ee) (* *mre* mp)))))
	   (setq *nz* 1 *szr* *s* *szi* 0.0)
	   (return 0)))
    (setq j (1+ j))
    (and (> j 10) (return 0))
    (cond ((not (or (< j 2)
		    (> (abs t1) (* 1e-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*)) 1e1 *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) 0.0)
	     (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 0.0)
    (and (> (abs kv) (* (abs (aref *hr-sl* *n*)) 10.0 *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*)) 1e2 *are*))
		  (> (abs *d*) (* (abs (aref *hr-sl* (1- *n*))) 1e2 *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) 0.0)
	 (setf (aref *hr-sl* 1) 0.0)
	 (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*)) 1e1 *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) 0.0)
	   (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*))
  (let ((a4 0.0) (a5 0.0)
	(b1 0.0) (b2 0.0)
	(c1 0.0) (c2 0.0) (c3 0.0) (c4 0.0))
    (cond ((= *my-type* 3)
	   (setq *ui* 0.0 *vi* 0.0))
	  (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* 0.0 *vi* 0.0)
		 (setq *ui* (- *u* (/ (+ (* *u* (+ c3 c2))
					 (* *v* (+ (* b1 *a1*) (* b2 *a7*))))
				      c1))
		       *vi* (* *v* (1+ (/ c4 c1)))))))
    nil))

;; Divide the polynomial in *pr-sl* by the quadratic x^2 + (*u*)*x +
;; (*v*).  Place the quotient in *qpr-sl* and the remainder in *a* and
;; *b*.  I (rtoy) think the remainder polynomial is (*b*)*x + (*a*).
(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)))

;; Compute the zeros of the quadratic a0*z^2+b1*z+c0.  The larger zero
;; is returned in *szr* and *szi*.  The smaller zero is in *lzr* and
;; *lzi*.
(defun quad-sl (a0 b1 c0)
  (setq *szr* 0.0 *szi* 0.0 *lzr* 0.0 *lzi* 0.0)
  (let ((b0 0.0)
	(l0 0.0)
	(*e* 0.0))
    ;; Handle the degenerate cases of a0 = 0 or c0 = 0 first.
    (cond ((zerop a0) (or (zerop b1) (setq *szr* (- (/ c0 b1)))))
	  ((zerop c0) (setq *lzr* (- (/ b1 a0))))
	  (t
	   ;; Quadratic formula.
	   (setq b0 (/ b1 2.0))
	   (cond ((< (abs b0) (abs c0))
		  (setq *e* a0)
		  (and (< c0 0.0) (setq *e* (- a0)))
		  (setq *e* (- (* b0 (/ b0 (abs c0))) *e*)
			l0 (* (sqrt (abs *e*)) (sqrt (abs c0)))))
		 (t (setq *e* (- 1.0 (* (/ a0 b0) (/ c0 b0)))
			  l0 (* (sqrt (abs *e*)) (abs b0)))))
	   (cond ((< *e* 0.0)
		  (setq *szr* (- (/ b0 a0))
			*lzr* *szr*
			*szi* (abs (/ l0 a0))
			*lzi* (- *szi*)))
		 (t (or (< b0 0.0) (setq l0 (- l0)))
		    (setq *lzr* (/ (- l0 b0) a0))
		    (or (zerop *lzr*) (setq *szr* (/ c0 *lzr* a0)))))))
    nil))

;; This is a very straightforward conversion of $allroots to use
;; bfloats instead of floats.

(defun bf-cpoly-err (expr)
  (merror (intl:gettext "bfallroots: expected a polynomial; found ~M") expr))

(defun fpzerop (x)
  (equal '(0 0) x))

;; (ar+%i*ai)/(br+%i*bi) -> cr+%i*ci.
(defun bf-cdivid-sl (ar ai br bi)
  (cond ((and (fpzerop br)
	      (fpzerop bi))
	 ;; Division by zero.  Should we do something else besides set
	 ;; both parts to be "infinity"?
	 (setq *cr* (setq *ci* *infin*)))
	((fpgreaterp (fpabs bi) (fpabs br))
	 (let ((r1 (fpquotient br bi)))
	   (setq bi (fpplus bi (fptimes* br r1))
		 br (fpplus ai (fptimes* ar r1))
		 *cr* (fpquotient br bi)
		 br (fpdifference (fptimes* ai r1) ar)
		 *ci* (fpquotient br bi))))
	(t
	 (let ((r1 (fpquotient bi br)))
	   (setq bi (fpplus br (fptimes* bi r1))
		 br (fpplus ar (fptimes* ai r1))
		 *cr* (fpquotient br bi)
		 br (fpdifference ai (fptimes* ar r1))
		 *ci* (fpquotient br bi))))))

(defun fpsqrt (x)
  (fproot (bcons x) 2))

(defun bf-cmod-sl (ar ai)
  (let ((ar (fpabs ar))
	(ai (fpabs ai)))
    (cond ((fpgreaterp ai ar)
	   (setq ar (fpquotient ar ai))
	   (fptimes* ai
		     (fpsqrt (fpplus (fpone) (fptimes* ar ar)))))
	  ((fpgreaterp ar ai)
	   (setq ai (fpquotient ai ar))
	   (fptimes* ar (fpsqrt (fpplus (fpone) (fptimes* ai ai)))))
	  ((fptimes* (fpsqrt (intofp 2))
		     ar)))))


(defun bf-calct-sl nil
  (do ((i 1 (1+ i))
       (tt)
       (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 (fpgreaterp (bf-cmod-sl hvr hvi)
			      (fptimes* (intofp 10)
					(fptimes* *are*
						  (bf-cmod-sl (aref *hr-sl* *n*)
							      (aref *hi-sl* *n*)))))))
       (cond ((not *bool*)
	      (bf-cdivid-sl (fpminus *pvr*) (fpminus *pvi*) hvr hvi)
	      (setq *tr* *cr*
		    *ti* *ci*))
	     (t (setq *tr* (intofp 0) *ti* (intofp 0))))
       nil)
    (setq tt (fpdifference (fpplus (aref *hr-sl* i)
				   (fptimes* hvr *sr*))
			   (fptimes* hvi *si*)))
    (setf (aref *qhi-sl* i)
	  (setq hvi (fpplus (aref *hi-sl* i)
			    (fpplus (fptimes* hvr *si*)
				    (fptimes* hvi *sr*)))))
    (setf (aref *qhr-sl* i) (setq hvr tt))))

(defun bf-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) (intofp 0))
	 (setf (aref *hi-sl* 0) (intofp 0)))
	(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)
		 (fpdifference (fpplus (aref *qpr-sl* j)
				       (fptimes* t1 *tr*))
			       (fptimes* t2 *ti*)))
	   (setf (aref *hi-sl* j)
		 (fpplus (aref *qpi-sl* j)
			 (fpplus (fptimes* t1 *ti*)
				 (fptimes* t2 *tr*)))))
	 (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
	 (setf (aref *hi-sl* 0) (aref *qpi-sl* 0))))
  nil)

(defun bf-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))
       (tt))
      ((> i *nn*))
    (setq tt (fpdifference (fpplus (aref *pr-sl* i) (fptimes* *pvr* *sr*))
			   (fptimes* *pvi* *si*)))
    (setf (aref *qpi-sl* i)
	  (setq *pvi* (fpplus (aref *pi-sl* i)
			      (fpplus (fptimes* *pvr* *si*)
				      (fptimes* *pvi* *sr*)))))
    (setf (aref *qpr-sl* i)
	  (setq *pvr* tt))))

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

(defun bf-cauchy-sl ()
  (let ((x (fpexp (fpquotient (fpdifference (fplog (aref *shr-sl* *nn*))
					    (fplog (aref *shr-sl* 0)))
			      (intofp *nn*))))
	(xm (intofp 0)))
    (setf (aref *shr-sl* *nn*) (fpminus (aref *shr-sl* *nn*)))
    (cond ((not (fpzerop (aref *shr-sl* *n*)))
	   (setq xm (fpminus (fpquotient (aref *shr-sl* *nn*)
					 (aref *shr-sl* *n*))))
	   (cond ((fpgreaterp x xm) (setq x xm)))))
    (do ((f))
	(nil)
      (setq xm (fptimes* (intofp 0.1) x)
	    f (aref *shr-sl* 0))
      (do ((i 1 (1+ i)))
	  ((> i *nn*))
	(setq f (fpplus (aref *shr-sl* i)
			(fptimes* f xm))))
      ;;(cond ((not (< 0.0 f)) (return t)))
      (when (fpgreaterp (intofp 0) f)
	(return t))
      (setq x xm))
    (do ((dx x)
	 (df)
	 (f))
	((fpgreaterp (intofp 5e-3)
		     (fpabs (fpquotient dx x)))
	 x)
      (setq f (aref *shr-sl* 0)
	    df f)
      (do ((i 1 (1+ i)))
	  ((> i *n*))
	(setq f (fpplus (fptimes* f x)
			(aref *shr-sl* i))
	      df (fpplus (fptimes* df x)
			 f)))
      (setq f (fpplus (fptimes* f x)
		      (aref *shr-sl* *nn*))
	    dx (fpquotient f df)
	    x (fpdifference x dx)))))

(defun bf-scale-float (bf scale)
  (destructuring-bind (mantissa exp)
      bf
    (if (zerop mantissa)
	(list mantissa 0)
	(list mantissa
	      (+ exp scale)))))

(defun bf-scale-sl ()
  (do ((i 0 (1+ i))
       (j 0)
       (x (intofp 0))
       (dx (intofp 0)))
      ((> i *nn*)
       (setq x (fpquotient x (intofp (- (1+ *nn*) j)))
	     dx (fpquotient (fpdifference (fplog (aref *shr-sl* *nn*))
					  (fplog (aref *shr-sl* 0)))
			    (intofp *nn*))
	     *polysc1* (fpentier (bcons (fpplus (cdr bfhalf)
						(fpquotient dx *logbas*))))
	     x (fpplus x (fptimes* (intofp (* *polysc1* *nn*))
				   (fptimes* *logbas*
					     (cdr bfhalf))))
	     *polysc* (fpentier (bcons (fpplus (cdr bfhalf) (fpquotient x *logbas*))))))
    (cond ((equalp (aref *shr-sl* i) (intofp 0))
	   (setq j (1+ j)))
	  (t
	   (setq x (fpplus x (fplog (aref *shr-sl* i)))))))
  (do ((i *nn* (1- i))
       (j (- *polysc*) (+ j *polysc1*)))
      ((< i 0))
    (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) j))
    (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) j)))
  nil)

(defun bf-noshft-sl (l1)
  (do ((i 0 (1+ i))
       (xni (intofp *nn*) (fpdifference xni (intofp 1)))
       (t1 (fpquotient (fpone) (intofp *nn*))))
      ((> i *n*))
    (setf (aref *hr-sl* i) (fptimes* (aref *pr-sl* i)
				     (fptimes* xni t1)))
    (setf (aref *hi-sl* i) (fptimes* (aref *pi-sl* i)
				     (fptimes* xni t1))))
  (do ((jj 1 (1+ jj)))
      ((> jj l1))
    (cond ((fpgreaterp (bf-cmod-sl (aref *hr-sl* *n*) (aref *hi-sl* *n*))
		       (fptimes* (intofp 10)
				 (fptimes* *are*
					   (bf-cmod-sl (aref *pr-sl* *n*)
						       (aref *pi-sl* *n*)))))
	   (bf-cdivid-sl (fpminus (aref *pr-sl* *nn*))
			 (fpminus (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) (fpdifference (fpplus (aref *pr-sl* j)
							  (fptimes* t1 *tr*))
						  (fptimes* t2 *ti*)))
	     (setf (aref *hi-sl* j) (fpplus (aref *pi-sl* j)
					    (fpplus (fptimes* t1 *ti*)
						    (fptimes* 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) (intofp 0))
	     (setf (aref *hi-sl* 0) (intofp 0))))))

(defun bf-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))
    (bf-polyev-sl)
    (setq mp (bf-cmod-sl *pvr* *pvi*)
	  ms (bf-cmod-sl *sr* *si*))
    (cond ((fpgreaterp (fptimes* (intofp 20) (bf-errev-sl ms mp))
		       mp)
	   (setq *conv* t
		 *zr* *sr*
		 *zi* *si*)
	   (return t)))
    (cond ((= i 1)
	   (setq omp mp))
	  ((or bool1
	       (fpgreaterp omp mp)
	       ;;(not (< relstp 0.05))
	       (fpgreaterp relstp (cdr bfhalf)))
	   (if (fpgreaterp (fptimes* (intofp 0.1) mp)
			   omp)
	       (return t)
	       (setq omp mp)))
	  (t
	   (setq tp relstp
		 bool1 t)
	   (when (fpgreaterp *are* relstp)
	     (setq tp *are*))
	   (setq r1 (fpsqrt tp)
		 *sr* (prog1
			  (fpdifference (fptimes* (fpplus (fpone) r1)
						  *sr*)
					(fptimes* r1 *si*))
			(setq *si* (fpplus (fptimes* (fpplus (fpone) r1)
						     *si*)
					   (fptimes* r1 *sr*)))))
	   (bf-polyev-sl)
	   (do ((j 1 (1+ j)))
	       ((> j 5))
	     (bf-calct-sl)
	     (bf-nexth-sl))
	   (setq omp *infin*)))
    (bf-calct-sl)
    (bf-nexth-sl)
    (bf-calct-sl)
    (or *bool*
	(setq relstp (fpquotient (bf-cmod-sl *tr* *ti*)
				 (bf-cmod-sl *sr* *si*))
	      *sr* (fpplus *sr* *tr*)
	      *si* (fpplus *si* *ti*)))))

(defun bf-fxshft-sl (l2)
  (let ((test t)
	(pasd nil)
	(otr (intofp 0))
	(oti (intofp 0))
	(svsr (intofp 0))
	(svsi (intofp 0))
	(*bool* nil)
	(*pvr* (intofp 0))
	(*pvi* (intofp 0)))
    (bf-polyev-sl)
    (setq *conv* nil)
    (bf-calct-sl)
    (do ((j 1 (1+ j)))
	((> j l2))
      (setq otr *tr*
	    oti *ti*)
      (bf-nexth-sl)
      (bf-calct-sl)
      (setq *zr* (fpplus *sr* *tr*)
	    *zi* (fpplus *si* *ti*))
      (cond ((and (not *bool*)
		  test
		  (not (= j l2)))
	     (cond ((fpgreaterp (fptimes* (cdr bfhalf) (bf-cmod-sl *zr* *zi*))
				(bf-cmod-sl (fpdifference *tr* otr)
					    (fpdifference *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*)
			   (bf-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)
			   (bf-polyev-sl)
			   (bf-calct-sl))
			  ((setq pasd t))))
		   ((setq pasd nil))))))
    (or *conv* (bf-vrshft-sl 10))
    nil))

(defun bf-cpoly-sl (degree)
  (let ( ;; Log of our floating point base.  (Do we need this much accuracy?)
	(*logbas* (fplog (intofp 2)))
	;; "Largest" bfloat.  What should we use?
	(*infin* (intofp most-positive-flonum))
	;; bfloat epsilon.  2^(-fpprec)
	(*are* (bf-scale-float (intofp 2) (- fpprec)))
	(*mre* (intofp 0))
	(xx (fproot bfhalf 2))
	(yy (intofp 0))
	;; cos(94deg).  Probably don't need full bfloat precision here.
	(cosr (intofp -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804))
	;; sin(94deg).  Probably don't need full bfloat precision here.
	(sinr (intofp 0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315))
	(*cr* (intofp 0))
	(*ci* (intofp 0))
	(*sr* (intofp 0))
	(*si* (intofp 0))
	(*tr* (intofp 0))
	(*ti* (intofp 0))
	(*zr* (intofp 0))
	(*zi* (intofp 0))
	(bnd (intofp 0))
	(*n* 0)
	(*polysc* 0)
	(*polysc1* 0)
	(*conv* nil))
    (setq *mre* (fptimes* (intofp 2)
			  (fptimes* (fpsqrt (intofp 2)) *are*))
	  yy (fpminus xx))
    (do ((i degree (1- i)))
	((not (and (equalp (intofp 0) (aref *pr-sl* i))
		   (equalp (intofp 0) (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) (bf-cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
    (if (> *nn* 0) (bf-scale-sl))
    (do ()
	((> 2 *nn*)
	 (bf-cdivid-sl (fpminus (aref *pr-sl* 1))
		       (fpminus (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) (bf-cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
      (setq bnd (bf-cauchy-sl))
      (catch 'newroot
	(do ((cnt1 1 (1+ cnt1)))
	    ((> cnt1 2))
	  (bf-noshft-sl 5)
	  (do ((cnt2 1 (1+ cnt2)))
	      ((> cnt2 9))
	    (setq xx (prog1
			 (fpdifference (fptimes* cosr xx)
				       (fptimes* sinr yy))
		       (setq yy (fpplus (fptimes* sinr xx)
					(fptimes* cosr yy))))
		  *sr* (fptimes* bnd xx)
		  *si* (fptimes* bnd yy))
	    (bf-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) (bf-scale-float (aref *pr-sl* i) *polysc1*))
      (setf (aref *pi-sl* i) (bf-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) (bf-scale-float (aref *pr-sl* i) j))
      (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) j)))
    *nn*))


(defmfun $bfallroots (expr)
  (prog (degree *nn* var res $partswitch
	 ($keepfloat t)
	 $demoivre
	 ($listconstvars t)
	 ($algebraic t) complex $ratfac den expr1)
     (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 (intl:gettext "bfallroots: expected a polynomial in one variable; found variables ~M") `((mlist) ,@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 (bf-cpoly-err expr1))))
		  (t (bf-cpoly-err expr1))))
	   (t (bf-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
			  (bf-cpoly-err expr1))))
		  (t (bf-cpoly-err expr1))))
	   (t (bf-cpoly-err expr1)))
     (cond ((= *nn* 0)
	    (cond ($polyfactor
		   (let ((*cr* (intofp 0))
			 (*ci* (intofp 0)))
		     (bf-cdivid-sl (cdr ($bfloat (car expr)))
				   (cdr ($bfloat (cadr expr)))
				   (cdr ($bfloat (car den)))
				   (cdr ($bfloat (cadr den))))
		     (return (add (bcons *cr*)
				  (mul '$%i (bcons *ci*))))))
		  (t (return (list '(mlist simp)))))))
     (setq degree (cadr expr) *nn* (1+ degree))
     (setq *pr-sl* (make-array *nn* :initial-element (intofp 0)))
     (setq *pi-sl* (make-array *nn* :initial-element (intofp 0)))
     (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) (cdr ($bfloat res))))
			   (t
			    (or (eq (car res) %i)
				(throw 'notpoly nil))
			    (setq res (cddr res))
			    (setf (aref *pi-sl* l) (cdr ($bfloat (car res))))
			    (setq res (caddr res))
			    (and res (setf (aref *pr-sl* l) (cdr ($bfloat res))))
			    (setq complex t))))))
	 ;; This should catch expressions like sin(x)-x
	 (bf-cpoly-err expr1))
     (setq *shr-sl* (make-array *nn* :initial-element (intofp 0)))
     (setq *shi-sl* (make-array *nn* :initial-element (intofp 0)))
     (setq *qpr-sl* (make-array *nn* :initial-element (intofp 0)))
     (setq *hr-sl*  (make-array degree :initial-element (intofp 0)))
     (setq *qhr-sl* (make-array degree :initial-element (intofp 0)))
     (setq *qpi-sl* (make-array *nn* :initial-element (intofp 0)))

     (when complex
       (setq *hi-sl*  (make-array degree :initial-element (intofp 0)))
       (setq *qhi-sl* (make-array degree :initial-element (intofp 0))))
     (setq *nn* degree)
     (if complex
	 (setq res (errset (bf-cpoly-sl degree)))
	 (setq res (bf-rpoly-sl degree)))
     (unless res
       (mtell (intl:gettext "bfallroots: unexpected error; treat results with caution.~%")))
     (when (= *nn* degree)
       (merror (intl:gettext "bfallroots: no roots found.")))
     (setq res nil)
     (cond ((not (zerop *nn*))
	    (mtell (intl:gettext "bfallroots: only ~S out of ~S roots found.~%") (- degree *nn*) degree)
	    (setq expr (bcons (intofp 0)))
	    (do ((i 0 (1+ i)))
		((> i *nn*))
	      (setq expr
		    (simplify
		     (list '(mplus) expr
			   (simplify (list '(mtimes)
					   (simplify (list '(mplus)
							   (simplify (list '(mtimes) '$%i
									   (bcons (aref *pi-sl* i))))
							   (bcons (aref *pr-sl* i))))
					   (simplify (list '(mexpt) var (- *nn* i)))))))))
	    (setq res (cons expr res)))
	   ($polyfactor
	    (setq expr (let ((*cr* (intofp 0))
			     (*ci* (intofp 0)))
			 (bf-cdivid-sl (aref *pr-sl* 0)
				       (aref *pi-sl* 0)
				       (cdr ($bfloat (car den)))
				       (cdr ($bfloat (cadr den))))
			 (add (bcons *cr*) (mul '$%i (bcons *ci*))))
		  res (cons expr res))))
     (do ((i degree (1- i)))
	 ((= i *nn*))
       ;; zr+%i*zi, where zr and zi parts of the root.
       (setq expr (add (bcons (aref *pr-sl* i))
		       (mul '$%i (bcons (aref *pi-sl* i)))))
       (setq res
	     (cond ($polyfactor
		    (cons (cond ((or complex (fpzerop (aref *pi-sl* i)))
				 (add var (mul -1 expr)))
				(t
				 (setq i (1- i))
				 (simplify (list '(mplus)
						 (simplify (list '(mexpt) var 2))
						 (simplify (list '(mtimes) var
								 (bcons (aref *pr-sl* i))))
						 (bcons (aref *pr-sl* (1+ i)))))))
			  res))
		   (t
		    (cons (let ((expr (simplify (list '(mequal) var expr))))
			    (if $programmode expr (displine expr)))
			  res)))))
     (return (simplify (if $polyfactor
			   (cons '(mtimes) res)
			   (cons '(mlist) (nreverse res)))))))

(defun bf-rpoly-sl (degree)
  (let ((*logbas* (fplog (intofp 2)))
	(*infin* (intofp most-positive-flonum))
	(*are* (bf-scale-float (intofp 2) (- fpprec)))
	(*mre* (intofp 0))
	(xx (fproot bfhalf 2))
	(yy (intofp 0))
	;; cos(94deg)
	(cosr (intofp
	       -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804))
	;; sin(94deg)
	(sinr (intofp
	       0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315))
	(aa (intofp 0))
	(cc (intofp 0))
	(bb (intofp 0))
	(bnd (intofp 0))
	(*sr* (intofp 0))
	(*u* (intofp 0))
	(*v* (intofp 0))
	(t1 (intofp 0))
	(*szr* (intofp 0))
	(*szi* (intofp 0))
	(*lzr* (intofp 0))
	(*lzi* (intofp 0))
	(*nz* 0)
	(*n* 0)
	(*polysc* 0)
	(*polysc1* 0)
	(zerok 0)
	(conv1 t))
    (setq *mre* *are*
	  yy (fpminus xx))
    (do ((i degree (1- i)))
	((not (fpzerop (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) (fpabs (aref *pr-sl* i))))
    (if (> *nn* 0) (bf-scale-sl))
    (do nil
	((< *nn* 3)
	 (cond ((= *nn* 2)
		(bf-quad-sl (aref *pr-sl* 0.) (aref *pr-sl* 1) (aref *pr-sl* 2))
		(cond ((and $polyfactor (not (fpzerop *szi*)))
		       (setf (aref *pr-sl* 2) (fpquotient (aref *pr-sl* 2)
							  (aref *pr-sl* 0)))
		       (setf (aref *pr-sl* 1) (fpquotient (aref *pr-sl* 1)
							  (aref *pr-sl* 0)))
		       (setf (aref *pi-sl* 2) (intofp 1)))
		      (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) (fpminus (fpquotient (aref *pr-sl* 1)
							    (aref *pr-sl* 0))))))
	 (setq *nn* 0))
      (do ((i 0 (1+ i)))
	  ((> i *nn*))
	(setf (aref *shr-sl* i) (fpabs (aref *pr-sl* i))))
      (setq bnd (bf-cauchy-sl))
      (do ((i 1 (1+ i)))
	  ((> i *n*))
	(setf (aref *hr-sl* i)
	      (fpquotient (fptimes* (intofp (- *n* i))
				    (aref *pr-sl* i))
			  (intofp *n*))))
      (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
      (setq aa (aref *pr-sl* *nn*)
	    bb (aref *pr-sl* *n*)
	    zerok (fpzerop (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) (intofp 0))
		     (setq zerok (fpzerop (aref *hr-sl* *n*))))
	      (t
	       (setq t1 (fpminus (fpquotient aa cc)))
	       (do ((j *n* (1- j)))
		   ((< j 1))
		 (setf (aref *hr-sl* j)
		       (fpplus (fptimes* t1 (aref *hr-sl* (1- j)))
			       (aref *pr-sl* j))))
	       (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
	       (setq zerok (not (fpgreaterp (fpabs (aref *hr-sl* *n*))
					    (fptimes* (fpabs bb)
						      (fptimes* *are* (intofp 10)))))))))
      (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
		     (fpdifference (fptimes* cosr xx)
				   (fptimes* sinr yy))
		   (setq yy (fpplus (fptimes* sinr xx)
				    (fptimes* cosr yy))))
	      *sr* (fptimes* bnd xx)
	      *u* (fptimes* (intofp -2) *sr*)
	      *v* bnd)
	(bf-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 (fpzerop *szi*)))
			     (setf (aref *pr-sl* *nn*) *v*)
			     (setf (aref *pr-sl* *n*) *u*)
			     (setf (aref *pi-sl* *nn*) (intofp 1))))))
	       (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 ((fpzerop (aref *pi-sl* i))
		    (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*)))
		   (t (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) (* 2 *polysc1*)))
		      (setq i (1- i))
		      (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*))))))
	  (t (do ((i (1+ *nn*) (1+ i)))
		 ((> i degree))
	       (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*))
	       (setf (aref *pi-sl* i) (bf-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) (bf-scale-float (aref *pr-sl* i) j)))
    t))

(defun bf-fxshfr-sl (l2)
  (let ((*my-type* 0)
	(*a* (intofp 0))
	(*b* (intofp 0))
	(*c* (intofp 0))
	(*d* (intofp 0))
	(*e* (intofp 0))
	(*f* (intofp 0))
	(*g* (intofp 0))
	(*h* (intofp 0))
	(*a1* (intofp 0))
	(*a3* (intofp 0))
	(*a7* (intofp 0)))
    (declare (special *my-type*))
    (setq *nz* 0)
    (bf-quadsd-sl)
    (bf-calcsc-sl)
    (do ((j 1 (1+ j))
	 (betav (intofp 0.25))
	 (betas (intofp 0.25))
	 (oss *sr*)
	 (ovv *v*)
	 (tvv) (tss) (ss) (vv) (tv) (ts) (ots) (otv)
	 (*ui*) (*vi*) (*s*) (svv) (svu) (iflag) (vpass) (spass) (vtry) (stry))
	((> j l2))
      (bf-nextk-sl)
      (bf-calcsc-sl)
      (bf-newest-sl)
      (setq vv *vi*
	    ss (intofp 0))
      (or (fpzerop (aref *hr-sl* *n*))
	  (setq ss (fpminus (fpquotient (aref *pr-sl* *nn*)
					(aref *hr-sl* *n*)))))
      (setq tv (intofp 1)
	    ts (intofp 1))
      (cond ((not (or (= j 1)
		      (= *my-type* 3)))
	     (or (fpzerop vv)
		 (setq tv (fpabs (fpquotient (fpdifference vv ovv)
					     vv))))
	     (or (fpzerop ss)
		 (setq ts (fpabs (fpquotient (fpdifference ss oss)
					     ss))))
	     (setq tvv (intofp 1))
	     (and (not (fpgreaterp tv otv))
		  (setq tvv (fptimes* tv otv)))
	     (setq tss (intofp 1))
	     (and (not (fpgreaterp ts ots))
		  (setq tss (fptimes* ts ots)))
	     (setq vpass (not (fpgreaterp tvv betav))
		   spass (not (fpgreaterp 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)
							(not (fpgreaterp tss tvv)))))
				    t)
			      (l50 nil nil))
			     (nil)
			   (cond (bool
				  (bf-quadit-sl)
				  (and (> *nz* 0) (return t))
				  (setq vtry t
					betav (fptimes* (intofp 0.25) 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 (bf-realit-sl))
				  (and (> *nz* 0) (return t))
				  (setq stry t
					betas (fptimes* (intofp 0.25) betas))
				  (cond ((zerop iflag)
					 (setq l50 t))
					(t
					 (setq *ui* (fpminus (fpplus *s* *s*))
					       *vi* (fptimes* *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))
		    (bf-quadsd-sl)
		    (bf-calcsc-sl)))))
      (setq ovv vv
	    oss ss
	    otv tv
	    ots ts))))

(defun bf-quadit-sl nil
  (setq *nz* 0 *u* *ui* *v* *vi*)
  (do ((tried)
       (j 0)
       (ee)
       (mp)
       (relstp)
       (omp)
       (ms))
      (nil)
    (bf-quad-sl (intofp 1) *u* *v*)
    (and (fpgreaterp (fpabs (fpdifference (fpabs *szr*)
					  (fpabs *lzr*)))
		     (fptimes* (intofp 1e-2) (fpabs *lzr*)))
	 (return nil))
    (bf-quadsd-sl)
    (setq mp (fpplus (fpabs (fpdifference *a* (fptimes* *szr* *b*)))
		     (fpabs (fptimes* *szi* *b*)))
	  ms (bf-cmod-sl *szr* *szi*)
	  ee (bf-errev-sl ms mp))
    (cond ((not (fpgreaterp mp (fptimes* (intofp 2e1) ee)))
	   (setq *nz* 2)
	   (return nil)))
    (setq j (1+ j))
    (and (> j 20) (return nil))
    (cond ((not (or (< j 2)
		    (fpgreaterp relstp (intofp 1e-2))
		    (not (fpgreaterp mp omp))
		    tried))
	   (and (not (fpgreaterp relstp *are*))
		(setq relstp *are*))
	   (setq relstp (fpsqrt relstp)
		 *u* (fpdifference *u* (fptimes* *u* relstp))
		 *v* (fpplus *v* (fptimes* *v* relstp)))
	   (bf-quadsd-sl)
	   (do ((i 1 (1+ i)))
	       ((> i 5))
	     (bf-calcsc-sl)
	     (bf-nextk-sl))
	   (setq tried t j 0)))
    (setq omp mp)
    (bf-calcsc-sl)
    (bf-nextk-sl)
    (bf-calcsc-sl)
    (bf-newest-sl)
    (and (fpzerop *vi*) (return nil))
    (setq relstp (fpabs (fpquotient (fpdifference *vi* *v*)
				    *vi*))
	  *u* *ui*
	  *v* *vi*)))

(defun bf-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 (fpplus (fptimes* pv *s*)
		       (aref *pr-sl* i)))
      (setf (aref *qpr-sl* i) pv))
    (setq mp (fpabs pv)
	  ms (fpabs *s*)
	  ee (fptimes* (fpquotient *mre* (fpplus *are* *mre*))
		       (fpabs (aref *qpr-sl* 0))))
    (do ((i 1 (1+ i)))
	((> i *nn*))
      (setq ee (fpplus (fptimes* ee ms)
		       (fpabs (aref *qpr-sl* i)))))
    (cond ((not (fpgreaterp mp
			    (fptimes* (intofp 2e1)
				      (fpdifference (fptimes* (fpplus *are* *mre*)
							      ee)
						    (fptimes* *mre* mp)))))
	   (setq *nz* 1 *szr* *s* *szi* (intofp 0))
	   (return 0)))
    (setq j (1+ j))
    (and (> j 10) (return 0))
    (cond ((not (or (< j 2)
		    (fpgreaterp (fpabs t1)
				(fptimes* (intofp 1e-3) (fpabs (fpdifference *s* t1))))
		    (not (fpgreaterp 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 (fpplus (fptimes* kv *s*)
		       (aref *hr-sl* i)))
      (setf (aref *qhr-sl* i) kv))
    (cond ((fpgreaterp (fpabs kv)
		       (fptimes* (fpabs (aref *hr-sl* *n*))
				 (fptimes* (intofp 1e1) *are*)))
	   (setq t1 (fpminus (fpquotient pv kv)))
	   (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
	   (do ((i 1 (1+ i)))
	       ((> i *n*))
	     (setf (aref *hr-sl* i)
		   (fpplus (fptimes* t1 (aref *qhr-sl* (1- i)))
			   (aref *qpr-sl* i)))))
	  (t
	   (setf (aref *hr-sl* 0) (intofp 0))
	   (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 (fpplus (fptimes* kv *s*)
		       (aref *hr-sl* i))))
    (setq t1 (intofp 0))
    (and (fpgreaterp (fpabs kv)
		     (fptimes* (fpabs (aref *hr-sl* *n*))
			       (fptimes* (intofp 10) *are*)))
	 (setq t1 (fpminus (fpquotient pv kv))))
    (setq *s* (fpplus *s* t1))))

(defun bf-calcsc-sl ()
  (declare (special *my-type*))
  (setq *d* (aref *hr-sl* 0))
  (setf (aref *qhr-sl* 0) *d*)
  (setq *c* (fpdifference (aref *hr-sl* 1)
			  (fptimes* *u* *d*)))
  (setf (aref *qhr-sl* 1) *c*)
  (do ((i 2 (1+ i))
       (c0))
      ((> i *n*))
    (setq c0 (fpdifference (fpdifference (aref *hr-sl* i)
					 (fptimes* *u* *c*))
			   (fptimes* *v* *d*)))
    (setf (aref *qhr-sl* i) c0)
    (setq *d* *c* *c* c0))
  (cond ((not (or (fpgreaterp (fpabs *c*)
			      (fptimes* (fpabs (aref *hr-sl* *n*))
					(fptimes* (intofp 100) *are*)))
		  (fpgreaterp (fpabs *d*)
			      (fptimes* (fpabs (aref *hr-sl* (1- *n*)))
					(fptimes* (intofp 100) *are*)))))
	 (setq *my-type* 3))
	((not (not (fpgreaterp (fpabs *d*) (fpabs *c*))))
	 (setq *my-type* 2
	       *e* (fpquotient *a* *d*)
	       *f* (fpquotient *c* *d*)
	       *g* (fptimes* *u* *b*)
	       *h* (fptimes* *v* *b*)
	       *a3* (fpplus (fptimes* (fpplus *a* *g*) *e*)
			    (fptimes* *h* (fpquotient *b* *d*)))
	       *a1* (fpdifference (fptimes* *b* *f*)
				  *a*)
	       *a7* (fpplus (fptimes* (fpplus *f* *u*) *a*)
			    *h*)))
	(t (setq *my-type* 1
		 *e* (fpquotient *a* *c*)
		 *f* (fpquotient *d* *c*)
		 *g* (fptimes* *u* *e*)
		 *h* (fptimes* *v* *b*)
		 *a3* (fpplus (fptimes* *a* *e*)
			      (fptimes* (fpplus (fpquotient *h* *c*)
						*g*)
					*b*))
		 *a1* (fpdifference *b*
				    (fptimes* *a* (fpquotient *d* *c*)))
		 *a7* (fpplus *a*
			      (fpplus (fptimes* *g* *d*)
				      (fptimes* *h* *f*))))))
  nil)

(defun bf-nextk-sl ()
  (declare (special *my-type*))
  (cond ((= *my-type* 3)
	 (setf (aref *hr-sl* 0) (intofp 0))
	 (setf (aref *hr-sl* 1) (intofp 0))
	 (do ((i 2 (1+ i)))
	     ((> i *n*))
	   (setf (aref *hr-sl* i) (aref *qhr-sl* (- i 2)))))
	((fpgreaterp (fpabs *a1*)
		     (fptimes* (fpabs (if (= *my-type* 1) *b* *a*))
			       (fptimes* (intofp 1e1) *are*)))
	 (setq *a7* (fpquotient *a7* *a1*)
	       *a3* (fpquotient *a3* *a1*))
	 (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
	 (setf (aref *hr-sl* 1)
	       (fpdifference (aref *qpr-sl* 1)
			     (fptimes* *a7* (aref *qpr-sl* 0))))
	 (do ((i 2 (1+ i)))
	     ((> i *n*))
	   (setf (aref *hr-sl* i)
		 (fpplus (fptimes* *a3* (aref *qhr-sl* (- i 2)))
			 (fpplus (fpminus (fptimes* *a7* (aref *qpr-sl* (1- i))))
				 (aref *qpr-sl* i))))))
	(t
	 (setf (aref *hr-sl* 0) (intofp 0))
	 (setf (aref *hr-sl* 1) (fpminus (fptimes* *a7* (aref *qpr-sl* 0))))
	 (do ((i 2 (1+ i)))
	     ((> i *n*))
	   (setf (aref *hr-sl* i)
		 (fpdifference (fptimes* *a3* (aref *qhr-sl* (- i 2)))
			       (fptimes* *a7* (aref *qpr-sl* (1- i))))))))
  nil)

(defun bf-newest-sl ()
  (declare (special *my-type*))
  (let ((a4 (intofp 0))
	(a5 (intofp 0))
	(b1 (intofp 0))
	(b2 (intofp 0))
	(c1 (intofp 0))
	(c2 (intofp 0))
	(c3 (intofp 0))
	(c4 (intofp 0)))
    (cond ((= *my-type* 3)
	   (setq *ui* (intofp 0)
		 *vi* (intofp 0)))
	  (t
	   (if (= *my-type* 2)
	       (setq a4 (fpplus (fptimes* (fpplus *a* *g*)
					  *f*)
				*h*)
		     a5 (fpplus (fptimes* (fpplus *f* *u*)
					  *c*)
				(fptimes* *v* *d*)))
	       (setq a4 (fpplus *a*
				(fpplus (fptimes* *u* *b*)
					(fptimes* *h* *f*)))
		     a5 (fpplus *c*
				(fptimes* (fpplus *u* (fptimes* *v* *f*))
					  *d*))))
	   (setq b1 (fpminus (fpquotient (aref *hr-sl* *n*)
					 (aref *pr-sl* *nn*)))
		 b2 (fpminus (fpquotient (fpplus (aref *hr-sl* (1- *n*))
						 (fptimes* b1 (aref *pr-sl* *n*)))
					 (aref *pr-sl* *nn*)))
		 c1 (fptimes* (fptimes* *v* b2) *a1*)
		 c2 (fptimes* b1 *a7*)
		 c3 (fptimes* (fptimes* b1 b1) *a3*)
		 c4 (fpdifference (fpdifference c1 c2) c3)
		 c1 (fpplus (fpplus a5 (fptimes* b1 a4))
			    (fpminus c4)))
	   (if (fpzerop c1)
	       (setq *ui* (intofp 0)
		     *vi* (intofp 0))
	       (setq *ui* (fpdifference
			   *u*
			   (fpquotient (fpplus (fptimes* *u* (fpplus c3 c2))
					       (fptimes* *v*
							 (fpplus (fptimes* b1 *a1*)
								 (fptimes* b2 *a7*))))
				       c1))
		     *vi* (fptimes* *v*
				    (fpplus (fpone) (fpquotient c4 c1)))))))
    nil))

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

(defun bf-quad-sl (a0 b1 c0)
  (setq *szr* (intofp 0)
	*szi* (intofp 0)
	*lzr* (intofp 0)
	*lzi* (intofp 0))
  (let ((b0 (intofp 0))
	(l0 (intofp 0))
	(*e* (intofp 0)))
    (cond ((fpzerop a0)
	   (or (fpzerop b1)
	       (setq *szr* (fpminus (fpquotient c0 b1)))))
	  ((fpzerop c0)
	   (setq *lzr* (fpminus (fpquotient b1 a0))))
	  (t
	   (setq b0 (fpquotient b1 (intofp 2)))
	   (cond ((not (fpgreaterp (fpabs b0) (fpabs c0)))
		  (setq *e* a0)
		  (and (not (fpgreaterp c0 (intofp 0)))
		       (setq *e* (fpminus a0)))
		  (setq *e* (fpdifference (fptimes* b0 (fpquotient b0 (fpabs c0)))
					  *e*)
			l0 (fptimes* (fpsqrt (fpabs *e*))
				     (fpsqrt (fpabs c0)))))
		 (t (setq *e* (fpdifference (intofp 1)
					    (fptimes* (fpquotient a0 b0)
						      (fpquotient c0 b0)))
			  l0 (fptimes* (fpsqrt (fpabs *e*))
				       (fpabs b0)))))
	   (cond ((not (fpgreaterp *e* (intofp 0)))
		  (setq *szr* (fpminus (fpquotient b0 a0))
			*lzr* *szr*
			*szi* (fpabs (fpquotient l0 a0))
			*lzi* (fpminus *szi*)))
		 (t (or (not (fpgreaterp b0 (intofp 0)))
			(setq l0 (fpminus l0)))
		    (setq *lzr* (fpquotient (fpdifference l0 b0) a0))
		    (or (fpzerop *lzr*)
			(setq *szr* (fpquotient (fpquotient c0 *lzr*) a0)))))))
    nil))