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