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 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)) 2021