src/cpoly.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancments.                    ;;;;;
   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 $partswitch $keepfloat $demoivre $listconstvars
  25                       $algebraic $ratfac $programmode))
  26
  27 (declare-top (special *logbas* *infin* *are* *mre* *cr* *ci* *sr* *si*
  28                       *tr* *ti* *zr* *zi* *n* *nn* *bool*
  29                       *conv* *pvr* *pvi* *polysc* *polysc1*))
  30
  31 (declare-top (special *pr-sl* *pi-sl* *shr-sl* *shi-sl* *qpr-sl* *qpi-sl* *hr-sl*
  32                       *hi-sl* *qhr-sl* *qhi-sl*))
  33
  34 (declare-top (special *u* *v* *a* *b* *c* *d* *a1* *a3* *a7* *e* *f* *g* *h*
  35                       *szr* *szi* *lzr* *lzi* *nz* *ui* *vi* *s*))
  36
  37 (defmvar $polyfactor nil
  38   "When T factor the polynomial over the real or complex numbers.")
  39
  40 (defmfun $allroots (expr)
  41   (prog (degree *nn* var res $partswitch $keepfloat $demoivre $listconstvars
  42          $algebraic complex $ratfac den expr1)
  43      (setq $keepfloat t $listconstvars t $algebraic t)
  44      (setq expr1 (setq expr (meqhk expr)))
  45      (setq var (delete '$%i (cdr ($listofvars expr)) :test #'eq))
  46      (or var (setq var (list (gensym))))
  47      (cond ((not (= (length var) 1))
  48             (merror (intl:gettext "allroots: expected a polynomial in one variable; found variables ~M") `((mlist) ,@var)))
  49            ((setq var (car var))))
  50      (setq expr ($rat expr '$%i var)
  51            res (reverse (car (cdddar expr))))
  52      (do ((i (- (length res) (length (caddar expr))) (1- i)))
  53          ((= i 0))
  54        (setq res (cdr res)))
  55      (setq den (cddr expr) expr (cadr expr))
  56 ;;;check denominator is a complex number
  57      (cond ((numberp den) (setq den (list den 0)))
  58            ((eq (car den) (cadr res))
  59             (setq den (cddr den))
  60             (cond ((numberp (car den))
  61                    (cond ((null (cddr den)) (setq den (list 0 (car den))))
  62                          ((numberp (caddr den))
  63                           (setq den (list (caddr den) (car den))))
  64                          (t (cpoly-err expr1))))
  65                   (t (cpoly-err expr1))))
  66            (t (cpoly-err expr1)))
  67 ;;;if the name variable has disappeared, this is caught here
  68      (setq *nn* 0)
  69      (cond ((numberp expr) (setq expr (list expr 0)))
  70            ((eq (car expr) (car res)) (setq *nn* 1))
  71            ((eq (car expr) (cadr res))
  72             (setq expr (cddr expr))
  73             (cond ((numberp (car expr))
  74                    (cond ((null (cddr expr)) (setq expr (list 0 (car expr))))
  75                          ((numberp (caddr expr))
  76                           (setq expr (list (caddr expr) (car expr))))
  77                          (t (cpoly-err expr1))))
  78                   (t (cpoly-err expr1))))
  79            (t (cpoly-err expr1)))
  80      (cond ((= *nn* 0)
  81             (cond ($polyfactor
  82                    (let ((*cr* 0.0) (*ci* 0.0))
  83                      (cdivid-sl (float (car expr)) (float (cadr expr))
  84                                 (float (car den)) (float (cadr den)))
  85                      (return (simplify (list '(mplus)
  86                                              (simplify (list '(mtimes) '$%i *ci*))
  87                                              *cr*)))))
  88                   (t (return (list '(mlist simp)))))))
  89      (setq degree (cadr expr) *nn* (1+ degree))
  90      (setq *pr-sl* (make-array *nn* :initial-element 0.0))
  91      (setq *pi-sl* (make-array *nn* :initial-element 0.0))
  92      (or (catch 'notpoly
  93            (errset (do ((expr (cdr expr) (cddr expr)) (l) (%i (cadr res)))
  94                        ((null expr))
  95                      (setq l (- degree (car expr)) res (cadr expr))
  96                      (cond ((numberp res) (setf (aref *pr-sl* l) (float res)))
  97                            (t
  98                             (or (eq (car res) %i)
  99                                 (throw 'notpoly nil))
 100                             (setq res (cddr res))
 101                             (setf (aref *pi-sl* l) (float (car res)))
 102                             (setq res (caddr res))
 103                             (and res (setf (aref *pr-sl* l) (float res)))
 104                             (setq complex t))))))
 105          ;;this should catch expressions like sin(x)-x
 106          (cpoly-err expr1))
 107      (setq *shr-sl* (make-array *nn* :initial-element 0.0))
 108      (setq *shi-sl* (make-array *nn* :initial-element 0.0))
 109      (setq *qpr-sl* (make-array *nn* :initial-element 0.0))
 110      (setq *hr-sl*  (make-array degree :initial-element 0.0))
 111      (setq *qhr-sl* (make-array degree :initial-element 0.0))
 112      (setq *qpi-sl* (make-array *nn* :initial-element 0.0))
 113      (when complex
 114        (setq *hi-sl*  (make-array degree :initial-element 0.0))
 115        (setq *qhi-sl* (make-array degree :initial-element 0.0)))
 116      (setq *nn* degree)
 117      (if complex
 118          (setq res (errset (cpoly-sl degree)))
 119          (setq res (errset (rpoly-sl degree))))
 120      (unless res
 121        (mtell (intl:gettext "allroots: unexpected error; treat results with caution.~%")))
 122      (when (= *nn* degree)
 123        (merror (intl:gettext "allroots: no roots found.")))
 124      (setq res nil)
 125      (cond ((not (zerop *nn*))
 126             (mtell (intl:gettext "allroots: only ~S out of ~S roots found.~%") (- degree *nn*) degree)
 127             (setq expr 0.0)
 128             (do ((i 0 (1+ i)))
 129                 ((> i *nn*))
 130               (setq expr
 131                     (simplify
 132                      (list '(mplus) expr
 133                            (simplify (list '(mtimes)
 134                                            (simplify (list '(mplus)
 135                                                            (simplify (list '(mtimes) '$%i (aref *pi-sl* i)))
 136                                                            (aref *pr-sl* i)))
 137                                            (simplify (list '(mexpt) var (- *nn* i)))))))))
 138             (setq res (cons expr res)))
 139            ($polyfactor
 140             (setq expr (let ((*cr* 0.0) (*ci* 0.0))
 141                          (cdivid-sl (aref *pr-sl* 0) (aref *pi-sl* 0)
 142                                     (float (car den))
 143                                     (float (cadr den)))
 144                          (simplify (list '(mplus) (simplify (list '(mtimes) '$%i *ci*)) *cr*)))
 145                   res (cons expr res))))
 146      (do ((i degree (1- i)))
 147          ((= i *nn*))
 148        (setq expr (simplify (list '(mplus)
 149                                   (simplify (list '(mtimes) '$%i (aref *pi-sl* i)))
 150                                   (aref *pr-sl* i))))
 151        (setq res
 152              (cond ($polyfactor (cons (cond ((or complex (zerop (aref *pi-sl* i)))
 153                                              (simplify (list '(mplus) var (simplify (list '(mminus) expr)))))
 154                                             (t (setq i (1- i))
 155                                                (simplify (list '(mplus)
 156                                                                (simplify (list '(mexpt) var 2))
 157                                                                (simplify (list '(mtimes) var
 158                                                                                (aref *pr-sl* i)))
 159                                                                (aref *pr-sl* (1+ i))))))
 160                                       res))
 161                    ((cons (let ((expr (simplify (list '(mequal) var expr))))
 162                             (if $programmode expr (displine expr)))
 163                           res)))))
 164      (return (simplify (if $polyfactor
 165                            (cons '(mtimes) res)
 166                            (cons '(mlist) (nreverse res)))))))
 167
 168 (defun cpoly-err (expr)
 169   (merror (intl:gettext "allroots: expected a polynomial; found ~M") expr))
 170
 171 (defun cpoly-sl (degree)
 172   (let ((*logbas* (log 2.0))
 173         (*infin* most-positive-flonum)
 174         (*are* flonum-epsilon)
 175         (*mre* 0.0)
 176         (xx (sqrt 0.5))
 177         (yy 0.0)
 178         (cosr (cos (float (* 94/180 pi))))
 179         (sinr (sin (float (* 94/180 pi))))
 180         (*cr* 0.0) (*ci* 0.0)
 181         (*sr* 0.0) (*si* 0.0)
 182         (*tr* 0.0) (*ti* 0.0)
 183         (*zr* 0.0) (*zi* 0.0)
 184         (bnd 0.0)
 185         (*n* 0)
 186         (*polysc* 0)
 187         (*polysc1* 0)
 188         (*conv* nil))
 189     (setq *mre* (* 2.0 (sqrt 2.0) *are*)
 190           yy (- xx))
 191     (do ((i degree (1- i)))
 192         ((not (and (zerop (aref *pr-sl* i)) (zerop (aref *pi-sl* i))))
 193          (setq *nn* i
 194                *n* (1- i))))
 195     (setq degree *nn*)
 196     (do ((i 0 (1+ i)))
 197         ((> i *nn*))
 198       (setf (aref *shr-sl* i) (cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
 199     (if (> *nn* 0) (scale-sl))
 200     (do ()
 201         ((> 2 *nn*)
 202          (cdivid-sl (- (aref *pr-sl* 1)) (- (aref *pi-sl* 1))
 203                     (aref *pr-sl* 0) (aref *pi-sl* 0))
 204          (setf (aref *pr-sl* 1) *cr*)
 205          (setf (aref *pi-sl* 1) *ci*)
 206          (setq *nn* 0))
 207       (do ((i 0 (1+ i)))
 208           ((> i *nn*))
 209         (setf (aref *shr-sl* i) (cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
 210       (setq bnd (cauchy-sl))
 211       (catch 'newroot
 212         (do ((cnt1 1 (1+ cnt1)))
 213             ((> cnt1 2))
 214           (noshft-sl 5)
 215           (do ((cnt2 1 (1+ cnt2)))
 216               ((> cnt2 9))
 217             (setq xx (prog1
 218                          (- (* cosr xx) (* sinr yy))
 219                        (setq yy (+ (* sinr xx) (* cosr yy))))
 220                   *sr* (* bnd xx)
 221                   *si* (* bnd yy))
 222             (fxshft-sl (* 10 cnt2))
 223             (cond (*conv* (setf (aref *pr-sl* *nn*) *zr*)
 224                           (setf (aref *pi-sl* *nn*) *zi*)
 225                           (setq *nn* *n* *n* (1- *n*))
 226                           (do ((i 0 (1+ i)))
 227                               ((> i *nn*))
 228                             (setf (aref *pr-sl* i) (aref *qpr-sl* i))
 229                             (setf (aref *pi-sl* i) (aref *qpi-sl* i)))
 230                           (throw 'newroot t))))))
 231       (or *conv* (return t)))
 232     (do ((i (1+ *nn*) (1+ i)))
 233         ((> i degree))
 234       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*))
 235       (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) *polysc1*)))
 236     (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
 237         ((> i *nn*))
 238       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))
 239       (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) j)))
 240     *nn*))
 241
 242 (defun noshft-sl (l1)
 243   (do ((i 0 (1+ i))
 244        (xni (float *nn*) (1- xni))
 245        (t1 (/ (float *nn*))))
 246       ((> i *n*))
 247     (setf (aref *hr-sl* i) (* (aref *pr-sl* i) xni t1))
 248     (setf (aref *hi-sl* i) (* (aref *pi-sl* i) xni t1)))
 249   (do ((jj 1 (1+ jj)))
 250       ((> jj l1))
 251     (cond ((> (cmod-sl (aref *hr-sl* *n*) (aref *hi-sl* *n*))
 252               (* 10.0 *are* (cmod-sl (aref *pr-sl* *n*) (aref *pi-sl* *n*))))
 253            (cdivid-sl (- (aref *pr-sl* *nn*)) (- (aref *pi-sl* *nn*)) (aref *hr-sl* *n*) (aref *hi-sl* *n*))
 254            (setq *tr* *cr* *ti* *ci*)
 255            (do ((j *n* (1- j)) (t1) (t2))
 256                ((> 1 j))
 257              (setq t1 (aref *hr-sl* (1- j)) t2 (aref *hi-sl* (1- j)))
 258              (setf (aref *hr-sl* j) (- (+ (aref *pr-sl* j) (* t1 *tr*)) (* t2 *ti*)))
 259              (setf (aref *hi-sl* j) (+ (aref *pi-sl* j) (* t1 *ti*) (* t2 *tr*))))
 260            (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
 261            (setf (aref *hi-sl* 0) (aref *pi-sl* 0)))
 262           (t (do ((j *n* (1- j)))
 263                  ((> 1 j))
 264                (setf (aref *hr-sl* j) (aref *hr-sl* (1- j)))
 265                (setf (aref *hi-sl* j) (aref *hi-sl* (1- j))))
 266              (setf (aref *hr-sl* 0) 0.0)
 267              (setf (aref *hi-sl* 0) 0.0)))))
 268
 269 (defun fxshft-sl (l2)
 270   (let ((test t)
 271         (pasd nil)
 272         (otr 0.0) (oti 0.0)
 273         (svsr 0.0) (svsi 0.0)
 274         (*bool* nil)
 275         (*pvr* 0.0) (*pvi* 0.0))
 276     (polyev-sl)
 277     (setq *conv* nil)
 278     (calct-sl)
 279     (do ((j 1 (1+ j)))
 280         ((> j l2))
 281       (setq otr *tr* oti *ti*)
 282       (nexth-sl)
 283       (calct-sl)
 284       (setq *zr* (+ *sr* *tr*) *zi* (+ *si* *ti*))
 285       (cond ((and (not *bool*) test (not (= j l2)))
 286              (cond ((> (* 0.5 (cmod-sl *zr* *zi*))
 287                        (cmod-sl (- *tr* otr) (- *ti* oti)))
 288                     (cond (pasd (do ((i 0 (1+ i)))
 289                                     ((> i *n*))
 290                                   (setf (aref *shr-sl* i) (aref *hr-sl* i))
 291                                   (setf (aref *shi-sl* i) (aref *hi-sl* i)))
 292                                 (setq svsr *sr* svsi *si*)
 293                                 (vrshft-sl 10.)
 294                                 (when *conv* (return nil))
 295                                 (setq test nil)
 296                                 (do ((i 0 (1+ i)))
 297                                     ((> i *n*))
 298                                   (setf (aref *hr-sl* i) (aref *shr-sl* i))
 299                                   (setf (aref *hi-sl* i) (aref *shi-sl* i)))
 300                                 (setq *sr* svsr *si* svsi)
 301                                 (polyev-sl)
 302                                 (calct-sl))
 303                           ((setq pasd t))))
 304                    ((setq pasd nil))))))
 305     (or *conv* (vrshft-sl 10))
 306     nil))
 307
 308 (defun vrshft-sl (l3)
 309   (setq *conv* nil *sr* *zr* *si* *zi*)
 310   (do ((i 1 (1+ i))
 311        (bool1 nil)
 312        (mp) (ms) (omp) (relstp) (tp) (r1))
 313       ((> i l3))
 314     (polyev-sl)
 315     (setq mp (cmod-sl *pvr* *pvi*) ms (cmod-sl *sr* *si*))
 316     (cond ((> (* 20.0 (errev-sl ms mp)) mp)
 317            (setq *conv* t *zr* *sr* *zi* *si*)
 318            (return t)))
 319     (cond ((= i 1) (setq omp mp))
 320           ((or bool1 (> omp mp) (not (< relstp 0.05)))
 321            (if (> (* 0.1 mp) omp)
 322                (return t)
 323                (setq omp mp)))
 324           (t (setq tp relstp bool1 t)
 325              (when (> *are* relstp) (setq tp *are*))
 326              (setq r1 (sqrt tp)
 327                    *sr* (prog1
 328                           (- (* (1+ r1) *sr*) (* r1 *si*))
 329                         (setq *si* (+ (* (1+ r1) *si*) (* r1 *sr*)))))
 330              (polyev-sl)
 331              (do ((j 1 (1+ j))) ((> j 5)) (calct-sl) (nexth-sl))
 332              (setq omp *infin*)))
 333     (calct-sl)
 334     (nexth-sl)
 335     (calct-sl)
 336     (or *bool*
 337         (setq relstp (/ (cmod-sl *tr* *ti*) (cmod-sl *sr* *si*))
 338               *sr* (+ *sr* *tr*)
 339               *si* (+ *si* *ti*)))))
 340
 341 (defun calct-sl nil
 342   (do ((i 1 (1+ i))
 343        ($t)
 344        (hvr (setf (aref *qhr-sl* 0) (aref *hr-sl* 0)))
 345        (hvi (setf (aref *qhi-sl* 0) (aref *hi-sl* 0))))
 346       ((> i *n*)
 347        (setq *bool* (not (> (cmod-sl hvr hvi) (* 10.0 *are* (cmod-sl (aref *hr-sl* *n*) (aref *hi-sl* *n*))))))
 348        (cond ((not *bool*) (cdivid-sl (- *pvr*) (- *pvi*) hvr hvi) (setq *tr* *cr* *ti* *ci*))
 349              (t (setq *tr* 0.0 *ti* 0.0)))
 350        nil)
 351     (setq $t (- (+ (aref *hr-sl* i) (* hvr *sr*)) (* hvi *si*)))
 352     (setf (aref *qhi-sl* i) (setq hvi (+ (aref *hi-sl* i) (* hvr *si*) (* hvi *sr*))))
 353     (setf (aref *qhr-sl* i) (setq hvr $t))))
 354
 355 (defun nexth-sl ()
 356   (cond (*bool* (do ((j 1 (1+ j)))
 357                   ((> j *n*))
 358                 (setf (aref *hr-sl* j) (aref *qhr-sl* (1- j)))
 359                 (setf (aref *hi-sl* j) (aref *qhi-sl* (1- j))))
 360               (setf (aref *hr-sl* 0) 0.0)
 361               (setf (aref *hi-sl* 0) 0.0))
 362         (t (do ((j 1. (1+ j)) (t1) (t2))
 363                ((> j *n*))
 364              (setq t1 (aref *qhr-sl* (1- j)) t2 (aref *qhi-sl* (1- j)))
 365              (setf (aref *hr-sl* j) (- (+ (aref *qpr-sl* j) (* t1 *tr*)) (* t2 *ti*)))
 366              (setf (aref *hi-sl* j) (+ (aref *qpi-sl* j) (* t1 *ti*) (* t2 *tr*))))
 367            (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
 368            (setf (aref *hi-sl* 0) (aref *qpi-sl* 0))))
 369   nil)
 370
 371 (defun polyev-sl ()
 372   (setq *pvr* (setf (aref *qpr-sl* 0) (aref *pr-sl* 0))
 373         *pvi* (setf (aref *qpi-sl* 0) (aref *pi-sl* 0)))
 374   (do ((i 1 (1+ i)) ($t))
 375       ((> i *nn*))
 376     (setq $t (- (+ (aref *pr-sl* i) (* *pvr* *sr*)) (* *pvi* *si*)))
 377     (setf (aref *qpi-sl* i) (setq *pvi* (+ (aref *pi-sl* i) (* *pvr* *si*) (* *pvi* *sr*))))
 378     (setf (aref *qpr-sl* i) (setq *pvr* $t))))
 379
 380 (defun errev-sl (ms mp)
 381   (- (* (do ((j 0 (1+ j))
 382                (*e* (/ (* (cmod-sl (aref *qpr-sl* 0) (aref *qpi-sl* 0)) *mre*) (+ *are* *mre*))))
 383               ((> j *nn*) *e*)
 384             (setq *e* (+ (cmod-sl (aref *qpr-sl* j) (aref *qpi-sl* j)) (* *e* ms))))
 385           (+ *are* *mre*))
 386       (* mp *mre*)))
 387
 388 ;; Compute a lower bound on the roots of the polynomial.  Let the
 389 ;; polynomial be sum(a[k]*x^k,k,0,n).  Then the lower bound is the
 390 ;; smallest real root of sum(|a[k]|*x^k,k,0,n-1)-a[n].  For our
 391 ;; purposes, this lower bound is computed to an accuracy of .005 or
 392 ;; so.
 393 (defun cauchy-sl ()
 394   (let ((x (exp (/ (- (log (aref *shr-sl* *nn*)) (log (aref *shr-sl* 0))) (float *nn*))))
 395         (xm 0.0))
 396     (setf (aref *shr-sl* *nn*) (- (aref *shr-sl* *nn*)))
 397     (cond ((not (zerop (aref *shr-sl* *n*)))
 398            (setq xm (- (/ (aref *shr-sl* *nn*) (aref *shr-sl* *n*))))
 399            (cond ((> x xm) (setq x xm)))))
 400     (do ((*f*))
 401         (nil)
 402       (setq xm (* 0.1 x) *f* (aref *shr-sl* 0))
 403       (do ((i 1 (1+ i))) ((> i *nn*)) (setq *f* (+ (aref *shr-sl* i) (* *f* xm))))
 404       (cond ((not (< 0.0 *f*)) (return t)))
 405       (setq x xm))
 406     (do ((dx x) (df) (*f*))
 407         ((> 5e-3 (abs (/ dx x))) x)
 408       (setq *f* (aref *shr-sl* 0) df *f*)
 409       (do ((i 1 (1+ i)))
 410           ((> i *n*))
 411         (setq *f* (+ (* *f* x) (aref *shr-sl* i)) df (+ (* df x) *f*)))
 412       (setq *f* (+ (* *f* x) (aref *shr-sl* *nn*)) dx (/ *f* df) x (- x dx)))))
 413
 414 ;; Scale the coefficients of the polynomial to reduce the chance of
 415 ;; overflow or underflow.  This is different from the algorithm given
 416 ;; in TOMS 419 and 493 which just scales the coefficients by some
 417 ;; fixed factor.  Instead, this algorithm computes a scale factor for
 418 ;; the roots and adjusts the coefficients of the polynomial
 419 ;; appropriately.
 420 ;;
 421 ;; The scale factor for the roots is in *polysc1*.  The roots are
 422 ;; scaled by 2^(*polysc1*).  There is an additional scaling *polysc*
 423 ;; such that the coefficients of the polynomial are all scaled by
 424 ;; 2^(*polysc*).
 425 (defun scale-sl ()
 426   (do ((i 0 (1+ i)) (j 0) (x 0.0) (dx 0.0))
 427       ((> i *nn*)
 428        (setq x (/ x (float (- (1+ *nn*) j)))
 429              dx (/ (- (log (aref *shr-sl* *nn*)) (log (aref *shr-sl* 0))) (float *nn*))
 430              *polysc1* (floor (+ 0.5 (/ dx *logbas*)))
 431              x (+ x (* (float (* *polysc1* *nn*)) *logbas* 0.5))
 432              *polysc* (floor (+ 0.5 (/ x *logbas*)))))
 433     (cond ((zerop (aref *shr-sl* i)) (setq j (1+ j)))
 434           (t (setq x (+ x (log (aref *shr-sl* i)))))))
 435   (do ((i *nn* (1- i)) (j (- *polysc*) (+ j *polysc1*)))
 436       ((< i 0))
 437     (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))
 438     (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) j))))
 439
 440 (defun cdivid-sl (ar ai br bi)
 441   (let ((r1 0.0))
 442     (cond ((and (zerop br) (zerop bi))
 443            (setq *cr* (setq *ci* *infin*)))
 444           ((> (abs bi) (abs br))
 445            (setq r1 (/ br bi)
 446                  bi (+ bi (* br r1))
 447                  br (+ ai (* ar r1))
 448                  *cr* (/ br bi)
 449                  br (- (* ai r1) ar)
 450                  *ci* (/ br bi)))
 451           ((setq r1 (/ bi br)
 452                  bi (+ br (* bi r1))
 453                  br (+ ar (* ai r1))
 454                  *cr* (/ br bi)
 455                  br (- ai (* ar r1))
 456                  *ci* (/ br bi)))))
 457   nil)
 458
 459 (defun cmod-sl (ar ai)
 460   (setq ar (abs ar)
 461         ai (abs ai))
 462   (cond ((> ai ar) (setq ar (/ ar ai)) (* ai (sqrt (1+ (* ar ar)))))
 463         ((> ar ai) (setq ai (/ ai ar)) (* ar (sqrt (1+ (* ai ai)))))
 464         ((* (sqrt 2.0)
 465           ar))))
 466
 467 ;;; This is the algorithm for doing real polynomials.  It is algorithm
 468 ;;; 493 from acm toms vol 1 p 178 (1975) by jenkins.  Note that array
 469 ;;; indexing starts from 0.  The names of the arrays have been changed
 470 ;;; to be the same as for cpoly.  The correspondence is: p - pr-sl, qp
 471 ;;; - qpr-sl, k - hr-sl, qk - qhr-sl, svk - shr-sl, temp - shi-sl.
 472 ;;; the roots *are* put in pr-sl and pi-sl.  The variable *si* appears
 473 ;;; not to be used here
 474
 475 (defun rpoly-sl (degree)
 476   (let ((*logbas* (log 2.0))
 477         (*infin* most-positive-flonum)
 478         (*are* flonum-epsilon)
 479         (*mre* 0.0)
 480         (xx (sqrt 0.5)) ;; sqrt(0.5)
 481         (yy 0.0)
 482         (cosr (cos (float (* 94/180 pi))))
 483         (sinr (sin (float (* 94/180 pi))))
 484         (aa 0.0)
 485         (cc 0.0)
 486         (bb 0.0)
 487         (bnd 0.0)
 488         (*sr* 0.0)
 489         (*u* 0.0)
 490         (*v* 0.0)
 491         (t1 0.0)
 492         (*szr* 0.0) (*szi* 0.0)
 493         (*lzr* 0.0) (*lzi* 0.0)
 494         (*nz* 0)
 495         (*n* 0)
 496         (*polysc* 0)
 497         (*polysc1* 0)
 498         (zerok 0)
 499         (conv1 t))
 500     (setq *mre* *are* yy (- xx))
 501     (do ((i degree (1- i))) ((not (zerop (aref *pr-sl* i))) (setq *nn* i *n* (1- i))))
 502     (setq degree *nn*)
 503     (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i))))
 504     (if (> *nn* 0) (scale-sl))
 505     ;; Loop to find all roots
 506     (do nil
 507         ((< *nn* 3)
 508          (cond ((= *nn* 2)
 509                 ;; Solve the final quadratic polynomial
 510                 (quad-sl (aref *pr-sl* 0.) (aref *pr-sl* 1) (aref *pr-sl* 2))
 511                 (cond ((and $polyfactor (not (zerop *szi*)))
 512                        (setf (aref *pr-sl* 2) (/ (aref *pr-sl* 2) (aref *pr-sl* 0)))
 513                        (setf (aref *pr-sl* 1) (/ (aref *pr-sl* 1) (aref *pr-sl* 0)))
 514                        (setf (aref *pi-sl* 2) 1.0))
 515                       (t (setf (aref *pr-sl* 2) *szr*)
 516                          (setf (aref *pi-sl* 2) *szi*)
 517                          (setf (aref *pr-sl* 1) *lzr*)
 518                          (setf (aref *pi-sl* 1) *lzi*))))
 519                (t
 520                 ;; Solve the final linear polynomial
 521                 (setf (aref *pr-sl* 1) (- (/ (aref *pr-sl* 1) (aref *pr-sl* 0))))))
 522          (setq *nn* 0))
 523       ;; Calculate a lower bound on the modulus of the zeros.
 524       (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i))))
 525       (setq bnd (cauchy-sl))
 526
 527       ;; Compute derivative of polynomial.  Save result in *hr-sl*.
 528       (do ((i 1 (1+ i)))
 529           ((> i *n*))
 530         (setf (aref *hr-sl* i) (/ (* (float (- *n* i)) (aref *pr-sl* i)) (float *n*))))
 531       (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
 532       (setq aa (aref *pr-sl* *nn*) bb (aref *pr-sl* *n*) zerok (zerop (aref *hr-sl* *n*)))
 533       ;; Do 5 steps with no shift.
 534       (do ((jj 1 (1+ jj)))
 535           ((> jj 5.))
 536         (setq cc (aref *hr-sl* *n*))
 537         (cond (zerok (do ((j *n* (1- j)))
 538                          ((< j 1))
 539                        (setf (aref *hr-sl* j) (aref *hr-sl* (1- j))))
 540                      (setf (aref *hr-sl* 0) 0.0)
 541                      (setq zerok (zerop (aref *hr-sl* *n*))))
 542               (t (setq t1 (- (/ aa cc)))
 543                  (do ((j *n* (1- j)))
 544                      ((< j 1))
 545                    (setf (aref *hr-sl* j) (+ (* t1 (aref *hr-sl* (1- j))) (aref *pr-sl* j))))
 546                  (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
 547                  (setq zerok (not (> (abs (aref *hr-sl* *n*))
 548                                      (* (abs bb) *are* 10.0)))))))
 549       (do ((i 0 (1+ i))) ((> i *n*)) (setf (aref *shi-sl* i) (aref *hr-sl* i)))
 550       (do ((cnt 1 (1+ cnt)))
 551           ((> cnt 20.) (setq conv1 nil))
 552         (setq xx (prog1
 553                      (- (* cosr xx) (* sinr yy))
 554                    (setq yy (+ (* sinr xx) (* cosr yy))))
 555               *sr* (* bnd xx)
 556               *u* (* -2.0 *sr*)
 557               *v* bnd)
 558         (fxshfr-sl (* 20 cnt))
 559         (cond ((> *nz* 0)
 560                (setf (aref *pr-sl* *nn*) *szr*)
 561                (setf (aref *pi-sl* *nn*) *szi*)
 562                (cond ((= *nz* 2)
 563                       (setf (aref *pr-sl* *n*) *lzr*)
 564                       (setf (aref *pi-sl* *n*) *lzi*)
 565                       (cond ((and $polyfactor (not (zerop *szi*)))
 566                              (setf (aref *pr-sl* *nn*) *v*)
 567                              (setf (aref *pr-sl* *n*) *u*)
 568                              (setf (aref *pi-sl* *nn*) 1.0)))))
 569                (setq *nn* (- *nn* *nz*) *n* (1- *nn*))
 570                (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *pr-sl* i) (aref *qpr-sl* i)))
 571                (return nil)))
 572         (do ((i 0 (1+ i))) ((> i *n*)) (setf (aref *hr-sl* i) (aref *shi-sl* i))))
 573       (or conv1 (return nil)))
 574     (cond ($polyfactor
 575            (do ((i degree (1- i)))
 576                ((= i *nn*))
 577              (cond ((zerop (aref *pi-sl* i))
 578                     (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*)))
 579                    (t (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) (* 2 *polysc1*)))
 580                       (setq i (1- i))
 581                       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*))))))
 582           (t (do ((i (1+ *nn*) (1+ i)))
 583                  ((> i degree))
 584                (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*))
 585                (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) *polysc1*)))))
 586     (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
 587         ((> i *nn*))
 588       (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j)))))
 589
 590 (defun fxshfr-sl (l2)
 591   (let ((*my-type* 0)
 592         (*a* 0.0) (*b* 0.0) (*c* 0.0) (*d* 0.0) (*e* 0.0) (*f* 0.0) (*g* 0.0) (*h* 0.0)
 593         (*a1* 0.0) (*a3* 0.0) (*a7* 0.0))
 594     (declare (special *my-type*))
 595     (setq *nz* 0)
 596     (quadsd-sl)
 597     (calcsc-sl)
 598     (do ((j 1 (1+ j))
 599          (betav 0.25)
 600          (betas 0.25)
 601          (oss *sr*)
 602          (ovv *v*)
 603          (tvv) (tss) (ss) (vv) (tv) (ts) (ots) (otv)
 604          (*ui*) (*vi*) (*s*) (svv) (svu) (iflag) (vpass) (spass) (vtry) (stry))
 605         ((> j l2))
 606       (nextk-sl)
 607       (calcsc-sl)
 608       (newest-sl)
 609       (setq vv *vi*
 610             ss 0.0)
 611       (or (zerop (aref *hr-sl* *n*))
 612           (setq ss (- (/ (aref *pr-sl* *nn*) (aref *hr-sl* *n*)))))
 613       (setq tv 1.0 ts 1.0)
 614       (cond ((not (or (= j 1) (= *my-type* 3)))
 615              (or (zerop vv) (setq tv (abs (/ (- vv ovv) vv))))
 616              (or (zerop ss) (setq ts (abs (/ (- ss oss) ss))))
 617              (setq tvv 1.0)
 618              (and (< tv otv) (setq tvv (* tv otv)))
 619              (setq tss 1.0)
 620              (and (< ts ots) (setq tss (* ts ots)))
 621              (setq vpass (< tvv betav) spass (< tss betas))
 622              (cond ((or spass vpass)
 623                     (setq svu *u* svv *v*)
 624                     (do ((i 0 (1+ i)))
 625                         ((> i *n*)) (setf (aref *shr-sl* i)
 626                                           (aref *hr-sl* i)))
 627                     (setq *s* ss vtry nil stry nil)
 628                     (and (do ((*bool* (not (and spass (or (not vpass) (< tss tvv)))) t)
 629                               (l50 nil nil))
 630                              (nil)
 631                            (cond (*bool* (quadit-sl)
 632                                          (and (> *nz* 0) (return t))
 633                                          (setq vtry t
 634                                                betav (* 0.25 betav))
 635                                          (cond ((or stry (not spass))
 636                                                 (setq l50 t))
 637                                                (t (do ((i 0 (1+ i)))
 638                                                       ((> i *n*))
 639                                                     (setf (aref *hr-sl* i)
 640                                                           (aref *shr-sl* i)))))))
 641                            (cond ((not l50)
 642                                   (setq iflag (realit-sl))
 643                                   (and (> *nz* 0) (return t))
 644                                   (setq stry t betas (* 0.25 betas))
 645                                   (cond ((zerop iflag) (setq l50 t))
 646                                         (t (setq *ui* (- (+ *s* *s*))
 647                                                  *vi* (* *s* *s*))))))
 648                            (cond (l50 (setq *u* svu *v* svv)
 649                                       (do ((i 0 (1+ i)))
 650                                           ((> i *n*))
 651                                         (setf (aref *hr-sl* i)
 652                                               (aref *shr-sl* i)))
 653                                       (and (or (not vpass) vtry)
 654                                            (return nil)))))
 655                          (return nil))
 656                     (quadsd-sl)
 657                     (calcsc-sl)))))
 658       (setq ovv vv
 659             oss ss
 660             otv tv
 661             ots ts))))
 662
 663 (defun quadit-sl nil
 664   (setq *nz* 0 *u* *ui* *v* *vi*)
 665   (do ((tried) (j 0) (ee) (mp) (relstp) (omp) (ms))
 666       (nil)
 667     (quad-sl 1.0 *u* *v*)
 668     (and (> (abs (- (abs *szr*) (abs *lzr*))) (* 1e-2 (abs *lzr*)))
 669          (return nil))
 670     (quadsd-sl)
 671     (setf mp (+ (abs (- *a* (* *szr* *b*))) (abs (* *szi* *b*))))
 672     (setf ms (cmod-sl *szr* *szi*))
 673     (setf ee (errev-sl ms mp))
 674     (cond ((not (> mp (* 2e1 ee))) (setq *nz* 2)
 675            (return nil)))
 676     (setq j (1+ j))
 677     (and (> j 20) (return nil))
 678     (cond ((not (or (< j 2) (> relstp 1e-2) (< mp omp) tried))
 679            (and (< relstp *are*) (setq relstp *are*))
 680            (setq relstp (sqrt relstp)
 681                  *u* (- *u* (* *u* relstp))
 682                  *v* (+ *v* (* *v* relstp)))
 683            (quadsd-sl)
 684            (do ((i 1 (1+ i)))
 685                ((> i 5)) (calcsc-sl) (nextk-sl))
 686            (setq tried t j 0)))
 687     (setq omp mp)
 688     (calcsc-sl)
 689     (nextk-sl)
 690     (calcsc-sl)
 691     (newest-sl)
 692     (and (zerop *vi*) (return nil))
 693     (setq relstp (abs (/ (- *vi* *v*) *vi*)) *u* *ui* *v* *vi*)))
 694
 695 (defun realit-sl ()
 696   (setq *nz* 0)
 697   (do ((j 0) (pv) (ee) (ms) (mp) (kv) (t1) (omp))
 698       (nil)
 699     (setq pv (aref *pr-sl* 0))
 700     (setf (aref *qpr-sl* 0) pv)
 701     (do ((i 1 (1+ i)))
 702         ((> i *nn*))
 703       (setq pv (+ (* pv *s*) (aref *pr-sl* i)))
 704       (setf (aref *qpr-sl* i) pv))
 705     (setq mp (abs pv)
 706           ms (abs *s*)
 707           ee (* (/ *mre* (+ *are* *mre*)) (abs (aref *qpr-sl* 0))))
 708     (do ((i 1 (1+ i)))
 709         ((> i *nn*)) (setq ee (+ (* ee ms) (abs (aref *qpr-sl* i)))))
 710     (cond ((not (> mp (* 2e1 (- (* (+ *are* *mre*) ee) (* *mre* mp)))))
 711            (setq *nz* 1 *szr* *s* *szi* 0.0)
 712            (return 0)))
 713     (setq j (1+ j))
 714     (and (> j 10) (return 0))
 715     (cond ((not (or (< j 2)
 716                     (> (abs t1) (* 1e-3 (abs (- *s* t1))))
 717                     (not (> mp omp))))
 718            (return 1)))
 719     (setq omp mp kv (aref *hr-sl* 0))
 720     (setf (aref *qhr-sl* 0) kv)
 721     (do ((i 1 (1+ i)))
 722         ((> i *n*))
 723       (setq kv (+ (* kv *s*) (aref *hr-sl* i)))
 724       (setf (aref *qhr-sl* i) kv))
 725     (cond ((> (abs kv) (* (abs (aref *hr-sl* *n*)) 1e1 *are*))
 726            (setq t1 (- (/ pv kv)))
 727            (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
 728            (do ((i 1 (1+ i)))
 729                ((> i *n*))
 730              (setf (aref *hr-sl* i)
 731                     (+ (* t1 (aref *qhr-sl* (1- i))) (aref *qpr-sl* i)))))
 732           (t (setf (aref *hr-sl* 0) 0.0)
 733              (do ((i 1 (1+ i)))
 734                  ((> i *n*)) (setf (aref *hr-sl* i) (aref *qhr-sl* (1- i))))))
 735     (setq kv (aref *hr-sl* 0))
 736     (do ((i 1 (1+ i)))
 737         ((> i *n*)) (setq kv (+ (* kv *s*) (aref *hr-sl* i))))
 738     (setq t1 0.0)
 739     (and (> (abs kv) (* (abs (aref *hr-sl* *n*)) 10.0 *are*))
 740          (setq t1 (- (/ pv kv))))
 741     (setq *s* (+ *s* t1))))
 742
 743 (defun calcsc-sl ()
 744   (declare (special *my-type*))
 745   (setq *d* (aref *hr-sl* 0))
 746   (setf (aref *qhr-sl* 0) *d*)
 747   (setq *c* (- (aref *hr-sl* 1) (* *u* *d*)))
 748   (setf (aref *qhr-sl* 1) *c*)
 749   (do ((i 2 (1+ i))
 750        (c0))
 751       ((> i *n*))
 752     (setq c0 (- (aref *hr-sl* i) (* *u* *c*) (* *v* *d*)))
 753     (setf (aref *qhr-sl* i) c0)
 754     (setq *d* *c* *c* c0))
 755   (cond ((not (or (> (abs *c*) (* (abs (aref *hr-sl* *n*)) 1e2 *are*))
 756                   (> (abs *d*) (* (abs (aref *hr-sl* (1- *n*))) 1e2 *are*))))
 757          (setq *my-type* 3))
 758         ((not (< (abs *d*) (abs *c*)))
 759          (setq *my-type* 2
 760                *e* (/ *a* *d*)
 761                *f* (/ *c* *d*)
 762                *g* (* *u* *b*)
 763                *h* (* *v* *b*)
 764                *a3* (+ (* (+ *a* *g*) *e*) (* *h* (/ *b* *d*)))
 765                *a1* (- (* *b* *f*) *a*)
 766                *a7* (+ (* (+ *f* *u*) *a*) *h*)))
 767         (t (setq *my-type* 1
 768                  *e* (/ *a* *c*)
 769                  *f* (/ *d* *c*)
 770                  *g* (* *u* *e*)
 771                  *h* (* *v* *b*)
 772                  *a3* (+ (* *a* *e*) (* (+ (/ *h* *c*) *g*) *b*))
 773                  *a1* (- *b* (* *a* (/ *d* *c*)))
 774                  *a7* (+ *a* (* *g* *d*) (* *h* *f*)))))
 775   nil)
 776
 777 (defun nextk-sl ()
 778   (declare (special *my-type*))
 779   (cond ((= *my-type* 3)
 780          (setf (aref *hr-sl* 0) 0.0)
 781          (setf (aref *hr-sl* 1) 0.0)
 782          (do ((i 2 (1+ i)))
 783              ((> i *n*)) (setf (aref *hr-sl* i) (aref *qhr-sl* (- i 2)))))
 784         ((> (abs *a1*) (* (abs (if (= *my-type* 1) *b* *a*)) 1e1 *are*))
 785          (setq *a7* (/ *a7* *a1*) *a3* (/ *a3* *a1*))
 786          (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
 787          (setf (aref *hr-sl* 1) (- (aref *qpr-sl* 1) (* *a7* (aref *qpr-sl* 0))))
 788          (do ((i 2 (1+ i)))
 789              ((> i *n*))
 790            (setf (aref *hr-sl* i)
 791                   (+ (* *a3* (aref *qhr-sl* (- i 2)))
 792                       (- (* *a7* (aref *qpr-sl* (1- i))))
 793                       (aref *qpr-sl* i)))))
 794         (t (setf (aref *hr-sl* 0) 0.0)
 795            (setf (aref *hr-sl* 1) (- (* *a7* (aref *qpr-sl* 0))))
 796            (do ((i 2 (1+ i)))
 797                ((> i *n*))
 798              (setf (aref *hr-sl* i)
 799                     (- (* *a3* (aref *qhr-sl* (- i 2)))
 800                         (* *a7* (aref *qpr-sl* (1- i))))))))
 801   nil)
 802
 803 (defun newest-sl ()
 804   (declare (special *my-type*))
 805   (let ((a4 0.0) (a5 0.0)
 806         (b1 0.0) (b2 0.0)
 807         (c1 0.0) (c2 0.0) (c3 0.0) (c4 0.0))
 808     (cond ((= *my-type* 3)
 809            (setq *ui* 0.0 *vi* 0.0))
 810           (t
 811            (if (= *my-type* 2)
 812                (setq a4 (+ (* (+ *a* *g*) *f*) *h*)
 813                      a5 (+ (* (+ *f* *u*) *c*) (* *v* *d*)))
 814                (setq a4 (+ *a* (* *u* *b*) (* *h* *f*))
 815                      a5 (+ *c* (* (+ *u* (* *v* *f*)) *d*))))
 816              (setq b1 (- (/ (aref *hr-sl* *n*) (aref *pr-sl* *nn*)))
 817                    b2 (- (/ (+ (aref *hr-sl* (1- *n*)) (* b1 (aref *pr-sl* *n*))) (aref *pr-sl* *nn*)))
 818                    c1 (* *v* b2 *a1*)
 819                    c2 (* b1 *a7*)
 820                    c3 (* b1 b1 *a3*)
 821                    c4 (- c1 c2 c3)
 822                    c1 (+ a5 (* b1 a4) (- c4)))
 823              (if (zerop c1)
 824                  (setq *ui* 0.0 *vi* 0.0)
 825                  (setq *ui* (- *u* (/ (+ (* *u* (+ c3 c2))
 826                                          (* *v* (+ (* b1 *a1*) (* b2 *a7*))))
 827                                       c1))
 828                        *vi* (* *v* (1+ (/ c4 c1)))))))
 829     nil))
 830
 831 ;; Divide the polynomial in *pr-sl* by the quadratic x^2 + (*u*)*x +
 832 ;; (*v*).  Place the quotient in *qpr-sl* and the remainder in *a* and
 833 ;; *b*.  I (rtoy) think the remainder polynomial is (*b*)*x + (*a*).
 834 (defun quadsd-sl ()
 835   (setq *b* (aref *pr-sl* 0))
 836   (setf (aref *qpr-sl* 0) *b*)
 837   (setq *a* (- (aref *pr-sl* 1) (* *u* *b*)))
 838   (setf (aref *qpr-sl* 1) *a*)
 839   (do ((i 2 (1+ i))
 840        (c0))
 841       ((> i *nn*))
 842     (setq c0 (- (aref *pr-sl* i) (* *u* *a*) (* *v* *b*)))
 843     (setf (aref *qpr-sl* i) c0)
 844     (setq *b* *a*
 845           *a* c0)))
 846
 847 ;; Compute the zeros of the quadratic a0*z^2+b1*z+c0.  The larger zero
 848 ;; is returned in *szr* and *szi*.  The smaller zero is in *lzr* and
 849 ;; *lzi*.
 850 (defun quad-sl (a0 b1 c0)
 851   (setq *szr* 0.0 *szi* 0.0 *lzr* 0.0 *lzi* 0.0)
 852   (let ((b0 0.0)
 853         (l0 0.0)
 854         (*e* 0.0))
 855     ;; Handle the degenerate cases of a0 = 0 or c0 = 0 first.
 856     (cond ((zerop a0) (or (zerop b1) (setq *szr* (- (/ c0 b1)))))
 857           ((zerop c0) (setq *lzr* (- (/ b1 a0))))
 858           (t
 859            ;; Quadratic formula.
 860            (setq b0 (/ b1 2.0))
 861            (cond ((< (abs b0) (abs c0))
 862                   (setq *e* a0)
 863                   (and (< c0 0.0) (setq *e* (- a0)))
 864                   (setq *e* (- (* b0 (/ b0 (abs c0))) *e*)
 865                         l0 (* (sqrt (abs *e*)) (sqrt (abs c0)))))
 866                  (t (setq *e* (- 1.0 (* (/ a0 b0) (/ c0 b0)))
 867                           l0 (* (sqrt (abs *e*)) (abs b0)))))
 868            (cond ((< *e* 0.0)
 869                   (setq *szr* (- (/ b0 a0))
 870                         *lzr* *szr*
 871                         *szi* (abs (/ l0 a0))
 872                         *lzi* (- *szi*)))
 873                  (t (or (< b0 0.0) (setq l0 (- l0)))
 874                     (setq *lzr* (/ (- l0 b0) a0))
 875                     (or (zerop *lzr*) (setq *szr* (/ c0 *lzr* a0)))))))
 876     nil))
 877 \f
 878 ;; This is a very straightforward conversion of $allroots to use
 879 ;; bfloats instead of floats.
 880
 881 (defun bf-cpoly-err (expr)
 882   (merror (intl:gettext "bfallroots: expected a polynomial; found ~M") expr))
 883
 884 (defun fpzerop (x)
 885   (equal '(0 0) x))
 886
 887 ;; (ar+%i*ai)/(br+%i*bi) -> cr+%i*ci.
 888 (defun bf-cdivid-sl (ar ai br bi)
 889   (cond ((and (fpzerop br)
 890               (fpzerop bi))
 891          ;; Division by zero.  Should we do something else besides set
 892          ;; both parts to be "infinity"?
 893          (setq *cr* (setq *ci* *infin*)))
 894         ((fpgreaterp (fpabs bi) (fpabs br))
 895          (let ((r1 (fpquotient br bi)))
 896            (setq bi (fpplus bi (fptimes* br r1))
 897                  br (fpplus ai (fptimes* ar r1))
 898                  *cr* (fpquotient br bi)
 899                  br (fpdifference (fptimes* ai r1) ar)
 900                  *ci* (fpquotient br bi))))
 901         (t
 902          (let ((r1 (fpquotient bi br)))
 903            (setq bi (fpplus br (fptimes* bi r1))
 904                  br (fpplus ar (fptimes* ai r1))
 905                  *cr* (fpquotient br bi)
 906                  br (fpdifference ai (fptimes* ar r1))
 907                  *ci* (fpquotient br bi))))))
 908
 909 (defun fpsqrt (x)
 910   (fproot (bcons x) 2))
 911
 912 (defun bf-cmod-sl (ar ai)
 913   (let ((ar (fpabs ar))
 914         (ai (fpabs ai)))
 915     (cond ((fpgreaterp ai ar)
 916            (setq ar (fpquotient ar ai))
 917            (fptimes* ai
 918                      (fpsqrt (fpplus (fpone) (fptimes* ar ar)))))
 919           ((fpgreaterp ar ai)
 920            (setq ai (fpquotient ai ar))
 921            (fptimes* ar (fpsqrt (fpplus (fpone) (fptimes* ai ai)))))
 922           ((fptimes* (fpsqrt (intofp 2))
 923                      ar)))))
 924
 925
 926 (defun bf-calct-sl nil
 927   (do ((i 1 (1+ i))
 928        (tt)
 929        (hvr (setf (aref *qhr-sl* 0) (aref *hr-sl* 0)))
 930        (hvi (setf (aref *qhi-sl* 0) (aref *hi-sl* 0))))
 931       ((> i *n*)
 932        (setq *bool*
 933              (not (fpgreaterp (bf-cmod-sl hvr hvi)
 934                               (fptimes* (intofp 10)
 935                                         (fptimes* *are*
 936                                                   (bf-cmod-sl (aref *hr-sl* *n*)
 937                                                               (aref *hi-sl* *n*)))))))
 938        (cond ((not *bool*)
 939               (bf-cdivid-sl (fpminus *pvr*) (fpminus *pvi*) hvr hvi)
 940               (setq *tr* *cr*
 941                     *ti* *ci*))
 942              (t (setq *tr* (intofp 0) *ti* (intofp 0))))
 943        nil)
 944     (setq tt (fpdifference (fpplus (aref *hr-sl* i)
 945                                    (fptimes* hvr *sr*))
 946                            (fptimes* hvi *si*)))
 947     (setf (aref *qhi-sl* i)
 948           (setq hvi (fpplus (aref *hi-sl* i)
 949                             (fpplus (fptimes* hvr *si*)
 950                                     (fptimes* hvi *sr*)))))
 951     (setf (aref *qhr-sl* i) (setq hvr tt))))
 952
 953 (defun bf-nexth-sl ()
 954   (cond (*bool*
 955          (do ((j 1 (1+ j)))
 956              ((> j *n*))
 957            (setf (aref *hr-sl* j) (aref *qhr-sl* (1- j)))
 958            (setf (aref *hi-sl* j) (aref *qhi-sl* (1- j))))
 959          (setf (aref *hr-sl* 0) (intofp 0))
 960          (setf (aref *hi-sl* 0) (intofp 0)))
 961         (t
 962          (do ((j 1. (1+ j))
 963               (t1)
 964               (t2))
 965              ((> j *n*))
 966            (setq t1 (aref *qhr-sl* (1- j))
 967                  t2 (aref *qhi-sl* (1- j)))
 968            (setf (aref *hr-sl* j)
 969                  (fpdifference (fpplus (aref *qpr-sl* j)
 970                                        (fptimes* t1 *tr*))
 971                                (fptimes* t2 *ti*)))
 972            (setf (aref *hi-sl* j)
 973                  (fpplus (aref *qpi-sl* j)
 974                          (fpplus (fptimes* t1 *ti*)
 975                                  (fptimes* t2 *tr*)))))
 976          (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
 977          (setf (aref *hi-sl* 0) (aref *qpi-sl* 0))))
 978   nil)
 979
 980 (defun bf-polyev-sl ()
 981   (setq *pvr* (setf (aref *qpr-sl* 0) (aref *pr-sl* 0))
 982         *pvi* (setf (aref *qpi-sl* 0) (aref *pi-sl* 0)))
 983   (do ((i 1 (1+ i))
 984        (tt))
 985       ((> i *nn*))
 986     (setq tt (fpdifference (fpplus (aref *pr-sl* i) (fptimes* *pvr* *sr*))
 987                            (fptimes* *pvi* *si*)))
 988     (setf (aref *qpi-sl* i)
 989           (setq *pvi* (fpplus (aref *pi-sl* i)
 990                               (fpplus (fptimes* *pvr* *si*)
 991                                       (fptimes* *pvi* *sr*)))))
 992     (setf (aref *qpr-sl* i)
 993           (setq *pvr* tt))))
 994
 995 (defun bf-errev-sl (ms mp)
 996   (fpdifference
 997    (fptimes* (do ((j 0 (1+ j))
 998                   (e (fpquotient (fptimes* (bf-cmod-sl (aref *qpr-sl* 0) (aref *qpi-sl* 0))
 999                                            *mre*)
1000                                  (fpplus *are* *mre*))))
1001                  ((> j *nn*) e)
1002                (setq e (fpplus (bf-cmod-sl (aref *qpr-sl* j) (aref *qpi-sl* j))
1003                                (fptimes* e ms))))
1004              (fpplus *are* *mre*))
1005    (fptimes* mp *mre*)))
1006
1007 (defun bf-cauchy-sl ()
1008   (let ((x (fpexp (fpquotient (fpdifference (fplog (aref *shr-sl* *nn*))
1009                                             (fplog (aref *shr-sl* 0)))
1010                               (intofp *nn*))))
1011         (xm (intofp 0)))
1012     (setf (aref *shr-sl* *nn*) (fpminus (aref *shr-sl* *nn*)))
1013     (cond ((not (fpzerop (aref *shr-sl* *n*)))
1014            (setq xm (fpminus (fpquotient (aref *shr-sl* *nn*)
1015                                          (aref *shr-sl* *n*))))
1016            (cond ((fpgreaterp x xm) (setq x xm)))))
1017     (do ((f))
1018         (nil)
1019       (setq xm (fptimes* (intofp 0.1) x)
1020             f (aref *shr-sl* 0))
1021       (do ((i 1 (1+ i)))
1022           ((> i *nn*))
1023         (setq f (fpplus (aref *shr-sl* i)
1024                         (fptimes* f xm))))
1025       ;;(cond ((not (< 0.0 f)) (return t)))
1026       (when (fpgreaterp (intofp 0) f)
1027         (return t))
1028       (setq x xm))
1029     (do ((dx x)
1030          (df)
1031          (f))
1032         ((fpgreaterp (intofp 5e-3)
1033                      (fpabs (fpquotient dx x)))
1034          x)
1035       (setq f (aref *shr-sl* 0)
1036             df f)
1037       (do ((i 1 (1+ i)))
1038           ((> i *n*))
1039         (setq f (fpplus (fptimes* f x)
1040                         (aref *shr-sl* i))
1041               df (fpplus (fptimes* df x)
1042                          f)))
1043       (setq f (fpplus (fptimes* f x)
1044                       (aref *shr-sl* *nn*))
1045             dx (fpquotient f df)
1046             x (fpdifference x dx)))))
1047
1048 (defun bf-scale-float (bf scale)
1049   (destructuring-bind (mantissa exp)
1050       bf
1051     (if (zerop mantissa)
1052         (list mantissa 0)
1053         (list mantissa
1054               (+ exp scale)))))
1055
1056 (defun bf-scale-sl ()
1057   (do ((i 0 (1+ i))
1058        (j 0)
1059        (x (intofp 0))
1060        (dx (intofp 0)))
1061       ((> i *nn*)
1062        (setq x (fpquotient x (intofp (- (1+ *nn*) j)))
1063              dx (fpquotient (fpdifference (fplog (aref *shr-sl* *nn*))
1064                                           (fplog (aref *shr-sl* 0)))
1065                             (intofp *nn*))
1066              *polysc1* (fpentier (bcons (fpplus (cdr bfhalf)
1067                                                 (fpquotient dx *logbas*))))
1068              x (fpplus x (fptimes* (intofp (* *polysc1* *nn*))
1069                                    (fptimes* *logbas*
1070                                              (cdr bfhalf))))
1071              *polysc* (fpentier (bcons (fpplus (cdr bfhalf) (fpquotient x *logbas*))))))
1072     (cond ((equalp (aref *shr-sl* i) (intofp 0))
1073            (setq j (1+ j)))
1074           (t
1075            (setq x (fpplus x (fplog (aref *shr-sl* i)))))))
1076   (do ((i *nn* (1- i))
1077        (j (- *polysc*) (+ j *polysc1*)))
1078       ((< i 0))
1079     (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) j))
1080     (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) j)))
1081   nil)
1082
1083 (defun bf-noshft-sl (l1)
1084   (do ((i 0 (1+ i))
1085        (xni (intofp *nn*) (fpdifference xni (intofp 1)))
1086        (t1 (fpquotient (fpone) (intofp *nn*))))
1087       ((> i *n*))
1088     (setf (aref *hr-sl* i) (fptimes* (aref *pr-sl* i)
1089                                      (fptimes* xni t1)))
1090     (setf (aref *hi-sl* i) (fptimes* (aref *pi-sl* i)
1091                                      (fptimes* xni t1))))
1092   (do ((jj 1 (1+ jj)))
1093       ((> jj l1))
1094     (cond ((fpgreaterp (bf-cmod-sl (aref *hr-sl* *n*) (aref *hi-sl* *n*))
1095                        (fptimes* (intofp 10)
1096                                  (fptimes* *are*
1097                                            (bf-cmod-sl (aref *pr-sl* *n*)
1098                                                        (aref *pi-sl* *n*)))))
1099            (bf-cdivid-sl (fpminus (aref *pr-sl* *nn*))
1100                          (fpminus (aref *pi-sl* *nn*))
1101                          (aref *hr-sl* *n*)
1102                          (aref *hi-sl* *n*))
1103            (setq *tr* *cr*
1104                  *ti* *ci*)
1105            (do ((j *n* (1- j)) (t1) (t2))
1106                ((> 1 j))
1107              (setq t1 (aref *hr-sl* (1- j))
1108                    t2 (aref *hi-sl* (1- j)))
1109              (setf (aref *hr-sl* j) (fpdifference (fpplus (aref *pr-sl* j)
1110                                                           (fptimes* t1 *tr*))
1111                                                   (fptimes* t2 *ti*)))
1112              (setf (aref *hi-sl* j) (fpplus (aref *pi-sl* j)
1113                                             (fpplus (fptimes* t1 *ti*)
1114                                                     (fptimes* t2 *tr*)))))
1115            (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
1116            (setf (aref *hi-sl* 0) (aref *pi-sl* 0)))
1117           (t (do ((j *n* (1- j)))
1118                  ((> 1 j))
1119                (setf (aref *hr-sl* j) (aref *hr-sl* (1- j)))
1120                (setf (aref *hi-sl* j) (aref *hi-sl* (1- j))))
1121              (setf (aref *hr-sl* 0) (intofp 0))
1122              (setf (aref *hi-sl* 0) (intofp 0))))))
1123
1124 (defun bf-vrshft-sl (l3)
1125   (setq *conv* nil
1126         *sr* *zr*
1127         *si* *zi*)
1128   (do ((i 1 (1+ i))
1129        (bool1 nil)
1130        (mp)
1131        (ms)
1132        (omp)
1133        (relstp)
1134        (tp)
1135        (r1))
1136       ((> i l3))
1137     (bf-polyev-sl)
1138     (setq mp (bf-cmod-sl *pvr* *pvi*)
1139           ms (bf-cmod-sl *sr* *si*))
1140     (cond ((fpgreaterp (fptimes* (intofp 20) (bf-errev-sl ms mp))
1141                        mp)
1142            (setq *conv* t
1143                  *zr* *sr*
1144                  *zi* *si*)
1145            (return t)))
1146     (cond ((= i 1)
1147            (setq omp mp))
1148           ((or bool1
1149                (fpgreaterp omp mp)
1150                ;;(not (< relstp 0.05))
1151                (fpgreaterp relstp (cdr bfhalf)))
1152            (if (fpgreaterp (fptimes* (intofp 0.1) mp)
1153                            omp)
1154                (return t)
1155                (setq omp mp)))
1156           (t
1157            (setq tp relstp
1158                  bool1 t)
1159            (when (fpgreaterp *are* relstp)
1160              (setq tp *are*))
1161            (setq r1 (fpsqrt tp)
1162                  *sr* (prog1
1163                           (fpdifference (fptimes* (fpplus (fpone) r1)
1164                                                   *sr*)
1165                                         (fptimes* r1 *si*))
1166                         (setq *si* (fpplus (fptimes* (fpplus (fpone) r1)
1167                                                      *si*)
1168                                            (fptimes* r1 *sr*)))))
1169            (bf-polyev-sl)
1170            (do ((j 1 (1+ j)))
1171                ((> j 5))
1172              (bf-calct-sl)
1173              (bf-nexth-sl))
1174            (setq omp *infin*)))
1175     (bf-calct-sl)
1176     (bf-nexth-sl)
1177     (bf-calct-sl)
1178     (or *bool*
1179         (setq relstp (fpquotient (bf-cmod-sl *tr* *ti*)
1180                                  (bf-cmod-sl *sr* *si*))
1181               *sr* (fpplus *sr* *tr*)
1182               *si* (fpplus *si* *ti*)))))
1183
1184 (defun bf-fxshft-sl (l2)
1185   (let ((test t)
1186         (pasd nil)
1187         (otr (intofp 0))
1188         (oti (intofp 0))
1189         (svsr (intofp 0))
1190         (svsi (intofp 0))
1191         (*bool* nil)
1192         (*pvr* (intofp 0))
1193         (*pvi* (intofp 0)))
1194     (bf-polyev-sl)
1195     (setq *conv* nil)
1196     (bf-calct-sl)
1197     (do ((j 1 (1+ j)))
1198         ((> j l2))
1199       (setq otr *tr*
1200             oti *ti*)
1201       (bf-nexth-sl)
1202       (bf-calct-sl)
1203       (setq *zr* (fpplus *sr* *tr*)
1204             *zi* (fpplus *si* *ti*))
1205       (cond ((and (not *bool*)
1206                   test
1207                   (not (= j l2)))
1208              (cond ((fpgreaterp (fptimes* (cdr bfhalf) (bf-cmod-sl *zr* *zi*))
1209                                 (bf-cmod-sl (fpdifference *tr* otr)
1210                                             (fpdifference *ti* oti)))
1211                     (cond (pasd
1212                            (do ((i 0 (1+ i)))
1213                                ((> i *n*))
1214                              (setf (aref *shr-sl* i) (aref *hr-sl* i))
1215                              (setf (aref *shi-sl* i) (aref *hi-sl* i)))
1216                            (setq svsr *sr* svsi *si*)
1217                            (bf-vrshft-sl 10.)
1218                            (when *conv* (return nil))
1219                            (setq test nil)
1220                            (do ((i 0 (1+ i)))
1221                                ((> i *n*))
1222                              (setf (aref *hr-sl* i) (aref *shr-sl* i))
1223                              (setf (aref *hi-sl* i) (aref *shi-sl* i)))
1224                            (setq *sr* svsr *si* svsi)
1225                            (bf-polyev-sl)
1226                            (bf-calct-sl))
1227                           ((setq pasd t))))
1228                    ((setq pasd nil))))))
1229     (or *conv* (bf-vrshft-sl 10))
1230     nil))
1231
1232 (defun bf-cpoly-sl (degree)
1233   (let ( ;; Log of our floating point base.  (Do we need this much accuracy?)
1234         (*logbas* (fplog (intofp 2)))
1235         ;; "Largest" bfloat.  What should we use?
1236         (*infin* (intofp most-positive-flonum))
1237         ;; bfloat epsilon.  2^(-fpprec)
1238         (*are* (bf-scale-float (intofp 2) (- fpprec)))
1239         (*mre* (intofp 0))
1240         (xx (fproot bfhalf 2))
1241         (yy (intofp 0))
1242         ;; cos(94deg).  Probably don't need full bfloat precision here.
1243         (cosr (intofp -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804))
1244         ;; sin(94deg).  Probably don't need full bfloat precision here.
1245         (sinr (intofp 0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315))
1246         (*cr* (intofp 0))
1247         (*ci* (intofp 0))
1248         (*sr* (intofp 0))
1249         (*si* (intofp 0))
1250         (*tr* (intofp 0))
1251         (*ti* (intofp 0))
1252         (*zr* (intofp 0))
1253         (*zi* (intofp 0))
1254         (bnd (intofp 0))
1255         (*n* 0)
1256         (*polysc* 0)
1257         (*polysc1* 0)
1258         (*conv* nil))
1259     (setq *mre* (fptimes* (intofp 2)
1260                           (fptimes* (fpsqrt (intofp 2)) *are*))
1261           yy (fpminus xx))
1262     (do ((i degree (1- i)))
1263         ((not (and (equalp (intofp 0) (aref *pr-sl* i))
1264                    (equalp (intofp 0) (aref *pi-sl* i))))
1265          (setq *nn* i
1266                *n* (1- i))))
1267     (setq degree *nn*)
1268     (do ((i 0 (1+ i)))
1269         ((> i *nn*))
1270       (setf (aref *shr-sl* i) (bf-cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
1271     (if (> *nn* 0) (bf-scale-sl))
1272     (do ()
1273         ((> 2 *nn*)
1274          (bf-cdivid-sl (fpminus (aref *pr-sl* 1))
1275                        (fpminus (aref *pi-sl* 1))
1276                        (aref *pr-sl* 0)
1277                        (aref *pi-sl* 0))
1278          (setf (aref *pr-sl* 1) *cr*)
1279          (setf (aref *pi-sl* 1) *ci*)
1280          (setq *nn* 0))
1281       (do ((i 0 (1+ i)))
1282           ((> i *nn*))
1283         (setf (aref *shr-sl* i) (bf-cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
1284       (setq bnd (bf-cauchy-sl))
1285       (catch 'newroot
1286         (do ((cnt1 1 (1+ cnt1)))
1287             ((> cnt1 2))
1288           (bf-noshft-sl 5)
1289           (do ((cnt2 1 (1+ cnt2)))
1290               ((> cnt2 9))
1291             (setq xx (prog1
1292                          (fpdifference (fptimes* cosr xx)
1293                                        (fptimes* sinr yy))
1294                        (setq yy (fpplus (fptimes* sinr xx)
1295                                         (fptimes* cosr yy))))
1296                   *sr* (fptimes* bnd xx)
1297                   *si* (fptimes* bnd yy))
1298             (bf-fxshft-sl (* 10 cnt2))
1299             (cond (*conv* (setf (aref *pr-sl* *nn*) *zr*)
1300                           (setf (aref *pi-sl* *nn*) *zi*)
1301                           (setq *nn* *n* *n* (1- *n*))
1302                           (do ((i 0 (1+ i)))
1303                               ((> i *nn*))
1304                             (setf (aref *pr-sl* i) (aref *qpr-sl* i))
1305                             (setf (aref *pi-sl* i) (aref *qpi-sl* i)))
1306                           (throw 'newroot t))))))
1307       (or *conv* (return t)))
1308     (do ((i (1+ *nn*) (1+ i)))
1309         ((> i degree))
1310       (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*))
1311       (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) *polysc1*)))
1312     (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
1313         ((> i *nn*))
1314       (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) j))
1315       (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) j)))
1316     *nn*))
1317
1318
1319 (defmfun $bfallroots (expr)
1320   (prog (degree *nn* var res $partswitch
1321          ($keepfloat t)
1322          $demoivre
1323          ($listconstvars t)
1324          ($algebraic t) complex $ratfac den expr1)
1325      (setq expr1 (setq expr (meqhk expr)))
1326      (setq var (delete '$%i (cdr ($listofvars expr)) :test #'eq))
1327      (or var (setq var (list (gensym))))
1328      (cond ((not (= (length var) 1))
1329             (merror (intl:gettext "bfallroots: expected a polynomial in one variable; found variables ~M") `((mlist) ,@var)))
1330            ((setq var (car var))))
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))