src/cpoly.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;;;     (c) Copyright 1981 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module cpoly)
  14
  15 ;;; This is a lisp version of algorithm 419 from the Communications of
  16 ;;; the ACM (p 97 vol 15 feb 1972) by Jenkins and Traub.  That
  17 ;;; algorithm is followed very closely.  Note the following
  18 ;;; modifications: arrays are indexed from 0 instead of 1.  This means
  19 ;;; that the variables n and nn are one less than the acm verson.  The
  20 ;;; zeros are put into the arrays pr-sl and pi-sl, rather than into
  21 ;;; their own arrays.  The algorithm seems to benefit be taking are
  22 ;;; mre 0.01 times the published values.
  23
  24 (declare-top (special *logbas* *infin* *are* *mre* *cr* *ci* *sr* *si*
  25                       *tr* *ti* *zr* *zi* *n* *nn* *bool*
  26                       *conv* *pvr* *pvi* *polysc* *polysc1*))
  27
  28 (declare-top (special *pr-sl* *pi-sl* *shr-sl* *shi-sl* *qpr-sl* *qpi-sl* *hr-sl*
  29                       *hi-sl* *qhr-sl* *qhi-sl*))
  30
  31 (declare-top (special *u* *v* *a* *b* *c* *d* *a1* *a3* *a7* *e* *f* *g* *h*
  32                       *szr* *szi* *lzr* *lzi* *nz* *ui* *vi* *s*))
  33
  34 (defmfun $allroots (expr)
  35   (prog (degree *nn* var res $partswitch $keepfloat $demoivre $listconstvars
  36          $algebraic complex $ratfac den expr1)
  37      (setq $keepfloat t $listconstvars t $algebraic t)
  38      (setq expr1 (setq expr (meqhk expr)))
  39      (setq var (delete '$%i (cdr ($listofvars expr)) :test #'eq))
  40      (or var (setq var (list (gensym))))
  41      (cond ((not (= (length var) 1))
  42             (merror (intl:gettext "allroots: expected a polynomial in one variable; found variables ~M") `((mlist) ,@var)))
  43            ((setq var (car var))))
  44      (setq expr ($rat expr '$%i var)
  45            res (reverse (car (cdddar expr))))
  46      (do ((i (- (length res) (length (caddar expr))) (1- i)))
  47          ((= i 0))
  48        (setq res (cdr res)))
  49      (setq den (cddr expr) expr (cadr expr))
  50 ;;;check denominator is a complex number
  51      (cond ((numberp den) (setq den (list den 0)))
  52            ((eq (car den) (cadr res))
  53             (setq den (cddr den))
  54             (cond ((numberp (car den))
  55                    (cond ((null (cddr den)) (setq den (list 0 (car den))))
  56                          ((numberp (caddr den))
  57                           (setq den (list (caddr den) (car den))))
  58                          (t (cpoly-err expr1))))
  59                   (t (cpoly-err expr1))))
  60            (t (cpoly-err expr1)))
  61 ;;;if the name variable has disappeared, this is caught here
  62      (setq *nn* 0)
  63      (cond ((numberp expr) (setq expr (list expr 0)))
  64            ((eq (car expr) (car res)) (setq *nn* 1))
  65            ((eq (car expr) (cadr res))
  66             (setq expr (cddr expr))
  67             (cond ((numberp (car expr))
  68                    (cond ((null (cddr expr)) (setq expr (list 0 (car expr))))
  69                          ((numberp (caddr expr))
  70                           (setq expr (list (caddr expr) (car expr))))
  71                          (t (cpoly-err expr1))))
  72                   (t (cpoly-err expr1))))
  73            (t (cpoly-err expr1)))
  74      (cond ((= *nn* 0)
  75             (cond ($polyfactor
  76                    (let ((*cr* 0.0) (*ci* 0.0))
  77                      (cdivid-sl (float (car expr)) (float (cadr expr))
  78                                 (float (car den)) (float (cadr den)))
  79                      (return (simplify (list '(mplus)
  80                                              (simplify (list '(mtimes) '$%i *ci*))
  81                                              *cr*)))))
  82                   (t (return (list '(mlist simp)))))))
  83      (setq degree (cadr expr) *nn* (1+ degree))
  84      (setq *pr-sl* (make-array *nn* :initial-element 0.0))
  85      (setq *pi-sl* (make-array *nn* :initial-element 0.0))
  86      (or (catch 'notpoly
  87            (errset (do ((expr (cdr expr) (cddr expr)) (l) (%i (cadr res)))
  88                        ((null expr))
  89                      (setq l (- degree (car expr)) res (cadr expr))
  90                      (cond ((numberp res) (setf (aref *pr-sl* l) (float res)))
  91                            (t
  92                             (or (eq (car res) %i)
  93                                 (throw 'notpoly nil))
  94                             (setq res (cddr res))
  95                             (setf (aref *pi-sl* l) (float (car res)))
  96                             (setq res (caddr res))
  97                             (and res (setf (aref *pr-sl* l) (float res)))
  98                             (setq complex t))))))
  99          ;;this should catch expressions like sin(x)-x
 100          (cpoly-err expr1))
 101      (setq *shr-sl* (make-array *nn* :initial-element 0.0))
 102      (setq *shi-sl* (make-array *nn* :initial-element 0.0))
 103      (setq *qpr-sl* (make-array *nn* :initial-element 0.0))
 104      (setq *hr-sl*  (make-array degree :initial-element 0.0))
 105      (setq *qhr-sl* (make-array degree :initial-element 0.0))
 106      (setq *qpi-sl* (make-array *nn* :initial-element 0.0))
 107      (when complex
 108        (setq *hi-sl*  (make-array degree :initial-element 0.0))
 109        (setq *qhi-sl* (make-array degree :initial-element 0.0)))
 110      (setq *nn* degree)
 111      (if complex
 112          (setq res (errset (cpoly-sl degree)))
 113          (setq res (errset (rpoly-sl degree))))
 114      (unless res
 115        (mtell (intl:gettext "allroots: unexpected error; treat results with caution.~%")))
 116      (when (= *nn* degree)
 117        (merror (intl:gettext "allroots: no roots found.")))
 118      (setq res nil)
 119      (cond ((not (zerop *nn*))
 120             (mtell (intl:gettext "allroots: only ~S out of ~S roots found.~%") (- degree *nn*) degree)
 121             (setq expr 0.0)
 122             (do ((i 0 (1+ i)))
 123                 ((> i *nn*))
 124               (setq expr
 125                     (simplify
 126                      (list '(mplus) expr
 127                            (simplify (list '(mtimes)
 128                                            (simplify (list '(mplus)
 129                                                            (simplify (list '(mtimes) '$%i (aref *pi-sl* i)))
 130                                                            (aref *pr-sl* i)))
 131                                            (simplify (list '(mexpt) var (- *nn* i)))))))))
 132             (setq res (cons expr res)))
 133            ($polyfactor
 134             (setq expr (let ((*cr* 0.0) (*ci* 0.0))
 135                          (cdivid-sl (aref *pr-sl* 0) (aref *pi-sl* 0)
 136                                     (float (car den))
 137                                     (float (cadr den)))
 138                          (simplify (list '(mplus) (simplify (list '(mtimes) '$%i *ci*)) *cr*)))
 139                   res (cons expr res))))
 140      (do ((i degree (1- i)))
 141          ((= i *nn*))
 142        (setq expr (simplify (list '(mplus)
 143                                   (simplify (list '(mtimes) '$%i (aref *pi-sl* i)))
 144                                   (aref *pr-sl* i))))
 145        (setq res
 146              (cond ($polyfactor (cons (cond ((or complex (zerop (aref *pi-sl* i)))
 147                                              (simplify (list '(mplus) var (simplify (list '(mminus) expr)))))
 148                                             (t (setq i (1- i))
 149                                                (simplify (list '(mplus)
 150                                                                (simplify (list '(mexpt) var 2))
 151                                                                (simplify (list '(mtimes) var
 152                                                                                (aref *pr-sl* i)))
 153                                                                (aref *pr-sl* (1+ i))))))
 154                                       res))
 155                    ((cons (let ((expr (simplify (list '(mequal) var expr))))
 156                             (if $programmode expr (displine expr)))
 157                           res)))))
 158      (return (simplify (if $polyfactor
 159                            (cons '(mtimes) res)
 160                            (cons '(mlist) (nreverse res)))))))
 161
 162 (defun cpoly-err (expr)
 163   (merror (intl:gettext "allroots: expected a polynomial; found ~M") expr))
 164
 165 (defun cpoly-sl (degree)
 166   (let ((*logbas* (log 2.0))
 167         (*infin* +most-positive-flonum+)
 168         (*are* +flonum-epsilon+)
 169         (*mre* 0.0)
 170         (xx (sqrt 0.5))
 171         (yy 0.0)
 172         (cosr (cos (float (* 94/180 pi))))
 173         (sinr (sin (float (* 94/180 pi))))
 174         (*cr* 0.0) (*ci* 0.0)
 175         (*sr* 0.0) (*si* 0.0)
 176         (*tr* 0.0) (*ti* 0.0)
 177         (*zr* 0.0) (*zi* 0.0)
 178         (bnd 0.0)
 179         (*n* 0)
 180         (*polysc* 0)
 181         (*polysc1* 0)
 182         (*conv* nil))
 183     (setq *mre* (* 2.0 (sqrt 2.0) *are*)
 184           yy (- xx))
 185     (do ((i degree (1- i)))
 186         ((not (and (zerop (aref *pr-sl* i)) (zerop (aref *pi-sl* i))))
 187          (setq *nn* i
 188                *n* (1- i))))
 189     (setq degree *nn*)
 190     (do ((i 0 (1+ i)))
 191         ((> i *nn*))
 192       (setf (aref *shr-sl* i) (cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
 193     (if (> *nn* 0) (scale-sl))
 194     (do ()
 195         ((> 2 *nn*)
 196          (cdivid-sl (- (aref *pr-sl* 1)) (- (aref *pi-sl* 1))
 197                     (aref *pr-sl* 0) (aref *pi-sl* 0))
 198          (setf (aref *pr-sl* 1) *cr*)
 199          (setf (aref *pi-sl* 1) *ci*)
 200          (setq *nn* 0))
 201       (do ((i 0 (1+ i)))
 202           ((> i *nn*))
 203         (setf (aref *shr-sl* i) (cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
 204       (setq bnd (cauchy-sl))
 205       (catch 'newroot
 206         (do ((cnt1 1 (1+ cnt1)))
 207             ((> cnt1 2))
 208           (noshft-sl 5)
 209           (do ((cnt2 1 (1+ cnt2)))
 210               ((> cnt2 9))
 211             (setq xx (prog1
 212                          (- (* cosr xx) (* sinr yy))
 213                        (setq yy (+ (* sinr xx) (* cosr yy))))
 214                   *sr* (* bnd xx)
 215                   *si* (* bnd yy))
 216             (fxshft-sl (* 10 cnt2))
 217             (cond (*conv* (setf (aref *pr-sl* *nn*) *zr*)
 218                           (setf (aref *pi-sl* *nn*) *zi*)
 219                           (setq *nn* *n* *n* (1- *n*))
 220                           (do ((i 0 (1+ i)))
 221                               ((> i *nn*))
 222                             (setf (aref *pr-sl* i) (aref *qpr-sl* i))
 223                             (setf (aref *pi-sl* i) (aref *qpi-sl* i)))
 224                           (throw 'newroot t))))))
 225       (or *conv* (return t)))
 226     (do ((i (1+ *nn*) (1+ i)))
 227         ((> i degree))
 228       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*))
 229       (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) *polysc1*)))
 230     (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
 231         ((> i *nn*))
 232       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))
 233       (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) j)))
 234     *nn*))
 235
 236 (defun noshft-sl (l1)
 237   (do ((i 0 (1+ i))
 238        (xni (float *nn*) (1- xni))
 239        (t1 (/ (float *nn*))))
 240       ((> i *n*))
 241     (setf (aref *hr-sl* i) (* (aref *pr-sl* i) xni t1))
 242     (setf (aref *hi-sl* i) (* (aref *pi-sl* i) xni t1)))
 243   (do ((jj 1 (1+ jj)))
 244       ((> jj l1))
 245     (cond ((> (cmod-sl (aref *hr-sl* *n*) (aref *hi-sl* *n*))
 246               (* 10.0 *are* (cmod-sl (aref *pr-sl* *n*) (aref *pi-sl* *n*))))
 247            (cdivid-sl (- (aref *pr-sl* *nn*)) (- (aref *pi-sl* *nn*)) (aref *hr-sl* *n*) (aref *hi-sl* *n*))
 248            (setq *tr* *cr* *ti* *ci*)
 249            (do ((j *n* (1- j)) (t1) (t2))
 250                ((> 1 j))
 251              (setq t1 (aref *hr-sl* (1- j)) t2 (aref *hi-sl* (1- j)))
 252              (setf (aref *hr-sl* j) (- (+ (aref *pr-sl* j) (* t1 *tr*)) (* t2 *ti*)))
 253              (setf (aref *hi-sl* j) (+ (aref *pi-sl* j) (* t1 *ti*) (* t2 *tr*))))
 254            (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
 255            (setf (aref *hi-sl* 0) (aref *pi-sl* 0)))
 256           (t (do ((j *n* (1- j)))
 257                  ((> 1 j))
 258                (setf (aref *hr-sl* j) (aref *hr-sl* (1- j)))
 259                (setf (aref *hi-sl* j) (aref *hi-sl* (1- j))))
 260              (setf (aref *hr-sl* 0) 0.0)
 261              (setf (aref *hi-sl* 0) 0.0)))))
 262
 263 (defun fxshft-sl (l2)
 264   (let ((test t)
 265         (pasd nil)
 266         (otr 0.0) (oti 0.0)
 267         (svsr 0.0) (svsi 0.0)
 268         (*bool* nil)
 269         (*pvr* 0.0) (*pvi* 0.0))
 270     (polyev-sl)
 271     (setq *conv* nil)
 272     (calct-sl)
 273     (do ((j 1 (1+ j)))
 274         ((> j l2))
 275       (setq otr *tr* oti *ti*)
 276       (nexth-sl)
 277       (calct-sl)
 278       (setq *zr* (+ *sr* *tr*) *zi* (+ *si* *ti*))
 279       (cond ((and (not *bool*) test (not (= j l2)))
 280              (cond ((> (* 0.5 (cmod-sl *zr* *zi*))
 281                        (cmod-sl (- *tr* otr) (- *ti* oti)))
 282                     (cond (pasd (do ((i 0 (1+ i)))
 283                                     ((> i *n*))
 284                                   (setf (aref *shr-sl* i) (aref *hr-sl* i))
 285                                   (setf (aref *shi-sl* i) (aref *hi-sl* i)))
 286                                 (setq svsr *sr* svsi *si*)
 287                                 (vrshft-sl 10.)
 288                                 (when *conv* (return nil))
 289                                 (setq test nil)
 290                                 (do ((i 0 (1+ i)))
 291                                     ((> i *n*))
 292                                   (setf (aref *hr-sl* i) (aref *shr-sl* i))
 293                                   (setf (aref *hi-sl* i) (aref *shi-sl* i)))
 294                                 (setq *sr* svsr *si* svsi)
 295                                 (polyev-sl)
 296                                 (calct-sl))
 297                           ((setq pasd t))))
 298                    ((setq pasd nil))))))
 299     (or *conv* (vrshft-sl 10))
 300     nil))
 301
 302 (defun vrshft-sl (l3)
 303   (setq *conv* nil *sr* *zr* *si* *zi*)
 304   (do ((i 1 (1+ i))
 305        (bool1 nil)
 306        (mp) (ms) (omp) (relstp) (tp) (r1))
 307       ((> i l3))
 308     (polyev-sl)
 309     (setq mp (cmod-sl *pvr* *pvi*) ms (cmod-sl *sr* *si*))
 310     (cond ((> (* 20.0 (errev-sl ms mp)) mp)
 311            (setq *conv* t *zr* *sr* *zi* *si*)
 312            (return t)))
 313     (cond ((= i 1) (setq omp mp))
 314           ((or bool1 (> omp mp) (not (< relstp 0.05)))
 315            (if (> (* 0.1 mp) omp)
 316                (return t)
 317                (setq omp mp)))
 318           (t (setq tp relstp bool1 t)
 319              (when (> *are* relstp) (setq tp *are*))
 320              (setq r1 (sqrt tp)
 321                    *sr* (prog1
 322                           (- (* (1+ r1) *sr*) (* r1 *si*))
 323                         (setq *si* (+ (* (1+ r1) *si*) (* r1 *sr*)))))
 324              (polyev-sl)
 325              (do ((j 1 (1+ j))) ((> j 5)) (calct-sl) (nexth-sl))
 326              (setq omp *infin*)))
 327     (calct-sl)
 328     (nexth-sl)
 329     (calct-sl)
 330     (or *bool*
 331         (setq relstp (/ (cmod-sl *tr* *ti*) (cmod-sl *sr* *si*))
 332               *sr* (+ *sr* *tr*)
 333               *si* (+ *si* *ti*)))))
 334
 335 (defun calct-sl nil
 336   (do ((i 1 (1+ i))
 337        ($t)
 338        (hvr (setf (aref *qhr-sl* 0) (aref *hr-sl* 0)))
 339        (hvi (setf (aref *qhi-sl* 0) (aref *hi-sl* 0))))
 340       ((> i *n*)
 341        (setq *bool* (not (> (cmod-sl hvr hvi) (* 10.0 *are* (cmod-sl (aref *hr-sl* *n*) (aref *hi-sl* *n*))))))
 342        (cond ((not *bool*) (cdivid-sl (- *pvr*) (- *pvi*) hvr hvi) (setq *tr* *cr* *ti* *ci*))
 343              (t (setq *tr* 0.0 *ti* 0.0)))
 344        nil)
 345     (setq $t (- (+ (aref *hr-sl* i) (* hvr *sr*)) (* hvi *si*)))
 346     (setf (aref *qhi-sl* i) (setq hvi (+ (aref *hi-sl* i) (* hvr *si*) (* hvi *sr*))))
 347     (setf (aref *qhr-sl* i) (setq hvr $t))))
 348
 349 (defun nexth-sl ()
 350   (cond (*bool* (do ((j 1 (1+ j)))
 351                   ((> j *n*))
 352                 (setf (aref *hr-sl* j) (aref *qhr-sl* (1- j)))
 353                 (setf (aref *hi-sl* j) (aref *qhi-sl* (1- j))))
 354               (setf (aref *hr-sl* 0) 0.0)
 355               (setf (aref *hi-sl* 0) 0.0))
 356         (t (do ((j 1. (1+ j)) (t1) (t2))
 357                ((> j *n*))
 358              (setq t1 (aref *qhr-sl* (1- j)) t2 (aref *qhi-sl* (1- j)))
 359              (setf (aref *hr-sl* j) (- (+ (aref *qpr-sl* j) (* t1 *tr*)) (* t2 *ti*)))
 360              (setf (aref *hi-sl* j) (+ (aref *qpi-sl* j) (* t1 *ti*) (* t2 *tr*))))
 361            (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
 362            (setf (aref *hi-sl* 0) (aref *qpi-sl* 0))))
 363   nil)
 364
 365 (defun polyev-sl ()
 366   (setq *pvr* (setf (aref *qpr-sl* 0) (aref *pr-sl* 0))
 367         *pvi* (setf (aref *qpi-sl* 0) (aref *pi-sl* 0)))
 368   (do ((i 1 (1+ i)) ($t))
 369       ((> i *nn*))
 370     (setq $t (- (+ (aref *pr-sl* i) (* *pvr* *sr*)) (* *pvi* *si*)))
 371     (setf (aref *qpi-sl* i) (setq *pvi* (+ (aref *pi-sl* i) (* *pvr* *si*) (* *pvi* *sr*))))
 372     (setf (aref *qpr-sl* i) (setq *pvr* $t))))
 373
 374 (defun errev-sl (ms mp)
 375   (- (* (do ((j 0 (1+ j))
 376                (*e* (/ (* (cmod-sl (aref *qpr-sl* 0) (aref *qpi-sl* 0)) *mre*) (+ *are* *mre*))))
 377               ((> j *nn*) *e*)
 378             (setq *e* (+ (cmod-sl (aref *qpr-sl* j) (aref *qpi-sl* j)) (* *e* ms))))
 379           (+ *are* *mre*))
 380       (* mp *mre*)))
 381
 382 ;; Compute a lower bound on the roots of the polynomial.  Let the
 383 ;; polynomial be sum(a[k]*x^k,k,0,n).  Then the lower bound is the
 384 ;; smallest real root of sum(|a[k]|*x^k,k,0,n-1)-a[n].  For our
 385 ;; purposes, this lower bound is computed to an accuracy of .005 or
 386 ;; so.
 387 (defun cauchy-sl ()
 388   (let ((x (exp (/ (- (log (aref *shr-sl* *nn*)) (log (aref *shr-sl* 0))) (float *nn*))))
 389         (xm 0.0))
 390     (setf (aref *shr-sl* *nn*) (- (aref *shr-sl* *nn*)))
 391     (cond ((not (zerop (aref *shr-sl* *n*)))
 392            (setq xm (- (/ (aref *shr-sl* *nn*) (aref *shr-sl* *n*))))
 393            (cond ((> x xm) (setq x xm)))))
 394     (do ((*f*))
 395         (nil)
 396       (setq xm (* 0.1 x) *f* (aref *shr-sl* 0))
 397       (do ((i 1 (1+ i))) ((> i *nn*)) (setq *f* (+ (aref *shr-sl* i) (* *f* xm))))
 398       (cond ((not (< 0.0 *f*)) (return t)))
 399       (setq x xm))
 400     (do ((dx x) (df) (*f*))
 401         ((> 5e-3 (abs (/ dx x))) x)
 402       (setq *f* (aref *shr-sl* 0) df *f*)
 403       (do ((i 1 (1+ i)))
 404           ((> i *n*))
 405         (setq *f* (+ (* *f* x) (aref *shr-sl* i)) df (+ (* df x) *f*)))
 406       (setq *f* (+ (* *f* x) (aref *shr-sl* *nn*)) dx (/ *f* df) x (- x dx)))))
 407
 408 ;; Scale the coefficients of the polynomial to reduce the chance of
 409 ;; overflow or underflow.  This is different from the algorithm given
 410 ;; in TOMS 419 and 493 which just scales the coefficients by some
 411 ;; fixed factor.  Instead, this algorithm computes a scale factor for
 412 ;; the roots and adjusts the coefficients of the polynomial
 413 ;; appropriately.
 414 ;;
 415 ;; The scale factor for the roots is in *polysc1*.  The roots are
 416 ;; scaled by 2^(*polysc1*).  There is an additional scaling *polysc*
 417 ;; such that the coefficients of the polynomial are all scaled by
 418 ;; 2^(*polysc*).
 419 (defun scale-sl ()
 420   (do ((i 0 (1+ i)) (j 0) (x 0.0) (dx 0.0))
 421       ((> i *nn*)
 422        (setq x (/ x (float (- (1+ *nn*) j)))
 423              dx (/ (- (log (aref *shr-sl* *nn*)) (log (aref *shr-sl* 0))) (float *nn*))
 424              *polysc1* (floor (+ 0.5 (/ dx *logbas*)))
 425              x (+ x (* (float (* *polysc1* *nn*)) *logbas* 0.5))
 426              *polysc* (floor (+ 0.5 (/ x *logbas*)))))
 427     (cond ((zerop (aref *shr-sl* i)) (setq j (1+ j)))
 428           (t (setq x (+ x (log (aref *shr-sl* i)))))))
 429   (do ((i *nn* (1- i)) (j (- *polysc*) (+ j *polysc1*)))
 430       ((< i 0))
 431     (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))
 432     (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) j))))
 433
 434 (defun cdivid-sl (ar ai br bi)
 435   (let ((r1 0.0))
 436     (cond ((and (zerop br) (zerop bi))
 437            (setq *cr* (setq *ci* *infin*)))
 438           ((> (abs bi) (abs br))
 439            (setq r1 (/ br bi)
 440                  bi (+ bi (* br r1))
 441                  br (+ ai (* ar r1))
 442                  *cr* (/ br bi)
 443                  br (- (* ai r1) ar)
 444                  *ci* (/ br bi)))
 445           ((setq r1 (/ bi br)
 446                  bi (+ br (* bi r1))
 447                  br (+ ar (* ai r1))
 448                  *cr* (/ br bi)
 449                  br (- ai (* ar r1))
 450                  *ci* (/ br bi)))))
 451   nil)
 452
 453 (defun cmod-sl (ar ai)
 454   (setq ar (abs ar)
 455         ai (abs ai))
 456   (cond ((> ai ar) (setq ar (/ ar ai)) (* ai (sqrt (1+ (* ar ar)))))
 457         ((> ar ai) (setq ai (/ ai ar)) (* ar (sqrt (1+ (* ai ai)))))
 458         ((* (sqrt 2.0)
 459           ar))))
 460
 461 ;;; This is the algorithm for doing real polynomials.  It is algorithm
 462 ;;; 493 from acm toms vol 1 p 178 (1975) by jenkins.  Note that array
 463 ;;; indexing starts from 0.  The names of the arrays have been changed
 464 ;;; to be the same as for cpoly.  The correspondence is: p - pr-sl, qp
 465 ;;; - qpr-sl, k - hr-sl, qk - qhr-sl, svk - shr-sl, temp - shi-sl.
 466 ;;; the roots *are* put in pr-sl and pi-sl.  The variable *si* appears
 467 ;;; not to be used here
 468
 469 (defun rpoly-sl (degree)
 470   (let ((*logbas* (log 2.0))
 471         (*infin* +most-positive-flonum+)
 472         (*are* +flonum-epsilon+)
 473         (*mre* 0.0)
 474         (xx (sqrt 0.5)) ;; sqrt(0.5)
 475         (yy 0.0)
 476         (cosr (cos (float (* 94/180 pi))))
 477         (sinr (sin (float (* 94/180 pi))))
 478         (aa 0.0)
 479         (cc 0.0)
 480         (bb 0.0)
 481         (bnd 0.0)
 482         (*sr* 0.0)
 483         (*u* 0.0)
 484         (*v* 0.0)
 485         (t1 0.0)
 486         (*szr* 0.0) (*szi* 0.0)
 487         (*lzr* 0.0) (*lzi* 0.0)
 488         (*nz* 0)
 489         (*n* 0)
 490         (*polysc* 0)
 491         (*polysc1* 0)
 492         (zerok 0)
 493         (conv1 t))
 494     (setq *mre* *are* yy (- xx))
 495     (do ((i degree (1- i))) ((not (zerop (aref *pr-sl* i))) (setq *nn* i *n* (1- i))))
 496     (setq degree *nn*)
 497     (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i))))
 498     (if (> *nn* 0) (scale-sl))
 499     ;; Loop to find all roots
 500     (do nil
 501         ((< *nn* 3)
 502          (cond ((= *nn* 2)
 503                 ;; Solve the final quadratic polynomial
 504                 (quad-sl (aref *pr-sl* 0.) (aref *pr-sl* 1) (aref *pr-sl* 2))
 505                 (cond ((and $polyfactor (not (zerop *szi*)))
 506                        (setf (aref *pr-sl* 2) (/ (aref *pr-sl* 2) (aref *pr-sl* 0)))
 507                        (setf (aref *pr-sl* 1) (/ (aref *pr-sl* 1) (aref *pr-sl* 0)))
 508                        (setf (aref *pi-sl* 2) 1.0))
 509                       (t (setf (aref *pr-sl* 2) *szr*)
 510                          (setf (aref *pi-sl* 2) *szi*)
 511                          (setf (aref *pr-sl* 1) *lzr*)
 512                          (setf (aref *pi-sl* 1) *lzi*))))
 513                (t
 514                 ;; Solve the final linear polynomial
 515                 (setf (aref *pr-sl* 1) (- (/ (aref *pr-sl* 1) (aref *pr-sl* 0))))))
 516          (setq *nn* 0))
 517       ;; Calculate a lower bound on the modulus of the zeros.
 518       (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i))))
 519       (setq bnd (cauchy-sl))
 520
 521       ;; Compute derivative of polynomial.  Save result in *hr-sl*.
 522       (do ((i 1 (1+ i)))
 523           ((> i *n*))
 524         (setf (aref *hr-sl* i) (/ (* (float (- *n* i)) (aref *pr-sl* i)) (float *n*))))
 525       (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
 526       (setq aa (aref *pr-sl* *nn*) bb (aref *pr-sl* *n*) zerok (zerop (aref *hr-sl* *n*)))
 527       ;; Do 5 steps with no shift.
 528       (do ((jj 1 (1+ jj)))
 529           ((> jj 5.))
 530         (setq cc (aref *hr-sl* *n*))
 531         (cond (zerok (do ((j *n* (1- j)))
 532                          ((< j 1))
 533                        (setf (aref *hr-sl* j) (aref *hr-sl* (1- j))))
 534                      (setf (aref *hr-sl* 0) 0.0)
 535                      (setq zerok (zerop (aref *hr-sl* *n*))))
 536               (t (setq t1 (- (/ aa cc)))
 537                  (do ((j *n* (1- j)))
 538                      ((< j 1))
 539                    (setf (aref *hr-sl* j) (+ (* t1 (aref *hr-sl* (1- j))) (aref *pr-sl* j))))
 540                  (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
 541                  (setq zerok (not (> (abs (aref *hr-sl* *n*))
 542                                      (* (abs bb) *are* 10.0)))))))
 543       (do ((i 0 (1+ i))) ((> i *n*)) (setf (aref *shi-sl* i) (aref *hr-sl* i)))
 544       (do ((cnt 1 (1+ cnt)))
 545           ((> cnt 20.) (setq conv1 nil))
 546         (setq xx (prog1
 547                      (- (* cosr xx) (* sinr yy))
 548                    (setq yy (+ (* sinr xx) (* cosr yy))))
 549               *sr* (* bnd xx)
 550               *u* (* -2.0 *sr*)
 551               *v* bnd)
 552         (fxshfr-sl (* 20 cnt))
 553         (cond ((> *nz* 0)
 554                (setf (aref *pr-sl* *nn*) *szr*)
 555                (setf (aref *pi-sl* *nn*) *szi*)
 556                (cond ((= *nz* 2)
 557                       (setf (aref *pr-sl* *n*) *lzr*)
 558                       (setf (aref *pi-sl* *n*) *lzi*)
 559                       (cond ((and $polyfactor (not (zerop *szi*)))
 560                              (setf (aref *pr-sl* *nn*) *v*)
 561                              (setf (aref *pr-sl* *n*) *u*)
 562                              (setf (aref *pi-sl* *nn*) 1.0)))))
 563                (setq *nn* (- *nn* *nz*) *n* (1- *nn*))
 564                (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *pr-sl* i) (aref *qpr-sl* i)))
 565                (return nil)))
 566         (do ((i 0 (1+ i))) ((> i *n*)) (setf (aref *hr-sl* i) (aref *shi-sl* i))))
 567       (or conv1 (return nil)))
 568     (cond ($polyfactor
 569            (do ((i degree (1- i)))
 570                ((= i *nn*))
 571              (cond ((zerop (aref *pi-sl* i))
 572                     (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*)))
 573                    (t (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) (* 2 *polysc1*)))
 574                       (setq i (1- i))
 575                       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*))))))
 576           (t (do ((i (1+ *nn*) (1+ i)))
 577                  ((> i degree))
 578                (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*))
 579                (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) *polysc1*)))))
 580     (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
 581         ((> i *nn*))
 582       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j)))))
 583
 584 (defun fxshfr-sl (l2)
 585   (let ((*my-type* 0)
 586         (*a* 0.0) (*b* 0.0) (*c* 0.0) (*d* 0.0) (*e* 0.0) (*f* 0.0) (*g* 0.0) (*h* 0.0)
 587         (*a1* 0.0) (*a3* 0.0) (*a7* 0.0))
 588     (declare (special *my-type*))
 589     (setq *nz* 0)
 590     (quadsd-sl)
 591     (calcsc-sl)
 592     (do ((j 1 (1+ j))
 593          (betav 0.25)
 594          (betas 0.25)
 595          (oss *sr*)
 596          (ovv *v*)
 597          (tvv) (tss) (ss) (vv) (tv) (ts) (ots) (otv)
 598          (*ui*) (*vi*) (*s*) (svv) (svu) (iflag) (vpass) (spass) (vtry) (stry))
 599         ((> j l2))
 600       (nextk-sl)
 601       (calcsc-sl)
 602       (newest-sl)
 603       (setq vv *vi*
 604             ss 0.0)
 605       (or (zerop (aref *hr-sl* *n*))
 606           (setq ss (- (/ (aref *pr-sl* *nn*) (aref *hr-sl* *n*)))))
 607       (setq tv 1.0 ts 1.0)
 608       (cond ((not (or (= j 1) (= *my-type* 3)))
 609              (or (zerop vv) (setq tv (abs (/ (- vv ovv) vv))))
 610              (or (zerop ss) (setq ts (abs (/ (- ss oss) ss))))
 611              (setq tvv 1.0)
 612              (and (< tv otv) (setq tvv (* tv otv)))
 613              (setq tss 1.0)
 614              (and (< ts ots) (setq tss (* ts ots)))
 615              (setq vpass (< tvv betav) spass (< tss betas))
 616              (cond ((or spass vpass)
 617                     (setq svu *u* svv *v*)
 618                     (do ((i 0 (1+ i)))
 619                         ((> i *n*)) (setf (aref *shr-sl* i)
 620                                           (aref *hr-sl* i)))
 621                     (setq *s* ss vtry nil stry nil)
 622                     (and (do ((*bool* (not (and spass (or (not vpass) (< tss tvv)))) t)
 623                               (l50 nil nil))
 624                              (nil)
 625                            (cond (*bool* (quadit-sl)
 626                                          (and (> *nz* 0) (return t))
 627                                          (setq vtry t
 628                                                betav (* 0.25 betav))
 629                                          (cond ((or stry (not spass))
 630                                                 (setq l50 t))
 631                                                (t (do ((i 0 (1+ i)))
 632                                                       ((> i *n*))
 633                                                     (setf (aref *hr-sl* i)
 634                                                           (aref *shr-sl* i)))))))
 635                            (cond ((not l50)
 636                                   (setq iflag (realit-sl))
 637                                   (and (> *nz* 0) (return t))
 638                                   (setq stry t betas (* 0.25 betas))
 639                                   (cond ((zerop iflag) (setq l50 t))
 640                                         (t (setq *ui* (- (+ *s* *s*))
 641                                                  *vi* (* *s* *s*))))))
 642                            (cond (l50 (setq *u* svu *v* svv)
 643                                       (do ((i 0 (1+ i)))
 644                                           ((> i *n*))
 645                                         (setf (aref *hr-sl* i)
 646                                               (aref *shr-sl* i)))
 647                                       (and (or (not vpass) vtry)
 648                                            (return nil)))))
 649                          (return nil))
 650                     (quadsd-sl)
 651                     (calcsc-sl)))))
 652       (setq ovv vv
 653             oss ss
 654             otv tv
 655             ots ts))))
 656
 657 (defun quadit-sl nil
 658   (setq *nz* 0 *u* *ui* *v* *vi*)
 659   (do ((tried) (j 0) (ee) (mp) (relstp) (omp) (ms))
 660       (nil)
 661     (quad-sl 1.0 *u* *v*)
 662     (and (> (abs (- (abs *szr*) (abs *lzr*))) (* 1e-2 (abs *lzr*)))
 663          (return nil))
 664     (quadsd-sl)
 665     (setf mp (+ (abs (- *a* (* *szr* *b*))) (abs (* *szi* *b*))))
 666     (setf ms (cmod-sl *szr* *szi*))
 667     (setf ee (errev-sl ms mp))
 668     (cond ((not (> mp (* 2e1 ee))) (setq *nz* 2)
 669            (return nil)))
 670     (setq j (1+ j))
 671     (and (> j 20) (return nil))
 672     (cond ((not (or (< j 2) (> relstp 1e-2) (< mp omp) tried))
 673            (and (< relstp *are*) (setq relstp *are*))
 674            (setq relstp (sqrt relstp)
 675                  *u* (- *u* (* *u* relstp))
 676                  *v* (+ *v* (* *v* relstp)))
 677            (quadsd-sl)
 678            (do ((i 1 (1+ i)))
 679                ((> i 5)) (calcsc-sl) (nextk-sl))
 680            (setq tried t j 0)))
 681     (setq omp mp)
 682     (calcsc-sl)
 683     (nextk-sl)
 684     (calcsc-sl)
 685     (newest-sl)
 686     (and (zerop *vi*) (return nil))
 687     (setq relstp (abs (/ (- *vi* *v*) *vi*)) *u* *ui* *v* *vi*)))
 688
 689 (defun realit-sl ()
 690   (setq *nz* 0)
 691   (do ((j 0) (pv) (ee) (ms) (mp) (kv) (t1) (omp))
 692       (nil)
 693     (setq pv (aref *pr-sl* 0))
 694     (setf (aref *qpr-sl* 0) pv)
 695     (do ((i 1 (1+ i)))
 696         ((> i *nn*))
 697       (setq pv (+ (* pv *s*) (aref *pr-sl* i)))
 698       (setf (aref *qpr-sl* i) pv))
 699     (setq mp (abs pv)
 700           ms (abs *s*)
 701           ee (* (/ *mre* (+ *are* *mre*)) (abs (aref *qpr-sl* 0))))
 702     (do ((i 1 (1+ i)))
 703         ((> i *nn*)) (setq ee (+ (* ee ms) (abs (aref *qpr-sl* i)))))
 704     (cond ((not (> mp (* 2e1 (- (* (+ *are* *mre*) ee) (* *mre* mp)))))
 705            (setq *nz* 1 *szr* *s* *szi* 0.0)
 706            (return 0)))
 707     (setq j (1+ j))
 708     (and (> j 10) (return 0))
 709     (cond ((not (or (< j 2)
 710                     (> (abs t1) (* 1e-3 (abs (- *s* t1))))
 711                     (not (> mp omp))))
 712            (return 1)))
 713     (setq omp mp kv (aref *hr-sl* 0))
 714     (setf (aref *qhr-sl* 0) kv)
 715     (do ((i 1 (1+ i)))
 716         ((> i *n*))
 717       (setq kv (+ (* kv *s*) (aref *hr-sl* i)))
 718       (setf (aref *qhr-sl* i) kv))
 719     (cond ((> (abs kv) (* (abs (aref *hr-sl* *n*)) 1e1 *are*))
 720            (setq t1 (- (/ pv kv)))
 721            (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
 722            (do ((i 1 (1+ i)))
 723                ((> i *n*))
 724              (setf (aref *hr-sl* i)
 725                     (+ (* t1 (aref *qhr-sl* (1- i))) (aref *qpr-sl* i)))))
 726           (t (setf (aref *hr-sl* 0) 0.0)
 727              (do ((i 1 (1+ i)))
 728                  ((> i *n*)) (setf (aref *hr-sl* i) (aref *qhr-sl* (1- i))))))
 729     (setq kv (aref *hr-sl* 0))
 730     (do ((i 1 (1+ i)))
 731         ((> i *n*)) (setq kv (+ (* kv *s*) (aref *hr-sl* i))))
 732     (setq t1 0.0)
 733     (and (> (abs kv) (* (abs (aref *hr-sl* *n*)) 10.0 *are*))
 734          (setq t1 (- (/ pv kv))))
 735     (setq *s* (+ *s* t1))))
 736
 737 (defun calcsc-sl ()
 738   (declare (special *my-type*))
 739   (setq *d* (aref *hr-sl* 0))
 740   (setf (aref *qhr-sl* 0) *d*)
 741   (setq *c* (- (aref *hr-sl* 1) (* *u* *d*)))
 742   (setf (aref *qhr-sl* 1) *c*)
 743   (do ((i 2 (1+ i))
 744        (c0))
 745       ((> i *n*))
 746     (setq c0 (- (aref *hr-sl* i) (* *u* *c*) (* *v* *d*)))
 747     (setf (aref *qhr-sl* i) c0)
 748     (setq *d* *c* *c* c0))
 749   (cond ((not (or (> (abs *c*) (* (abs (aref *hr-sl* *n*)) 1e2 *are*))
 750                   (> (abs *d*) (* (abs (aref *hr-sl* (1- *n*))) 1e2 *are*))))
 751          (setq *my-type* 3))
 752         ((not (< (abs *d*) (abs *c*)))
 753          (setq *my-type* 2
 754                *e* (/ *a* *d*)
 755                *f* (/ *c* *d*)
 756                *g* (* *u* *b*)
 757                *h* (* *v* *b*)
 758                *a3* (+ (* (+ *a* *g*) *e*) (* *h* (/ *b* *d*)))
 759                *a1* (- (* *b* *f*) *a*)
 760                *a7* (+ (* (+ *f* *u*) *a*) *h*)))
 761         (t (setq *my-type* 1
 762                  *e* (/ *a* *c*)
 763                  *f* (/ *d* *c*)
 764                  *g* (* *u* *e*)
 765                  *h* (* *v* *b*)
 766                  *a3* (+ (* *a* *e*) (* (+ (/ *h* *c*) *g*) *b*))
 767                  *a1* (- *b* (* *a* (/ *d* *c*)))
 768                  *a7* (+ *a* (* *g* *d*) (* *h* *f*)))))
 769   nil)
 770
 771 (defun nextk-sl ()
 772   (declare (special *my-type*))
 773   (cond ((= *my-type* 3)
 774          (setf (aref *hr-sl* 0) 0.0)
 775          (setf (aref *hr-sl* 1) 0.0)
 776          (do ((i 2 (1+ i)))
 777              ((> i *n*)) (setf (aref *hr-sl* i) (aref *qhr-sl* (- i 2)))))
 778         ((> (abs *a1*) (* (abs (if (= *my-type* 1) *b* *a*)) 1e1 *are*))
 779          (setq *a7* (/ *a7* *a1*) *a3* (/ *a3* *a1*))
 780          (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
 781          (setf (aref *hr-sl* 1) (- (aref *qpr-sl* 1) (* *a7* (aref *qpr-sl* 0))))
 782          (do ((i 2 (1+ i)))
 783              ((> i *n*))
 784            (setf (aref *hr-sl* i)
 785                   (+ (* *a3* (aref *qhr-sl* (- i 2)))
 786                       (- (* *a7* (aref *qpr-sl* (1- i))))
 787                       (aref *qpr-sl* i)))))
 788         (t (setf (aref *hr-sl* 0) 0.0)
 789            (setf (aref *hr-sl* 1) (- (* *a7* (aref *qpr-sl* 0))))
 790            (do ((i 2 (1+ i)))
 791                ((> i *n*))
 792              (setf (aref *hr-sl* i)
 793                     (- (* *a3* (aref *qhr-sl* (- i 2)))
 794                         (* *a7* (aref *qpr-sl* (1- i))))))))
 795   nil)
 796
 797 (defun newest-sl ()
 798   (declare (special *my-type*))
 799   (let ((a4 0.0) (a5 0.0)
 800         (b1 0.0) (b2 0.0)
 801         (c1 0.0) (c2 0.0) (c3 0.0) (c4 0.0))
 802     (cond ((= *my-type* 3)
 803            (setq *ui* 0.0 *vi* 0.0))
 804           (t
 805            (if (= *my-type* 2)
 806                (setq a4 (+ (* (+ *a* *g*) *f*) *h*)
 807                      a5 (+ (* (+ *f* *u*) *c*) (* *v* *d*)))
 808                (setq a4 (+ *a* (* *u* *b*) (* *h* *f*))
 809                      a5 (+ *c* (* (+ *u* (* *v* *f*)) *d*))))
 810              (setq b1 (- (/ (aref *hr-sl* *n*) (aref *pr-sl* *nn*)))
 811                    b2 (- (/ (+ (aref *hr-sl* (1- *n*)) (* b1 (aref *pr-sl* *n*))) (aref *pr-sl* *nn*)))
 812                    c1 (* *v* b2 *a1*)
 813                    c2 (* b1 *a7*)
 814                    c3 (* b1 b1 *a3*)
 815                    c4 (- c1 c2 c3)
 816                    c1 (+ a5 (* b1 a4) (- c4)))
 817              (if (zerop c1)
 818                  (setq *ui* 0.0 *vi* 0.0)
 819                  (setq *ui* (- *u* (/ (+ (* *u* (+ c3 c2))
 820                                          (* *v* (+ (* b1 *a1*) (* b2 *a7*))))
 821                                       c1))
 822                        *vi* (* *v* (1+ (/ c4 c1)))))))
 823     nil))
 824
 825 ;; Divide the polynomial in *pr-sl* by the quadratic x^2 + (*u*)*x +
 826 ;; (*v*).  Place the quotient in *qpr-sl* and the remainder in *a* and
 827 ;; *b*.  I (rtoy) think the remainder polynomial is (*b*)*x + (*a*).
 828 (defun quadsd-sl ()
 829   (setq *b* (aref *pr-sl* 0))
 830   (setf (aref *qpr-sl* 0) *b*)
 831   (setq *a* (- (aref *pr-sl* 1) (* *u* *b*)))
 832   (setf (aref *qpr-sl* 1) *a*)
 833   (do ((i 2 (1+ i))
 834        (c0))
 835       ((> i *nn*))
 836     (setq c0 (- (aref *pr-sl* i) (* *u* *a*) (* *v* *b*)))
 837     (setf (aref *qpr-sl* i) c0)
 838     (setq *b* *a*
 839           *a* c0)))
 840
 841 ;; Compute the zeros of the quadratic a0*z^2+b1*z+c0.  The larger zero
 842 ;; is returned in *szr* and *szi*.  The smaller zero is in *lzr* and
 843 ;; *lzi*.
 844 (defun quad-sl (a0 b1 c0)
 845   (setq *szr* 0.0 *szi* 0.0 *lzr* 0.0 *lzi* 0.0)
 846   (let ((b0 0.0)
 847         (l0 0.0)
 848         (*e* 0.0))
 849     ;; Handle the degenerate cases of a0 = 0 or c0 = 0 first.
 850     (cond ((zerop a0) (or (zerop b1) (setq *szr* (- (/ c0 b1)))))
 851           ((zerop c0) (setq *lzr* (- (/ b1 a0))))
 852           (t
 853            ;; Quadratic formula.
 854            (setq b0 (/ b1 2.0))
 855            (cond ((< (abs b0) (abs c0))
 856                   (setq *e* a0)
 857                   (and (< c0 0.0) (setq *e* (- a0)))
 858                   (setq *e* (- (* b0 (/ b0 (abs c0))) *e*)
 859                         l0 (* (sqrt (abs *e*)) (sqrt (abs c0)))))
 860                  (t (setq *e* (- 1.0 (* (/ a0 b0) (/ c0 b0)))
 861                           l0 (* (sqrt (abs *e*)) (abs b0)))))
 862            (cond ((< *e* 0.0)
 863                   (setq *szr* (- (/ b0 a0))
 864                         *lzr* *szr*
 865                         *szi* (abs (/ l0 a0))
 866                         *lzi* (- *szi*)))
 867                  (t (or (< b0 0.0) (setq l0 (- l0)))
 868                     (setq *lzr* (/ (- l0 b0) a0))
 869                     (or (zerop *lzr*) (setq *szr* (/ c0 *lzr* a0)))))))
 870     nil))
 871 \f
 872 ;; This is a very straightforward conversion of $allroots to use
 873 ;; bfloats instead of floats.
 874
 875 (defun bf-cpoly-err (expr)
 876   (merror (intl:gettext "bfallroots: expected a polynomial; found ~M") expr))
 877
 878 (defun fpzerop (x)
 879   (equal '(0 0) x))
 880
 881 ;; (ar+%i*ai)/(br+%i*bi) -> cr+%i*ci.
 882 (defun bf-cdivid-sl (ar ai br bi)
 883   (cond ((and (fpzerop br)
 884               (fpzerop bi))
 885          ;; Division by zero.  Should we do something else besides set
 886          ;; both parts to be "infinity"?
 887          (setq *cr* (setq *ci* *infin*)))
 888         ((fpgreaterp (fpabs bi) (fpabs br))
 889          (let ((r1 (fpquotient br bi)))
 890            (setq bi (fpplus bi (fptimes* br r1))
 891                  br (fpplus ai (fptimes* ar r1))
 892                  *cr* (fpquotient br bi)
 893                  br (fpdifference (fptimes* ai r1) ar)
 894                  *ci* (fpquotient br bi))))
 895         (t
 896          (let ((r1 (fpquotient bi br)))
 897            (setq bi (fpplus br (fptimes* bi r1))
 898                  br (fpplus ar (fptimes* ai r1))
 899                  *cr* (fpquotient br bi)
 900                  br (fpdifference ai (fptimes* ar r1))
 901                  *ci* (fpquotient br bi))))))
 902
 903 (defun fpsqrt (x)
 904   (fproot (bcons x) 2))
 905
 906 (defun bf-cmod-sl (ar ai)
 907   (let ((ar (fpabs ar))
 908         (ai (fpabs ai)))
 909     (cond ((fpgreaterp ai ar)
 910            (setq ar (fpquotient ar ai))
 911            (fptimes* ai
 912                      (fpsqrt (fpplus (fpone) (fptimes* ar ar)))))
 913           ((fpgreaterp ar ai)
 914            (setq ai (fpquotient ai ar))
 915            (fptimes* ar (fpsqrt (fpplus (fpone) (fptimes* ai ai)))))
 916           ((fptimes* (fpsqrt (intofp 2))
 917                      ar)))))
 918
 919
 920 (defun bf-calct-sl nil
 921   (do ((i 1 (1+ i))
 922        (tt)
 923        (hvr (setf (aref *qhr-sl* 0) (aref *hr-sl* 0)))
 924        (hvi (setf (aref *qhi-sl* 0) (aref *hi-sl* 0))))
 925       ((> i *n*)
 926        (setq *bool*
 927              (not (fpgreaterp (bf-cmod-sl hvr hvi)
 928                               (fptimes* (intofp 10)
 929                                         (fptimes* *are*
 930                                                   (bf-cmod-sl (aref *hr-sl* *n*)
 931                                                               (aref *hi-sl* *n*)))))))
 932        (cond ((not *bool*)
 933               (bf-cdivid-sl (fpminus *pvr*) (fpminus *pvi*) hvr hvi)
 934               (setq *tr* *cr*
 935                     *ti* *ci*))
 936              (t (setq *tr* (intofp 0) *ti* (intofp 0))))
 937        nil)
 938     (setq tt (fpdifference (fpplus (aref *hr-sl* i)
 939                                    (fptimes* hvr *sr*))
 940                            (fptimes* hvi *si*)))
 941     (setf (aref *qhi-sl* i)
 942           (setq hvi (fpplus (aref *hi-sl* i)
 943                             (fpplus (fptimes* hvr *si*)
 944                                     (fptimes* hvi *sr*)))))
 945     (setf (aref *qhr-sl* i) (setq hvr tt))))
 946
 947 (defun bf-nexth-sl ()
 948   (cond (*bool*
 949          (do ((j 1 (1+ j)))
 950              ((> j *n*))
 951            (setf (aref *hr-sl* j) (aref *qhr-sl* (1- j)))
 952            (setf (aref *hi-sl* j) (aref *qhi-sl* (1- j))))
 953          (setf (aref *hr-sl* 0) (intofp 0))
 954          (setf (aref *hi-sl* 0) (intofp 0)))
 955         (t
 956          (do ((j 1. (1+ j))
 957               (t1)
 958               (t2))
 959              ((> j *n*))
 960            (setq t1 (aref *qhr-sl* (1- j))
 961                  t2 (aref *qhi-sl* (1- j)))
 962            (setf (aref *hr-sl* j)
 963                  (fpdifference (fpplus (aref *qpr-sl* j)
 964                                        (fptimes* t1 *tr*))
 965                                (fptimes* t2 *ti*)))
 966            (setf (aref *hi-sl* j)
 967                  (fpplus (aref *qpi-sl* j)
 968                          (fpplus (fptimes* t1 *ti*)
 969                                  (fptimes* t2 *tr*)))))
 970          (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
 971          (setf (aref *hi-sl* 0) (aref *qpi-sl* 0))))
 972   nil)
 973
 974 (defun bf-polyev-sl ()
 975   (setq *pvr* (setf (aref *qpr-sl* 0) (aref *pr-sl* 0))
 976         *pvi* (setf (aref *qpi-sl* 0) (aref *pi-sl* 0)))
 977   (do ((i 1 (1+ i))
 978        (tt))
 979       ((> i *nn*))
 980     (setq tt (fpdifference (fpplus (aref *pr-sl* i) (fptimes* *pvr* *sr*))
 981                            (fptimes* *pvi* *si*)))
 982     (setf (aref *qpi-sl* i)
 983           (setq *pvi* (fpplus (aref *pi-sl* i)
 984                               (fpplus (fptimes* *pvr* *si*)
 985                                       (fptimes* *pvi* *sr*)))))
 986     (setf (aref *qpr-sl* i)
 987           (setq *pvr* tt))))
 988
 989 (defun bf-errev-sl (ms mp)
 990   (fpdifference
 991    (fptimes* (do ((j 0 (1+ j))
 992                   (e (fpquotient (fptimes* (bf-cmod-sl (aref *qpr-sl* 0) (aref *qpi-sl* 0))
 993                                            *mre*)
 994                                  (fpplus *are* *mre*))))
 995                  ((> j *nn*) e)
 996                (setq e (fpplus (bf-cmod-sl (aref *qpr-sl* j) (aref *qpi-sl* j))
 997                                (fptimes* e ms))))
 998              (fpplus *are* *mre*))
 999    (fptimes* mp *mre*)))
1000
1001 (defun bf-cauchy-sl ()
1002   (let ((x (fpexp (fpquotient (fpdifference (fplog (aref *shr-sl* *nn*))
1003                                             (fplog (aref *shr-sl* 0)))
1004                               (intofp *nn*))))
1005         (xm (intofp 0)))
1006     (setf (aref *shr-sl* *nn*) (fpminus (aref *shr-sl* *nn*)))
1007     (cond ((not (fpzerop (aref *shr-sl* *n*)))
1008            (setq xm (fpminus (fpquotient (aref *shr-sl* *nn*)
1009                                          (aref *shr-sl* *n*))))
1010            (cond ((fpgreaterp x xm) (setq x xm)))))
1011     (do ((f))
1012         (nil)
1013       (setq xm (fptimes* (intofp 0.1) x)
1014             f (aref *shr-sl* 0))
1015       (do ((i 1 (1+ i)))
1016           ((> i *nn*))
1017         (setq f (fpplus (aref *shr-sl* i)
1018                         (fptimes* f xm))))
1019       ;;(cond ((not (< 0.0 f)) (return t)))
1020       (when (fpgreaterp (intofp 0) f)
1021         (return t))
1022       (setq x xm))
1023     (do ((dx x)
1024          (df)
1025          (f))
1026         ((fpgreaterp (intofp 5e-3)
1027                      (fpabs (fpquotient dx x)))
1028          x)
1029       (setq f (aref *shr-sl* 0)
1030             df f)
1031       (do ((i 1 (1+ i)))
1032           ((> i *n*))
1033         (setq f (fpplus (fptimes* f x)
1034                         (aref *shr-sl* i))
1035               df (fpplus (fptimes* df x)
1036                          f)))
1037       (setq f (fpplus (fptimes* f x)
1038                       (aref *shr-sl* *nn*))
1039             dx (fpquotient f df)
1040             x (fpdifference x dx)))))
1041
1042 (defun bf-scale-float (bf scale)
1043   (destructuring-bind (mantissa exp)
1044       bf
1045     (if (zerop mantissa)
1046         (list mantissa 0)
1047         (list mantissa
1048               (+ exp scale)))))
1049
1050 (defun bf-scale-sl ()
1051   (do ((i 0 (1+ i))
1052        (j 0)
1053        (x (intofp 0))
1054        (dx (intofp 0)))
1055       ((> i *nn*)
1056        (setq x (fpquotient x (intofp (- (1+ *nn*) j)))
1057              dx (fpquotient (fpdifference (fplog (aref *shr-sl* *nn*))
1058                                           (fplog (aref *shr-sl* 0)))
1059                             (intofp *nn*))
1060              *polysc1* (fpentier (bcons (fpplus (cdr bfhalf)
1061                                                 (fpquotient dx *logbas*))))
1062              x (fpplus x (fptimes* (intofp (* *polysc1* *nn*))
1063                                    (fptimes* *logbas*
1064                                              (cdr bfhalf))))
1065              *polysc* (fpentier (bcons (fpplus (cdr bfhalf) (fpquotient x *logbas*))))))
1066     (cond ((equalp (aref *shr-sl* i) (intofp 0))
1067            (setq j (1+ j)))
1068           (t
1069            (setq x (fpplus x (fplog (aref *shr-sl* i)))))))
1070   (do ((i *nn* (1- i))
1071        (j (- *polysc*) (+ j *polysc1*)))
1072       ((< i 0))
1073     (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) j))
1074     (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) j)))
1075   nil)
1076
1077 (defun bf-noshft-sl (l1)
1078   (do ((i 0 (1+ i))
1079        (xni (intofp *nn*) (fpdifference xni (intofp 1)))
1080        (t1 (fpquotient (fpone) (intofp *nn*))))
1081       ((> i *n*))
1082     (setf (aref *hr-sl* i) (fptimes* (aref *pr-sl* i)
1083                                      (fptimes* xni t1)))
1084     (setf (aref *hi-sl* i) (fptimes* (aref *pi-sl* i)
1085                                      (fptimes* xni t1))))
1086   (do ((jj 1 (1+ jj)))
1087       ((> jj l1))
1088     (cond ((fpgreaterp (bf-cmod-sl (aref *hr-sl* *n*) (aref *hi-sl* *n*))
1089                        (fptimes* (intofp 10)
1090                                  (fptimes* *are*
1091                                            (bf-cmod-sl (aref *pr-sl* *n*)
1092                                                        (aref *pi-sl* *n*)))))
1093            (bf-cdivid-sl (fpminus (aref *pr-sl* *nn*))
1094                          (fpminus (aref *pi-sl* *nn*))
1095                          (aref *hr-sl* *n*)
1096                          (aref *hi-sl* *n*))
1097            (setq *tr* *cr*
1098                  *ti* *ci*)
1099            (do ((j *n* (1- j)) (t1) (t2))
1100                ((> 1 j))
1101              (setq t1 (aref *hr-sl* (1- j))
1102                    t2 (aref *hi-sl* (1- j)))
1103              (setf (aref *hr-sl* j) (fpdifference (fpplus (aref *pr-sl* j)
1104                                                           (fptimes* t1 *tr*))
1105                                                   (fptimes* t2 *ti*)))
1106              (setf (aref *hi-sl* j) (fpplus (aref *pi-sl* j)
1107                                             (fpplus (fptimes* t1 *ti*)
1108                                                     (fptimes* t2 *tr*)))))
1109            (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
1110            (setf (aref *hi-sl* 0) (aref *pi-sl* 0)))
1111           (t (do ((j *n* (1- j)))
1112                  ((> 1 j))
1113                (setf (aref *hr-sl* j) (aref *hr-sl* (1- j)))
1114                (setf (aref *hi-sl* j) (aref *hi-sl* (1- j))))
1115              (setf (aref *hr-sl* 0) (intofp 0))
1116              (setf (aref *hi-sl* 0) (intofp 0))))))
1117
1118 (defun bf-vrshft-sl (l3)
1119   (setq *conv* nil
1120         *sr* *zr*
1121         *si* *zi*)
1122   (do ((i 1 (1+ i))
1123        (bool1 nil)
1124        (mp)
1125        (ms)
1126        (omp)
1127        (relstp)
1128        (tp)
1129        (r1))
1130       ((> i l3))
1131     (bf-polyev-sl)
1132     (setq mp (bf-cmod-sl *pvr* *pvi*)
1133           ms (bf-cmod-sl *sr* *si*))
1134     (cond ((fpgreaterp (fptimes* (intofp 20) (bf-errev-sl ms mp))
1135                        mp)
1136            (setq *conv* t
1137                  *zr* *sr*
1138                  *zi* *si*)
1139            (return t)))
1140     (cond ((= i 1)
1141            (setq omp mp))
1142           ((or bool1
1143                (fpgreaterp omp mp)
1144                ;;(not (< relstp 0.05))
1145                (fpgreaterp relstp (cdr bfhalf)))
1146            (if (fpgreaterp (fptimes* (intofp 0.1) mp)
1147                            omp)
1148                (return t)
1149                (setq omp mp)))
1150           (t
1151            (setq tp relstp
1152                  bool1 t)
1153            (when (fpgreaterp *are* relstp)
1154              (setq tp *are*))
1155            (setq r1 (fpsqrt tp)
1156                  *sr* (prog1
1157                           (fpdifference (fptimes* (fpplus (fpone) r1)
1158                                                   *sr*)
1159                                         (fptimes* r1 *si*))
1160                         (setq *si* (fpplus (fptimes* (fpplus (fpone) r1)
1161                                                      *si*)
1162                                            (fptimes* r1 *sr*)))))
1163            (bf-polyev-sl)
1164            (do ((j 1 (1+ j)))
1165                ((> j 5))
1166              (bf-calct-sl)
1167              (bf-nexth-sl))
1168            (setq omp *infin*)))
1169     (bf-calct-sl)
1170     (bf-nexth-sl)
1171     (bf-calct-sl)
1172     (or *bool*
1173         (setq relstp (fpquotient (bf-cmod-sl *tr* *ti*)
1174                                  (bf-cmod-sl *sr* *si*))
1175               *sr* (fpplus *sr* *tr*)
1176               *si* (fpplus *si* *ti*)))))
1177
1178 (defun bf-fxshft-sl (l2)
1179   (let ((test t)
1180         (pasd nil)
1181         (otr (intofp 0))
1182         (oti (intofp 0))
1183         (svsr (intofp 0))
1184         (svsi (intofp 0))
1185         (*bool* nil)
1186         (*pvr* (intofp 0))
1187         (*pvi* (intofp 0)))
1188     (bf-polyev-sl)
1189     (setq *conv* nil)
1190     (bf-calct-sl)
1191     (do ((j 1 (1+ j)))
1192         ((> j l2))
1193       (setq otr *tr*
1194             oti *ti*)
1195       (bf-nexth-sl)
1196       (bf-calct-sl)
1197       (setq *zr* (fpplus *sr* *tr*)
1198             *zi* (fpplus *si* *ti*))
1199       (cond ((and (not *bool*)
1200                   test
1201                   (not (= j l2)))
1202              (cond ((fpgreaterp (fptimes* (cdr bfhalf) (bf-cmod-sl *zr* *zi*))
1203                                 (bf-cmod-sl (fpdifference *tr* otr)
1204                                             (fpdifference *ti* oti)))
1205                     (cond (pasd
1206                            (do ((i 0 (1+ i)))
1207                                ((> i *n*))
1208                              (setf (aref *shr-sl* i) (aref *hr-sl* i))
1209                              (setf (aref *shi-sl* i) (aref *hi-sl* i)))
1210                            (setq svsr *sr* svsi *si*)
1211                            (bf-vrshft-sl 10.)
1212                            (when *conv* (return nil))
1213                            (setq test nil)
1214                            (do ((i 0 (1+ i)))
1215                                ((> i *n*))
1216                              (setf (aref *hr-sl* i) (aref *shr-sl* i))
1217                              (setf (aref *hi-sl* i) (aref *shi-sl* i)))
1218                            (setq *sr* svsr *si* svsi)
1219                            (bf-polyev-sl)
1220                            (bf-calct-sl))
1221                           ((setq pasd t))))
1222                    ((setq pasd nil))))))
1223     (or *conv* (bf-vrshft-sl 10))
1224     nil))
1225
1226 (defun bf-cpoly-sl (degree)
1227   (let ( ;; Log of our floating point base.  (Do we need this much accuracy?)
1228         (*logbas* (fplog (intofp 2)))
1229         ;; "Largest" bfloat.  What should we use?
1230         (*infin* (intofp +most-positive-flonum+))
1231         ;; bfloat epsilon.  2^(-fpprec)
1232         (*are* (bf-scale-float (intofp 2) (- fpprec)))
1233         (*mre* (intofp 0))
1234         (xx (fproot bfhalf 2))
1235         (yy (intofp 0))
1236         ;; cos(94deg).  Probably don't need full bfloat precision here.
1237         (cosr (intofp -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804))
1238         ;; sin(94deg).  Probably don't need full bfloat precision here.
1239         (sinr (intofp 0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315))
1240         (*cr* (intofp 0))
1241         (*ci* (intofp 0))
1242         (*sr* (intofp 0))
1243         (*si* (intofp 0))
1244         (*tr* (intofp 0))
1245         (*ti* (intofp 0))
1246         (*zr* (intofp 0))
1247         (*zi* (intofp 0))
1248         (bnd (intofp 0))
1249         (*n* 0)
1250         (*polysc* 0)
1251         (*polysc1* 0)
1252         (*conv* nil))
1253     (setq *mre* (fptimes* (intofp 2)
1254                           (fptimes* (fpsqrt (intofp 2)) *are*))
1255           yy (fpminus xx))
1256     (do ((i degree (1- i)))
1257         ((not (and (equalp (intofp 0) (aref *pr-sl* i))
1258                    (equalp (intofp 0) (aref *pi-sl* i))))
1259          (setq *nn* i
1260                *n* (1- i))))
1261     (setq degree *nn*)
1262     (do ((i 0 (1+ i)))
1263         ((> i *nn*))
1264       (setf (aref *shr-sl* i) (bf-cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
1265     (if (> *nn* 0) (bf-scale-sl))
1266     (do ()
1267         ((> 2 *nn*)
1268          (bf-cdivid-sl (fpminus (aref *pr-sl* 1))
1269                        (fpminus (aref *pi-sl* 1))
1270                        (aref *pr-sl* 0)
1271                        (aref *pi-sl* 0))
1272          (setf (aref *pr-sl* 1) *cr*)
1273          (setf (aref *pi-sl* 1) *ci*)
1274          (setq *nn* 0))
1275       (do ((i 0 (1+ i)))
1276           ((> i *nn*))
1277         (setf (aref *shr-sl* i) (bf-cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
1278       (setq bnd (bf-cauchy-sl))
1279       (catch 'newroot
1280         (do ((cnt1 1 (1+ cnt1)))
1281             ((> cnt1 2))
1282           (bf-noshft-sl 5)
1283           (do ((cnt2 1 (1+ cnt2)))
1284               ((> cnt2 9))
1285             (setq xx (prog1
1286                          (fpdifference (fptimes* cosr xx)
1287                                        (fptimes* sinr yy))
1288                        (setq yy (fpplus (fptimes* sinr xx)
1289                                         (fptimes* cosr yy))))
1290                   *sr* (fptimes* bnd xx)
1291                   *si* (fptimes* bnd yy))
1292             (bf-fxshft-sl (* 10 cnt2))
1293             (cond (*conv* (setf (aref *pr-sl* *nn*) *zr*)
1294                           (setf (aref *pi-sl* *nn*) *zi*)
1295                           (setq *nn* *n* *n* (1- *n*))
1296                           (do ((i 0 (1+ i)))
1297                               ((> i *nn*))
1298                             (setf (aref *pr-sl* i) (aref *qpr-sl* i))
1299                             (setf (aref *pi-sl* i) (aref *qpi-sl* i)))
1300                           (throw 'newroot t))))))
1301       (or *conv* (return t)))
1302     (do ((i (1+ *nn*) (1+ i)))
1303         ((> i degree))
1304       (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*))
1305       (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) *polysc1*)))
1306     (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
1307         ((> i *nn*))
1308       (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) j))
1309       (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) j)))
1310     *nn*))
1311
1312
1313 (defmfun $bfallroots (expr)
1314   (prog (degree *nn* var res $partswitch
1315          ($keepfloat t)
1316          $demoivre
1317          ($listconstvars t)
1318          ($algebraic t) complex $ratfac den expr1)
1319      (setq expr1 (setq expr (meqhk expr)))
1320      (setq var (delete '$%i (cdr ($listofvars expr)) :test #'eq))
1321      (or var (setq var (list (gensym))))
1322      (cond ((not (= (length var) 1))
1323             (merror (intl:gettext "bfallroots: expected a polynomial in one variable; found variables ~M") `((mlist) ,@var)))
1324            ((setq var (car var))))
1325      ;; It's VERY important to set $bftorat to T so that we preserve
1326      ;; the precision of the bfloats when converting them to rationals
1327      ;; for $rat.  Might as well turn off printing of the rat messages
1328      ;; too.
1329      (let (($ratprint nil)
1330            ($bftorat t))
1331        (setq expr ($rat expr '$%i var)
1332              res (reverse (car (cdddar expr)))))
1333      (do ((i (- (length res)
1334                 (length (caddar expr)))
1335              (1- i)))
1336          ((= i 0))
1337        (setq res (cdr res)))
1338      (setq den (cddr expr)
1339            expr (cadr expr))
1340      ;; Check denominator is a complex number
1341      (cond ((numberp den) (setq den (list den 0)))
1342            ((eq (car den) (cadr res))
1343             (setq den (cddr den))
1344             (cond ((numberp (car den))
1345                    (cond ((null (cddr den))
1346                           (setq den (list 0 (car den))))
1347                          ((numberp (caddr den))
1348                           (setq den (list (caddr den) (car den))))
1349                          (t (bf-cpoly-err expr1))))
1350                   (t (bf-cpoly-err expr1))))
1351            (t (bf-cpoly-err expr1)))
1352      ;; If the name variable has disappeared, this is caught here
1353      (setq *nn* 0)
1354      (cond ((numberp expr)
1355             (setq expr (list expr 0)))
1356            ((eq (car expr) (car res))
1357             (setq *nn* 1))
1358            ((eq (car expr) (cadr res))
1359             (setq expr (cddr expr))
1360             (cond ((numberp (car expr))
1361                    (cond ((null (cddr expr))
1362                           (setq expr (list 0 (car expr))))
1363                          ((numberp (caddr expr))
1364                           (setq expr (list (caddr expr) (car expr))))
1365                          (t
1366                           (bf-cpoly-err expr1))))
1367                   (t (bf-cpoly-err expr1))))
1368            (t (bf-cpoly-err expr1)))
1369      (cond ((= *nn* 0)
1370             (cond ($polyfactor
1371                    (let ((*cr* (intofp 0))
1372                          (*ci* (intofp 0)))
1373                      (bf-cdivid-sl (cdr ($bfloat (car expr)))
1374                                    (cdr ($bfloat (cadr expr)))
1375                                    (cdr ($bfloat (car den)))
1376                                    (cdr ($bfloat (cadr den))))
1377                      (return (add (bcons *cr*)
1378                                   (mul '$%i (bcons *ci*))))))
1379                   (t (return (list '(mlist simp)))))))
1380      (setq degree (cadr expr) *nn* (1+ degree))
1381      (setq *pr-sl* (make-array *nn* :initial-element (intofp 0)))
1382      (setq *pi-sl* (make-array *nn* :initial-element (intofp 0)))
1383      (or (catch 'notpoly
1384            (errset (do ((expr (cdr expr) (cddr expr)) (l) (%i (cadr res)))
1385                        ((null expr))
1386                      (setq l (- degree (car expr)) res (cadr expr))
1387                      (cond ((numberp res)
1388                             (setf (aref *pr-sl* l) (cdr ($bfloat res))))
1389                            (t
1390                             (or (eq (car res) %i)
1391                                 (throw 'notpoly nil))
1392                             (setq res (cddr res))
1393                             (setf (aref *pi-sl* l) (cdr ($bfloat (car res))))
1394                             (setq res (caddr res))
1395                             (and res (setf (aref *pr-sl* l) (cdr ($bfloat res))))
1396                             (setq complex t))))))
1397          ;; This should catch expressions like sin(x)-x
1398          (bf-cpoly-err expr1))
1399      (setq *shr-sl* (make-array *nn* :initial-element (intofp 0)))
1400      (setq *shi-sl* (make-array *nn* :initial-element (intofp 0)))
1401      (setq *qpr-sl* (make-array *nn* :initial-element (intofp 0)))
1402      (setq *hr-sl*  (make-array degree :initial-element (intofp 0)))
1403      (setq *qhr-sl* (make-array degree :initial-element (intofp 0)))
1404      (setq *qpi-sl* (make-array *nn* :initial-element (intofp 0)))
1405
1406      (when complex
1407        (setq *hi-sl*  (make-array degree :initial-element (intofp 0)))
1408        (setq *qhi-sl* (make-array degree :initial-element (intofp 0))))
1409      (setq *nn* degree)
1410      (if complex
1411          (setq res (errset (bf-cpoly-sl degree)))
1412          (setq res (bf-rpoly-sl degree)))
1413      (unless res
1414        (mtell (intl:gettext "bfallroots: unexpected error; treat results with caution.~%")))
1415      (when (= *nn* degree)
1416        (merror (intl:gettext "bfallroots: no roots found.")))
1417      (setq res nil)
1418      (cond ((not (zerop *nn*))
1419             (mtell (intl:gettext "bfallroots: only ~S out of ~S roots found.~%") (- degree *nn*) degree)
1420             (setq expr (bcons (intofp 0)))
1421             (do ((i 0 (1+ i)))
1422                 ((> i *nn*))
1423               (setq expr
1424                     (simplify
1425                      (list '(mplus) expr
1426                            (simplify (list '(mtimes)
1427                                            (simplify (list '(mplus)
1428                                                            (simplify (list '(mtimes) '$%i
1429                                                                            (bcons (aref *pi-sl* i))))
1430                                                            (bcons (aref *pr-sl* i))))
1431                                            (simplify (list '(mexpt) var (- *nn* i)))))))))
1432             (setq res (cons expr res)))
1433            ($polyfactor
1434             (setq expr (let ((*cr* (intofp 0))
1435                              (*ci* (intofp 0)))
1436                          (bf-cdivid-sl (aref *pr-sl* 0)
1437                                        (aref *pi-sl* 0)
1438                                        (cdr ($bfloat (car den)))
1439                                        (cdr ($bfloat (cadr den))))
1440                          (add (bcons *cr*) (mul '$%i (bcons *ci*))))
1441                   res (cons expr res))))
1442      (do ((i degree (1- i)))
1443          ((= i *nn*))
1444        ;; zr+%i*zi, where zr and zi parts of the root.
1445        (setq expr (add (bcons (aref *pr-sl* i))
1446                        (mul '$%i (bcons (aref *pi-sl* i)))))
1447        (setq res
1448              (cond ($polyfactor
1449                     (cons (cond ((or complex (fpzerop (aref *pi-sl* i)))
1450                                  (add var (mul -1 expr)))
1451                                 (t
1452                                  (setq i (1- i))
1453                                  (simplify (list '(mplus)
1454                                                  (simplify (list '(mexpt) var 2))
1455                                                  (simplify (list '(mtimes) var
1456                                                                  (bcons (aref *pr-sl* i))))
1457                                                  (bcons (aref *pr-sl* (1+ i)))))))
1458                           res))
1459                    (t
1460                     (cons (let ((expr (simplify (list '(mequal) var expr))))
1461                             (if $programmode expr (displine expr)))
1462                           res)))))
1463      (return (simplify (if $polyfactor
1464                            (cons '(mtimes) res)
1465                            (cons '(mlist) (nreverse res)))))))
1466
1467 (defun bf-rpoly-sl (degree)
1468   (let ((*logbas* (fplog (intofp 2)))
1469         (*infin* (intofp +most-positive-flonum+))
1470         (*are* (bf-scale-float (intofp 2) (- fpprec)))
1471         (*mre* (intofp 0))
1472         (xx (fproot bfhalf 2))
1473         (yy (intofp 0))
1474         ;; cos(94deg)
1475         (cosr (intofp
1476                -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804))
1477         ;; sin(94deg)
1478         (sinr (intofp
1479                0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315))
1480         (aa (intofp 0))
1481         (cc (intofp 0))
1482         (bb (intofp 0))
1483         (bnd (intofp 0))
1484         (*sr* (intofp 0))
1485         (*u* (intofp 0))
1486         (*v* (intofp 0))
1487         (t1 (intofp 0))
1488         (*szr* (intofp 0))
1489         (*szi* (intofp 0))
1490         (*lzr* (intofp 0))
1491         (*lzi* (intofp 0))
1492         (*nz* 0)
1493         (*n* 0)
1494         (*polysc* 0)
1495         (*polysc1* 0)
1496         (zerok 0)
1497         (conv1 t))
1498     (setq *mre* *are*
1499           yy (fpminus xx))
1500     (do ((i degree (1- i)))
1501         ((not (fpzerop (aref *pr-sl* i)))
1502          (setq *nn* i *n* (1- i))))
1503     (setq degree *nn*)
1504     (do ((i 0 (1+ i)))
1505         ((> i *nn*))
1506       (setf (aref *shr-sl* i) (fpabs (aref *pr-sl* i))))
1507     (if (> *nn* 0) (bf-scale-sl))
1508     (do nil
1509         ((< *nn* 3)
1510          (cond ((= *nn* 2)
1511                 (bf-quad-sl (aref *pr-sl* 0.) (aref *pr-sl* 1) (aref *pr-sl* 2))
1512                 (cond ((and $polyfactor (not (fpzerop *szi*)))
1513                        (setf (aref *pr-sl* 2) (fpquotient (aref *pr-sl* 2)
1514                                                           (aref *pr-sl* 0)))
1515                        (setf (aref *pr-sl* 1) (fpquotient (aref *pr-sl* 1)
1516                                                           (aref *pr-sl* 0)))
1517                        (setf (aref *pi-sl* 2) (intofp 1)))
1518                       (t (setf (aref *pr-sl* 2) *szr*)
1519                          (setf (aref *pi-sl* 2) *szi*)
1520                          (setf (aref *pr-sl* 1) *lzr*)
1521                          (setf (aref *pi-sl* 1) *lzi*))))
1522                (t
1523                 (setf (aref *pr-sl* 1) (fpminus (fpquotient (aref *pr-sl* 1)
1524                                                             (aref *pr-sl* 0))))))
1525          (setq *nn* 0))
1526       (do ((i 0 (1+ i)))
1527           ((> i *nn*))
1528         (setf (aref *shr-sl* i) (fpabs (aref *pr-sl* i))))
1529       (setq bnd (bf-cauchy-sl))
1530       (do ((i 1 (1+ i)))
1531           ((> i *n*))
1532         (setf (aref *hr-sl* i)
1533               (fpquotient (fptimes* (intofp (- *n* i))
1534                                     (aref *pr-sl* i))
1535                           (intofp *n*))))
1536       (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
1537       (setq aa (aref *pr-sl* *nn*)
1538             bb (aref *pr-sl* *n*)
1539             zerok (fpzerop (aref *hr-sl* *n*)))
1540       (do ((jj 1 (1+ jj)))
1541           ((> jj 5.))
1542         (setq cc (aref *hr-sl* *n*))
1543         (cond (zerok (do ((j *n* (1- j)))
1544                          ((< j 1))
1545                        (setf (aref *hr-sl* j) (aref *hr-sl* (1- j))))
1546                      (setf (aref *hr-sl* 0) (intofp 0))
1547                      (setq zerok (fpzerop (aref *hr-sl* *n*))))
1548               (t
1549                (setq t1 (fpminus (fpquotient aa cc)))
1550                (do ((j *n* (1- j)))
1551                    ((< j 1))
1552                  (setf (aref *hr-sl* j)
1553                        (fpplus (fptimes* t1 (aref *hr-sl* (1- j)))
1554                                (aref *pr-sl* j))))
1555                (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
1556                (setq zerok (not (fpgreaterp (fpabs (aref *hr-sl* *n*))
1557                                             (fptimes* (fpabs bb)
1558                                                       (fptimes* *are* (intofp 10)))))))))
1559       (do ((i 0 (1+ i)))
1560           ((> i *n*))
1561         (setf (aref *shi-sl* i) (aref *hr-sl* i)))
1562       (do ((cnt 1 (1+ cnt)))
1563           ((> cnt 20.)
1564            (setq conv1 nil))
1565         (setq xx (prog1
1566                      (fpdifference (fptimes* cosr xx)
1567                                    (fptimes* sinr yy))
1568                    (setq yy (fpplus (fptimes* sinr xx)
1569                                     (fptimes* cosr yy))))
1570               *sr* (fptimes* bnd xx)
1571               *u* (fptimes* (intofp -2) *sr*)
1572               *v* bnd)
1573         (bf-fxshfr-sl (* 20 cnt))
1574         (cond ((> *nz* 0)
1575                (setf (aref *pr-sl* *nn*) *szr*)
1576                (setf (aref *pi-sl* *nn*) *szi*)
1577                (cond ((= *nz* 2)
1578                       (setf (aref *pr-sl* *n*) *lzr*)
1579                       (setf (aref *pi-sl* *n*) *lzi*)
1580                       (cond ((and $polyfactor (not (fpzerop *szi*)))
1581                              (setf (aref *pr-sl* *nn*) *v*)
1582                              (setf (aref *pr-sl* *n*) *u*)
1583                              (setf (aref *pi-sl* *nn*) (intofp 1))))))
1584                (setq *nn* (- *nn* *nz*) *n* (1- *nn*))
1585                (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *pr-sl* i) (aref *qpr-sl* i)))
1586                (return nil)))
1587         (do ((i 0 (1+ i))) ((> i *n*)) (setf (aref *hr-sl* i) (aref *shi-sl* i))))
1588       (or conv1 (return nil)))
1589     (cond ($polyfactor
1590            (do ((i degree (1- i)))
1591                ((= i *nn*))
1592              (cond ((fpzerop (aref *pi-sl* i))
1593                     (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*)))
1594                    (t (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) (* 2 *polysc1*)))
1595                       (setq i (1- i))
1596                       (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*))))))
1597           (t (do ((i (1+ *nn*) (1+ i)))
1598                  ((> i degree))
1599                (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*))
1600                (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) *polysc1*)))))
1601     (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
1602         ((> i *nn*))
1603       (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) j)))
1604     t))
1605
1606 (defun bf-fxshfr-sl (l2)
1607   (let ((*my-type* 0)
1608         (*a* (intofp 0))
1609         (*b* (intofp 0))
1610         (*c* (intofp 0))
1611         (*d* (intofp 0))
1612         (*e* (intofp 0))
1613         (*f* (intofp 0))
1614         (*g* (intofp 0))
1615         (*h* (intofp 0))
1616         (*a1* (intofp 0))
1617         (*a3* (intofp 0))
1618         (*a7* (intofp 0)))
1619     (declare (special *my-type*))
1620     (setq *nz* 0)
1621     (bf-quadsd-sl)
1622     (bf-calcsc-sl)
1623     (do ((j 1 (1+ j))
1624          (betav (intofp 0.25))
1625          (betas (intofp 0.25))
1626          (oss *sr*)
1627          (ovv *v*)
1628          (tvv) (tss) (ss) (vv) (tv) (ts) (ots) (otv)
1629          (*ui*) (*vi*) (*s*) (svv) (svu) (iflag) (vpass) (spass) (vtry) (stry))
1630         ((> j l2))
1631       (bf-nextk-sl)
1632       (bf-calcsc-sl)
1633       (bf-newest-sl)
1634       (setq vv *vi*
1635             ss (intofp 0))
1636       (or (fpzerop (aref *hr-sl* *n*))
1637           (setq ss (fpminus (fpquotient (aref *pr-sl* *nn*)
1638                                         (aref *hr-sl* *n*)))))
1639       (setq tv (intofp 1)
1640             ts (intofp 1))
1641       (cond ((not (or (= j 1)
1642                       (= *my-type* 3)))
1643              (or (fpzerop vv)
1644                  (setq tv (fpabs (fpquotient (fpdifference vv ovv)
1645                                              vv))))
1646              (or (fpzerop ss)
1647                  (setq ts (fpabs (fpquotient (fpdifference ss oss)
1648                                              ss))))
1649              (setq tvv (intofp 1))
1650              (and (not (fpgreaterp tv otv))
1651                   (setq tvv (fptimes* tv otv)))
1652              (setq tss (intofp 1))
1653              (and (not (fpgreaterp ts ots))
1654                   (setq tss (fptimes* ts ots)))
1655              (setq vpass (not (fpgreaterp tvv betav))
1656                    spass (not (fpgreaterp tss betas)))
1657              (cond ((or spass vpass)
1658                     (setq svu *u* svv *v*)
1659                     (do ((i 0 (1+ i)))
1660                         ((> i *n*))
1661                       (setf (aref *shr-sl* i)
1662                             (aref *hr-sl* i)))
1663                     (setq *s* ss vtry nil stry nil)
1664                     (and (do ((bool (not (and spass (or (not vpass)
1665                                                         (not (fpgreaterp tss tvv)))))
1666                                     t)
1667                               (l50 nil nil))
1668                              (nil)
1669                            (cond (bool
1670                                   (bf-quadit-sl)
1671                                   (and (> *nz* 0) (return t))
1672                                   (setq vtry t
1673                                         betav (fptimes* (intofp 0.25) betav))
1674                                   (cond ((or stry (not spass))
1675                                          (setq l50 t))
1676                                         (t (do ((i 0 (1+ i)))
1677                                                ((> i *n*))
1678                                              (setf (aref *hr-sl* i)
1679                                                    (aref *shr-sl* i)))))))
1680                            (cond ((not l50)
1681                                   (setq iflag (bf-realit-sl))
1682                                   (and (> *nz* 0) (return t))
1683                                   (setq stry t
1684                                         betas (fptimes* (intofp 0.25) betas))
1685                                   (cond ((zerop iflag)
1686                                          (setq l50 t))
1687                                         (t
1688                                          (setq *ui* (fpminus (fpplus *s* *s*))
1689                                                *vi* (fptimes* *s* *s*))))))
1690                            (cond (l50
1691                                   (setq *u* svu *v* svv)
1692                                   (do ((i 0 (1+ i)))
1693                                       ((> i *n*))
1694                                     (setf (aref *hr-sl* i)
1695                                           (aref *shr-sl* i)))
1696                                   (and (or (not vpass) vtry)
1697                                        (return nil)))))
1698                          (return nil))
1699                     (bf-quadsd-sl)
1700                     (bf-calcsc-sl)))))
1701       (setq ovv vv
1702             oss ss
1703             otv tv
1704             ots ts))))
1705
1706 (defun bf-quadit-sl nil
1707   (setq *nz* 0 *u* *ui* *v* *vi*)
1708   (do ((tried)
1709        (j 0)
1710        (ee)
1711        (mp)
1712        (relstp)
1713        (omp)
1714        (ms))
1715       (nil)
1716     (bf-quad-sl (intofp 1) *u* *v*)
1717     (and (fpgreaterp (fpabs (fpdifference (fpabs *szr*)
1718                                           (fpabs *lzr*)))
1719                      (fptimes* (intofp 1e-2) (fpabs *lzr*)))
1720          (return nil))
1721     (bf-quadsd-sl)
1722     (setq mp (fpplus (fpabs (fpdifference *a* (fptimes* *szr* *b*)))
1723                      (fpabs (fptimes* *szi* *b*)))
1724           ms (bf-cmod-sl *szr* *szi*)
1725           ee (bf-errev-sl ms mp))
1726     (cond ((not (fpgreaterp mp (fptimes* (intofp 2e1) ee)))
1727            (setq *nz* 2)
1728            (return nil)))
1729     (setq j (1+ j))
1730     (and (> j 20) (return nil))
1731     (cond ((not (or (< j 2)
1732                     (fpgreaterp relstp (intofp 1e-2))
1733                     (not (fpgreaterp mp omp))
1734                     tried))
1735            (and (not (fpgreaterp relstp *are*))
1736                 (setq relstp *are*))
1737            (setq relstp (fpsqrt relstp)
1738                  *u* (fpdifference *u* (fptimes* *u* relstp))
1739                  *v* (fpplus *v* (fptimes* *v* relstp)))
1740            (bf-quadsd-sl)
1741            (do ((i 1 (1+ i)))
1742                ((> i 5))
1743              (bf-calcsc-sl)
1744              (bf-nextk-sl))
1745            (setq tried t j 0)))
1746     (setq omp mp)
1747     (bf-calcsc-sl)
1748     (bf-nextk-sl)
1749     (bf-calcsc-sl)
1750     (bf-newest-sl)
1751     (and (fpzerop *vi*) (return nil))
1752     (setq relstp (fpabs (fpquotient (fpdifference *vi* *v*)
1753                                     *vi*))
1754           *u* *ui*
1755           *v* *vi*)))
1756
1757 (defun bf-realit-sl ()
1758   (setq *nz* 0)
1759   (do ((j 0)
1760        (pv)
1761        (ee)
1762        (ms)
1763        (mp)
1764        (kv)
1765        (t1)
1766        (omp))
1767       (nil)
1768     (setq pv (aref *pr-sl* 0))
1769     (setf (aref *qpr-sl* 0) pv)
1770     (do ((i 1 (1+ i)))
1771         ((> i *nn*))
1772       (setq pv (fpplus (fptimes* pv *s*)
1773                        (aref *pr-sl* i)))
1774       (setf (aref *qpr-sl* i) pv))
1775     (setq mp (fpabs pv)
1776           ms (fpabs *s*)
1777           ee (fptimes* (fpquotient *mre* (fpplus *are* *mre*))
1778                        (fpabs (aref *qpr-sl* 0))))
1779     (do ((i 1 (1+ i)))
1780         ((> i *nn*))
1781       (setq ee (fpplus (fptimes* ee ms)
1782                        (fpabs (aref *qpr-sl* i)))))
1783     (cond ((not (fpgreaterp mp
1784                             (fptimes* (intofp 2e1)
1785                                       (fpdifference (fptimes* (fpplus *are* *mre*)
1786                                                               ee)
1787                                                     (fptimes* *mre* mp)))))
1788            (setq *nz* 1 *szr* *s* *szi* (intofp 0))
1789            (return 0)))
1790     (setq j (1+ j))
1791     (and (> j 10) (return 0))
1792     (cond ((not (or (< j 2)
1793                     (fpgreaterp (fpabs t1)
1794                                 (fptimes* (intofp 1e-3) (fpabs (fpdifference *s* t1))))
1795                     (not (fpgreaterp mp omp))))
1796            (return 1)))
1797     (setq omp mp kv (aref *hr-sl* 0))
1798     (setf (aref *qhr-sl* 0) kv)
1799     (do ((i 1 (1+ i)))
1800         ((> i *n*))
1801       (setq kv (fpplus (fptimes* kv *s*)
1802                        (aref *hr-sl* i)))
1803       (setf (aref *qhr-sl* i) kv))
1804     (cond ((fpgreaterp (fpabs kv)
1805                        (fptimes* (fpabs (aref *hr-sl* *n*))
1806                                  (fptimes* (intofp 1e1) *are*)))
1807            (setq t1 (fpminus (fpquotient pv kv)))
1808            (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
1809            (do ((i 1 (1+ i)))
1810                ((> i *n*))
1811              (setf (aref *hr-sl* i)
1812                    (fpplus (fptimes* t1 (aref *qhr-sl* (1- i)))
1813                            (aref *qpr-sl* i)))))
1814           (t
1815            (setf (aref *hr-sl* 0) (intofp 0))
1816            (do ((i 1 (1+ i)))
1817                ((> i *n*))
1818              (setf (aref *hr-sl* i) (aref *qhr-sl* (1- i))))))
1819     (setq kv (aref *hr-sl* 0))
1820     (do ((i 1 (1+ i)))
1821         ((> i *n*))
1822       (setq kv (fpplus (fptimes* kv *s*)
1823                        (aref *hr-sl* i))))
1824     (setq t1 (intofp 0))
1825     (and (fpgreaterp (fpabs kv)
1826                      (fptimes* (fpabs (aref *hr-sl* *n*))
1827                                (fptimes* (intofp 10) *are*)))
1828          (setq t1 (fpminus (fpquotient pv kv))))
1829     (setq *s* (fpplus *s* t1))))
1830
1831 (defun bf-calcsc-sl ()
1832   (declare (special *my-type*))
1833   (setq *d* (aref *hr-sl* 0))
1834   (setf (aref *qhr-sl* 0) *d*)
1835   (setq *c* (fpdifference (aref *hr-sl* 1)
1836                           (fptimes* *u* *d*)))
1837   (setf (aref *qhr-sl* 1) *c*)
1838   (do ((i 2 (1+ i))
1839        (c0))
1840       ((> i *n*))
1841     (setq c0 (fpdifference (fpdifference (aref *hr-sl* i)
1842                                          (fptimes* *u* *c*))
1843                            (fptimes* *v* *d*)))
1844     (setf (aref *qhr-sl* i) c0)
1845     (setq *d* *c* *c* c0))
1846   (cond ((not (or (fpgreaterp (fpabs *c*)
1847                               (fptimes* (fpabs (aref *hr-sl* *n*))
1848                                         (fptimes* (intofp 100) *are*)))
1849                   (fpgreaterp (fpabs *d*)
1850                               (fptimes* (fpabs (aref *hr-sl* (1- *n*)))
1851                                         (fptimes* (intofp 100) *are*)))))
1852          (setq *my-type* 3))
1853         ((not (not (fpgreaterp (fpabs *d*) (fpabs *c*))))
1854          (setq *my-type* 2
1855                *e* (fpquotient *a* *d*)
1856                *f* (fpquotient *c* *d*)
1857                *g* (fptimes* *u* *b*)
1858                *h* (fptimes* *v* *b*)
1859                *a3* (fpplus (fptimes* (fpplus *a* *g*) *e*)
1860                             (fptimes* *h* (fpquotient *b* *d*)))
1861                *a1* (fpdifference (fptimes* *b* *f*)
1862                                   *a*)
1863                *a7* (fpplus (fptimes* (fpplus *f* *u*) *a*)
1864                             *h*)))
1865         (t (setq *my-type* 1
1866                  *e* (fpquotient *a* *c*)
1867                  *f* (fpquotient *d* *c*)
1868                  *g* (fptimes* *u* *e*)
1869                  *h* (fptimes* *v* *b*)
1870                  *a3* (fpplus (fptimes* *a* *e*)
1871                               (fptimes* (fpplus (fpquotient *h* *c*)
1872                                                 *g*)
1873                                         *b*))
1874                  *a1* (fpdifference *b*
1875                                     (fptimes* *a* (fpquotient *d* *c*)))
1876                  *a7* (fpplus *a*
1877                               (fpplus (fptimes* *g* *d*)
1878                                       (fptimes* *h* *f*))))))
1879   nil)
1880
1881 (defun bf-nextk-sl ()
1882   (declare (special *my-type*))
1883   (cond ((= *my-type* 3)
1884          (setf (aref *hr-sl* 0) (intofp 0))
1885          (setf (aref *hr-sl* 1) (intofp 0))
1886          (do ((i 2 (1+ i)))
1887              ((> i *n*))
1888            (setf (aref *hr-sl* i) (aref *qhr-sl* (- i 2)))))
1889         ((fpgreaterp (fpabs *a1*)
1890                      (fptimes* (fpabs (if (= *my-type* 1) *b* *a*))
1891                                (fptimes* (intofp 1e1) *are*)))
1892          (setq *a7* (fpquotient *a7* *a1*)
1893                *a3* (fpquotient *a3* *a1*))
1894          (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
1895          (setf (aref *hr-sl* 1)
1896                (fpdifference (aref *qpr-sl* 1)
1897                              (fptimes* *a7* (aref *qpr-sl* 0))))
1898          (do ((i 2 (1+ i)))
1899              ((> i *n*))
1900            (setf (aref *hr-sl* i)
1901                  (fpplus (fptimes* *a3* (aref *qhr-sl* (- i 2)))
1902                          (fpplus (fpminus (fptimes* *a7* (aref *qpr-sl* (1- i))))
1903                                  (aref *qpr-sl* i))))))
1904         (t
1905          (setf (aref *hr-sl* 0) (intofp 0))
1906          (setf (aref *hr-sl* 1) (fpminus (fptimes* *a7* (aref *qpr-sl* 0))))
1907          (do ((i 2 (1+ i)))
1908              ((> i *n*))
1909            (setf (aref *hr-sl* i)
1910                  (fpdifference (fptimes* *a3* (aref *qhr-sl* (- i 2)))
1911                                (fptimes* *a7* (aref *qpr-sl* (1- i))))))))
1912   nil)
1913
1914 (defun bf-newest-sl ()
1915   (declare (special *my-type*))
1916   (let ((a4 (intofp 0))
1917         (a5 (intofp 0))
1918         (b1 (intofp 0))
1919         (b2 (intofp 0))
1920         (c1 (intofp 0))
1921         (c2 (intofp 0))
1922         (c3 (intofp 0))
1923         (c4 (intofp 0)))
1924     (cond ((= *my-type* 3)
1925            (setq *ui* (intofp 0)
1926                  *vi* (intofp 0)))
1927           (t
1928            (if (= *my-type* 2)
1929                (setq a4 (fpplus (fptimes* (fpplus *a* *g*)
1930                                           *f*)
1931                                 *h*)
1932                      a5 (fpplus (fptimes* (fpplus *f* *u*)
1933                                           *c*)
1934                                 (fptimes* *v* *d*)))
1935                (setq a4 (fpplus *a*
1936                                 (fpplus (fptimes* *u* *b*)
1937                                         (fptimes* *h* *f*)))
1938                      a5 (fpplus *c*
1939                                 (fptimes* (fpplus *u* (fptimes* *v* *f*))
1940                                           *d*))))
1941            (setq b1 (fpminus (fpquotient (aref *hr-sl* *n*)
1942                                          (aref *pr-sl* *nn*)))
1943                  b2 (fpminus (fpquotient (fpplus (aref *hr-sl* (1- *n*))
1944                                                  (fptimes* b1 (aref *pr-sl* *n*)))
1945                                          (aref *pr-sl* *nn*)))
1946                  c1 (fptimes* (fptimes* *v* b2) *a1*)
1947                  c2 (fptimes* b1 *a7*)
1948                  c3 (fptimes* (fptimes* b1 b1) *a3*)
1949                  c4 (fpdifference (fpdifference c1 c2) c3)
1950                  c1 (fpplus (fpplus a5 (fptimes* b1 a4))
1951                             (fpminus c4)))
1952            (if (fpzerop c1)
1953                (setq *ui* (intofp 0)
1954                      *vi* (intofp 0))
1955                (setq *ui* (fpdifference
1956                            *u*
1957                            (fpquotient (fpplus (fptimes* *u* (fpplus c3 c2))
1958                                                (fptimes* *v*
1959                                                          (fpplus (fptimes* b1 *a1*)
1960                                                                  (fptimes* b2 *a7*))))
1961                                        c1))
1962                      *vi* (fptimes* *v*
1963                                     (fpplus (fpone) (fpquotient c4 c1)))))))
1964     nil))
1965
1966 (defun bf-quadsd-sl ()
1967   (setq *b* (aref *pr-sl* 0))
1968   (setf (aref *qpr-sl* 0) *b*)
1969   (setq *a* (fpdifference (aref *pr-sl* 1)
1970                           (fptimes* *u* *b*)))
1971   (setf (aref *qpr-sl* 1) *a*)
1972   (do ((i 2 (1+ i))
1973        (c0))
1974       ((> i *nn*))
1975     (setq c0 (fpdifference (fpdifference (aref *pr-sl* i)
1976                                          (fptimes* *u* *a*))
1977                            (fptimes* *v* *b*)))
1978     (setf (aref *qpr-sl* i) c0)
1979     (setq *b* *a*
1980           *a* c0)))
1981
1982 (defun bf-quad-sl (a0 b1 c0)
1983   (setq *szr* (intofp 0)
1984         *szi* (intofp 0)
1985         *lzr* (intofp 0)
1986         *lzi* (intofp 0))
1987   (let ((b0 (intofp 0))
1988         (l0 (intofp 0))
1989         (*e* (intofp 0)))
1990     (cond ((fpzerop a0)
1991            (or (fpzerop b1)
1992                (setq *szr* (fpminus (fpquotient c0 b1)))))
1993           ((fpzerop c0)
1994            (setq *lzr* (fpminus (fpquotient b1 a0))))
1995           (t
1996            (setq b0 (fpquotient b1 (intofp 2)))
1997            (cond ((not (fpgreaterp (fpabs b0) (fpabs c0)))
1998                   (setq *e* a0)
1999                   (and (not (fpgreaterp c0 (intofp 0)))
2000                        (setq *e* (fpminus a0)))
2001                   (setq *e* (fpdifference (fptimes* b0 (fpquotient b0 (fpabs c0)))
2002                                           *e*)
2003                         l0 (fptimes* (fpsqrt (fpabs *e*))
2004                                      (fpsqrt (fpabs c0)))))
2005                  (t (setq *e* (fpdifference (intofp 1)
2006                                             (fptimes* (fpquotient a0 b0)
2007                                                       (fpquotient c0 b0)))
2008                           l0 (fptimes* (fpsqrt (fpabs *e*))
2009                                        (fpabs b0)))))
2010            (cond ((not (fpgreaterp *e* (intofp 0)))
2011                   (setq *szr* (fpminus (fpquotient b0 a0))
2012                         *lzr* *szr*
2013                         *szi* (fpabs (fpquotient l0 a0))
2014                         *lzi* (fpminus *szi*)))
2015                  (t (or (not (fpgreaterp b0 (intofp 0)))
2016                         (setq l0 (fpminus l0)))
2017                     (setq *lzr* (fpquotient (fpdifference l0 b0) a0))
2018                     (or (fpzerop *lzr*)
2019                         (setq *szr* (fpquotient (fpquotient c0 *lzr*) a0)))))))
2020     nil))