/maxima-5.27.0/src/cpoly.lisp
Lisp | 2020 lines | 1884 code | 62 blank | 74 comment | 0 complexity | 9664558b65e2882d388dcc44600e75cf MD5 | raw file
Possible License(s): GPL-2.0
- ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- ;;; The data in this file contains enhancments. ;;;;;
- ;;; ;;;;;
- ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
- ;;; All rights reserved ;;;;;
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- ;;; (c) Copyright 1981 Massachusetts Institute of Technology ;;;
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- (in-package :maxima)
- (macsyma-module cpoly)
- ;;; This is a lisp version of algorithm 419 from the Communications of
- ;;; the ACM (p 97 vol 15 feb 1972) by Jenkins and Traub. That
- ;;; algorithm is followed very closely. Note the following
- ;;; modifications: arrays are indexed from 0 instead of 1. This means
- ;;; that the variables n and nn are one less than the acm verson. The
- ;;; zeros are put into the arrays pr-sl and pi-sl, rather than into
- ;;; their own arrays. The algorithm seems to benefit be taking are
- ;;; mre 0.01 times the published values.
- (declare-top (special $partswitch $keepfloat $demoivre $listconstvars
- $algebraic $ratfac $programmode))
- (declare-top (special *logbas* *infin* *are* *mre* *cr* *ci* *sr* *si*
- *tr* *ti* *zr* *zi* *n* *nn* *bool*
- *conv* *pvr* *pvi* *polysc* *polysc1*))
- (declare-top (special *pr-sl* *pi-sl* *shr-sl* *shi-sl* *qpr-sl* *qpi-sl* *hr-sl*
- *hi-sl* *qhr-sl* *qhi-sl*))
- (declare-top (special *u* *v* *a* *b* *c* *d* *a1* *a3* *a7* *e* *f* *g* *h*
- *szr* *szi* *lzr* *lzi* *nz* *ui* *vi* *s*))
- (defmvar $polyfactor nil
- "When T factor the polynomial over the real or complex numbers.")
- (defmfun $allroots (expr)
- (prog (degree *nn* var res $partswitch $keepfloat $demoivre $listconstvars
- $algebraic complex $ratfac den expr1)
- (setq $keepfloat t $listconstvars t $algebraic t)
- (setq expr1 (setq expr (meqhk expr)))
- (setq var (delete '$%i (cdr ($listofvars expr)) :test #'eq))
- (or var (setq var (list (gensym))))
- (cond ((not (= (length var) 1))
- (merror (intl:gettext "allroots: expected a polynomial in one variable; found variables ~M") `((mlist) ,@var)))
- ((setq var (car var))))
- (setq expr ($rat expr '$%i var)
- res (reverse (car (cdddar expr))))
- (do ((i (- (length res) (length (caddar expr))) (1- i)))
- ((= i 0))
- (setq res (cdr res)))
- (setq den (cddr expr) expr (cadr expr))
- ;;;check denominator is a complex number
- (cond ((numberp den) (setq den (list den 0)))
- ((eq (car den) (cadr res))
- (setq den (cddr den))
- (cond ((numberp (car den))
- (cond ((null (cddr den)) (setq den (list 0 (car den))))
- ((numberp (caddr den))
- (setq den (list (caddr den) (car den))))
- (t (cpoly-err expr1))))
- (t (cpoly-err expr1))))
- (t (cpoly-err expr1)))
- ;;;if the name variable has disappeared, this is caught here
- (setq *nn* 0)
- (cond ((numberp expr) (setq expr (list expr 0)))
- ((eq (car expr) (car res)) (setq *nn* 1))
- ((eq (car expr) (cadr res))
- (setq expr (cddr expr))
- (cond ((numberp (car expr))
- (cond ((null (cddr expr)) (setq expr (list 0 (car expr))))
- ((numberp (caddr expr))
- (setq expr (list (caddr expr) (car expr))))
- (t (cpoly-err expr1))))
- (t (cpoly-err expr1))))
- (t (cpoly-err expr1)))
- (cond ((= *nn* 0)
- (cond ($polyfactor
- (let ((*cr* 0.0) (*ci* 0.0))
- (cdivid-sl (float (car expr)) (float (cadr expr))
- (float (car den)) (float (cadr den)))
- (return (simplify (list '(mplus)
- (simplify (list '(mtimes) '$%i *ci*))
- *cr*)))))
- (t (return (list '(mlist simp)))))))
- (setq degree (cadr expr) *nn* (1+ degree))
- (setq *pr-sl* (make-array *nn* :initial-element 0.0))
- (setq *pi-sl* (make-array *nn* :initial-element 0.0))
- (or (catch 'notpoly
- (errset (do ((expr (cdr expr) (cddr expr)) (l) (%i (cadr res)))
- ((null expr))
- (setq l (- degree (car expr)) res (cadr expr))
- (cond ((numberp res) (setf (aref *pr-sl* l) (float res)))
- (t
- (or (eq (car res) %i)
- (throw 'notpoly nil))
- (setq res (cddr res))
- (setf (aref *pi-sl* l) (float (car res)))
- (setq res (caddr res))
- (and res (setf (aref *pr-sl* l) (float res)))
- (setq complex t))))))
- ;;this should catch expressions like sin(x)-x
- (cpoly-err expr1))
- (setq *shr-sl* (make-array *nn* :initial-element 0.0))
- (setq *shi-sl* (make-array *nn* :initial-element 0.0))
- (setq *qpr-sl* (make-array *nn* :initial-element 0.0))
- (setq *hr-sl* (make-array degree :initial-element 0.0))
- (setq *qhr-sl* (make-array degree :initial-element 0.0))
- (setq *qpi-sl* (make-array *nn* :initial-element 0.0))
- (when complex
- (setq *hi-sl* (make-array degree :initial-element 0.0))
- (setq *qhi-sl* (make-array degree :initial-element 0.0)))
- (setq *nn* degree)
- (if complex
- (setq res (errset (cpoly-sl degree)))
- (setq res (errset (rpoly-sl degree))))
- (unless res
- (mtell (intl:gettext "allroots: unexpected error; treat results with caution.~%")))
- (when (= *nn* degree)
- (merror (intl:gettext "allroots: no roots found.")))
- (setq res nil)
- (cond ((not (zerop *nn*))
- (mtell (intl:gettext "allroots: only ~S out of ~S roots found.~%") (- degree *nn*) degree)
- (setq expr 0.0)
- (do ((i 0 (1+ i)))
- ((> i *nn*))
- (setq expr
- (simplify
- (list '(mplus) expr
- (simplify (list '(mtimes)
- (simplify (list '(mplus)
- (simplify (list '(mtimes) '$%i (aref *pi-sl* i)))
- (aref *pr-sl* i)))
- (simplify (list '(mexpt) var (- *nn* i)))))))))
- (setq res (cons expr res)))
- ($polyfactor
- (setq expr (let ((*cr* 0.0) (*ci* 0.0))
- (cdivid-sl (aref *pr-sl* 0) (aref *pi-sl* 0)
- (float (car den))
- (float (cadr den)))
- (simplify (list '(mplus) (simplify (list '(mtimes) '$%i *ci*)) *cr*)))
- res (cons expr res))))
- (do ((i degree (1- i)))
- ((= i *nn*))
- (setq expr (simplify (list '(mplus)
- (simplify (list '(mtimes) '$%i (aref *pi-sl* i)))
- (aref *pr-sl* i))))
- (setq res
- (cond ($polyfactor (cons (cond ((or complex (zerop (aref *pi-sl* i)))
- (simplify (list '(mplus) var (simplify (list '(mminus) expr)))))
- (t (setq i (1- i))
- (simplify (list '(mplus)
- (simplify (list '(mexpt) var 2))
- (simplify (list '(mtimes) var
- (aref *pr-sl* i)))
- (aref *pr-sl* (1+ i))))))
- res))
- ((cons (let ((expr (simplify (list '(mequal) var expr))))
- (if $programmode expr (displine expr)))
- res)))))
- (return (simplify (if $polyfactor
- (cons '(mtimes) res)
- (cons '(mlist) (nreverse res)))))))
- (defun cpoly-err (expr)
- (merror (intl:gettext "allroots: expected a polynomial; found ~M") expr))
- (defun cpoly-sl (degree)
- (let ((*logbas* (log 2.0))
- (*infin* most-positive-flonum)
- (*are* flonum-epsilon)
- (*mre* 0.0)
- (xx (sqrt 0.5))
- (yy 0.0)
- (cosr (cos (float (* 94/180 pi))))
- (sinr (sin (float (* 94/180 pi))))
- (*cr* 0.0) (*ci* 0.0)
- (*sr* 0.0) (*si* 0.0)
- (*tr* 0.0) (*ti* 0.0)
- (*zr* 0.0) (*zi* 0.0)
- (bnd 0.0)
- (*n* 0)
- (*polysc* 0)
- (*polysc1* 0)
- (*conv* nil))
- (setq *mre* (* 2.0 (sqrt 2.0) *are*)
- yy (- xx))
- (do ((i degree (1- i)))
- ((not (and (zerop (aref *pr-sl* i)) (zerop (aref *pi-sl* i))))
- (setq *nn* i
- *n* (1- i))))
- (setq degree *nn*)
- (do ((i 0 (1+ i)))
- ((> i *nn*))
- (setf (aref *shr-sl* i) (cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
- (if (> *nn* 0) (scale-sl))
- (do ()
- ((> 2 *nn*)
- (cdivid-sl (- (aref *pr-sl* 1)) (- (aref *pi-sl* 1))
- (aref *pr-sl* 0) (aref *pi-sl* 0))
- (setf (aref *pr-sl* 1) *cr*)
- (setf (aref *pi-sl* 1) *ci*)
- (setq *nn* 0))
- (do ((i 0 (1+ i)))
- ((> i *nn*))
- (setf (aref *shr-sl* i) (cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
- (setq bnd (cauchy-sl))
- (catch 'newroot
- (do ((cnt1 1 (1+ cnt1)))
- ((> cnt1 2))
- (noshft-sl 5)
- (do ((cnt2 1 (1+ cnt2)))
- ((> cnt2 9))
- (setq xx (prog1
- (- (* cosr xx) (* sinr yy))
- (setq yy (+ (* sinr xx) (* cosr yy))))
- *sr* (* bnd xx)
- *si* (* bnd yy))
- (fxshft-sl (* 10 cnt2))
- (cond (*conv* (setf (aref *pr-sl* *nn*) *zr*)
- (setf (aref *pi-sl* *nn*) *zi*)
- (setq *nn* *n* *n* (1- *n*))
- (do ((i 0 (1+ i)))
- ((> i *nn*))
- (setf (aref *pr-sl* i) (aref *qpr-sl* i))
- (setf (aref *pi-sl* i) (aref *qpi-sl* i)))
- (throw 'newroot t))))))
- (or *conv* (return t)))
- (do ((i (1+ *nn*) (1+ i)))
- ((> i degree))
- (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*))
- (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) *polysc1*)))
- (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
- ((> i *nn*))
- (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))
- (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) j)))
- *nn*))
- (defun noshft-sl (l1)
- (do ((i 0 (1+ i))
- (xni (float *nn*) (1- xni))
- (t1 (/ (float *nn*))))
- ((> i *n*))
- (setf (aref *hr-sl* i) (* (aref *pr-sl* i) xni t1))
- (setf (aref *hi-sl* i) (* (aref *pi-sl* i) xni t1)))
- (do ((jj 1 (1+ jj)))
- ((> jj l1))
- (cond ((> (cmod-sl (aref *hr-sl* *n*) (aref *hi-sl* *n*))
- (* 10.0 *are* (cmod-sl (aref *pr-sl* *n*) (aref *pi-sl* *n*))))
- (cdivid-sl (- (aref *pr-sl* *nn*)) (- (aref *pi-sl* *nn*)) (aref *hr-sl* *n*) (aref *hi-sl* *n*))
- (setq *tr* *cr* *ti* *ci*)
- (do ((j *n* (1- j)) (t1) (t2))
- ((> 1 j))
- (setq t1 (aref *hr-sl* (1- j)) t2 (aref *hi-sl* (1- j)))
- (setf (aref *hr-sl* j) (- (+ (aref *pr-sl* j) (* t1 *tr*)) (* t2 *ti*)))
- (setf (aref *hi-sl* j) (+ (aref *pi-sl* j) (* t1 *ti*) (* t2 *tr*))))
- (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
- (setf (aref *hi-sl* 0) (aref *pi-sl* 0)))
- (t (do ((j *n* (1- j)))
- ((> 1 j))
- (setf (aref *hr-sl* j) (aref *hr-sl* (1- j)))
- (setf (aref *hi-sl* j) (aref *hi-sl* (1- j))))
- (setf (aref *hr-sl* 0) 0.0)
- (setf (aref *hi-sl* 0) 0.0)))))
- (defun fxshft-sl (l2)
- (let ((test t)
- (pasd nil)
- (otr 0.0) (oti 0.0)
- (svsr 0.0) (svsi 0.0)
- (*bool* nil)
- (*pvr* 0.0) (*pvi* 0.0))
- (polyev-sl)
- (setq *conv* nil)
- (calct-sl)
- (do ((j 1 (1+ j)))
- ((> j l2))
- (setq otr *tr* oti *ti*)
- (nexth-sl)
- (calct-sl)
- (setq *zr* (+ *sr* *tr*) *zi* (+ *si* *ti*))
- (cond ((and (not *bool*) test (not (= j l2)))
- (cond ((> (* 0.5 (cmod-sl *zr* *zi*))
- (cmod-sl (- *tr* otr) (- *ti* oti)))
- (cond (pasd (do ((i 0 (1+ i)))
- ((> i *n*))
- (setf (aref *shr-sl* i) (aref *hr-sl* i))
- (setf (aref *shi-sl* i) (aref *hi-sl* i)))
- (setq svsr *sr* svsi *si*)
- (vrshft-sl 10.)
- (when *conv* (return nil))
- (setq test nil)
- (do ((i 0 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i) (aref *shr-sl* i))
- (setf (aref *hi-sl* i) (aref *shi-sl* i)))
- (setq *sr* svsr *si* svsi)
- (polyev-sl)
- (calct-sl))
- ((setq pasd t))))
- ((setq pasd nil))))))
- (or *conv* (vrshft-sl 10))
- nil))
- (defun vrshft-sl (l3)
- (setq *conv* nil *sr* *zr* *si* *zi*)
- (do ((i 1 (1+ i))
- (bool1 nil)
- (mp) (ms) (omp) (relstp) (tp) (r1))
- ((> i l3))
- (polyev-sl)
- (setq mp (cmod-sl *pvr* *pvi*) ms (cmod-sl *sr* *si*))
- (cond ((> (* 20.0 (errev-sl ms mp)) mp)
- (setq *conv* t *zr* *sr* *zi* *si*)
- (return t)))
- (cond ((= i 1) (setq omp mp))
- ((or bool1 (> omp mp) (not (< relstp 0.05)))
- (if (> (* 0.1 mp) omp)
- (return t)
- (setq omp mp)))
- (t (setq tp relstp bool1 t)
- (when (> *are* relstp) (setq tp *are*))
- (setq r1 (sqrt tp)
- *sr* (prog1
- (- (* (1+ r1) *sr*) (* r1 *si*))
- (setq *si* (+ (* (1+ r1) *si*) (* r1 *sr*)))))
- (polyev-sl)
- (do ((j 1 (1+ j))) ((> j 5)) (calct-sl) (nexth-sl))
- (setq omp *infin*)))
- (calct-sl)
- (nexth-sl)
- (calct-sl)
- (or *bool*
- (setq relstp (/ (cmod-sl *tr* *ti*) (cmod-sl *sr* *si*))
- *sr* (+ *sr* *tr*)
- *si* (+ *si* *ti*)))))
- (defun calct-sl nil
- (do ((i 1 (1+ i))
- ($t)
- (hvr (setf (aref *qhr-sl* 0) (aref *hr-sl* 0)))
- (hvi (setf (aref *qhi-sl* 0) (aref *hi-sl* 0))))
- ((> i *n*)
- (setq *bool* (not (> (cmod-sl hvr hvi) (* 10.0 *are* (cmod-sl (aref *hr-sl* *n*) (aref *hi-sl* *n*))))))
- (cond ((not *bool*) (cdivid-sl (- *pvr*) (- *pvi*) hvr hvi) (setq *tr* *cr* *ti* *ci*))
- (t (setq *tr* 0.0 *ti* 0.0)))
- nil)
- (setq $t (- (+ (aref *hr-sl* i) (* hvr *sr*)) (* hvi *si*)))
- (setf (aref *qhi-sl* i) (setq hvi (+ (aref *hi-sl* i) (* hvr *si*) (* hvi *sr*))))
- (setf (aref *qhr-sl* i) (setq hvr $t))))
- (defun nexth-sl ()
- (cond (*bool* (do ((j 1 (1+ j)))
- ((> j *n*))
- (setf (aref *hr-sl* j) (aref *qhr-sl* (1- j)))
- (setf (aref *hi-sl* j) (aref *qhi-sl* (1- j))))
- (setf (aref *hr-sl* 0) 0.0)
- (setf (aref *hi-sl* 0) 0.0))
- (t (do ((j 1. (1+ j)) (t1) (t2))
- ((> j *n*))
- (setq t1 (aref *qhr-sl* (1- j)) t2 (aref *qhi-sl* (1- j)))
- (setf (aref *hr-sl* j) (- (+ (aref *qpr-sl* j) (* t1 *tr*)) (* t2 *ti*)))
- (setf (aref *hi-sl* j) (+ (aref *qpi-sl* j) (* t1 *ti*) (* t2 *tr*))))
- (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
- (setf (aref *hi-sl* 0) (aref *qpi-sl* 0))))
- nil)
- (defun polyev-sl ()
- (setq *pvr* (setf (aref *qpr-sl* 0) (aref *pr-sl* 0))
- *pvi* (setf (aref *qpi-sl* 0) (aref *pi-sl* 0)))
- (do ((i 1 (1+ i)) ($t))
- ((> i *nn*))
- (setq $t (- (+ (aref *pr-sl* i) (* *pvr* *sr*)) (* *pvi* *si*)))
- (setf (aref *qpi-sl* i) (setq *pvi* (+ (aref *pi-sl* i) (* *pvr* *si*) (* *pvi* *sr*))))
- (setf (aref *qpr-sl* i) (setq *pvr* $t))))
- (defun errev-sl (ms mp)
- (- (* (do ((j 0 (1+ j))
- (*e* (/ (* (cmod-sl (aref *qpr-sl* 0) (aref *qpi-sl* 0)) *mre*) (+ *are* *mre*))))
- ((> j *nn*) *e*)
- (setq *e* (+ (cmod-sl (aref *qpr-sl* j) (aref *qpi-sl* j)) (* *e* ms))))
- (+ *are* *mre*))
- (* mp *mre*)))
- ;; Compute a lower bound on the roots of the polynomial. Let the
- ;; polynomial be sum(a[k]*x^k,k,0,n). Then the lower bound is the
- ;; smallest real root of sum(|a[k]|*x^k,k,0,n-1)-a[n]. For our
- ;; purposes, this lower bound is computed to an accuracy of .005 or
- ;; so.
- (defun cauchy-sl ()
- (let ((x (exp (/ (- (log (aref *shr-sl* *nn*)) (log (aref *shr-sl* 0))) (float *nn*))))
- (xm 0.0))
- (setf (aref *shr-sl* *nn*) (- (aref *shr-sl* *nn*)))
- (cond ((not (zerop (aref *shr-sl* *n*)))
- (setq xm (- (/ (aref *shr-sl* *nn*) (aref *shr-sl* *n*))))
- (cond ((> x xm) (setq x xm)))))
- (do ((*f*))
- (nil)
- (setq xm (* 0.1 x) *f* (aref *shr-sl* 0))
- (do ((i 1 (1+ i))) ((> i *nn*)) (setq *f* (+ (aref *shr-sl* i) (* *f* xm))))
- (cond ((not (< 0.0 *f*)) (return t)))
- (setq x xm))
- (do ((dx x) (df) (*f*))
- ((> 5e-3 (abs (/ dx x))) x)
- (setq *f* (aref *shr-sl* 0) df *f*)
- (do ((i 1 (1+ i)))
- ((> i *n*))
- (setq *f* (+ (* *f* x) (aref *shr-sl* i)) df (+ (* df x) *f*)))
- (setq *f* (+ (* *f* x) (aref *shr-sl* *nn*)) dx (/ *f* df) x (- x dx)))))
- ;; Scale the coefficients of the polynomial to reduce the chance of
- ;; overflow or underflow. This is different from the algorithm given
- ;; in TOMS 419 and 493 which just scales the coefficients by some
- ;; fixed factor. Instead, this algorithm computes a scale factor for
- ;; the roots and adjusts the coefficients of the polynomial
- ;; appropriately.
- ;;
- ;; The scale factor for the roots is in *polysc1*. The roots are
- ;; scaled by 2^(*polysc1*). There is an additional scaling *polysc*
- ;; such that the coefficients of the polynomial are all scaled by
- ;; 2^(*polysc*).
- (defun scale-sl ()
- (do ((i 0 (1+ i)) (j 0) (x 0.0) (dx 0.0))
- ((> i *nn*)
- (setq x (/ x (float (- (1+ *nn*) j)))
- dx (/ (- (log (aref *shr-sl* *nn*)) (log (aref *shr-sl* 0))) (float *nn*))
- *polysc1* (floor (+ 0.5 (/ dx *logbas*)))
- x (+ x (* (float (* *polysc1* *nn*)) *logbas* 0.5))
- *polysc* (floor (+ 0.5 (/ x *logbas*)))))
- (cond ((zerop (aref *shr-sl* i)) (setq j (1+ j)))
- (t (setq x (+ x (log (aref *shr-sl* i)))))))
- (do ((i *nn* (1- i)) (j (- *polysc*) (+ j *polysc1*)))
- ((< i 0))
- (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j))
- (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) j))))
- (defun cdivid-sl (ar ai br bi)
- (let ((r1 0.0))
- (cond ((and (zerop br) (zerop bi))
- (setq *cr* (setq *ci* *infin*)))
- ((> (abs bi) (abs br))
- (setq r1 (/ br bi)
- bi (+ bi (* br r1))
- br (+ ai (* ar r1))
- *cr* (/ br bi)
- br (- (* ai r1) ar)
- *ci* (/ br bi)))
- ((setq r1 (/ bi br)
- bi (+ br (* bi r1))
- br (+ ar (* ai r1))
- *cr* (/ br bi)
- br (- ai (* ar r1))
- *ci* (/ br bi)))))
- nil)
- (defun cmod-sl (ar ai)
- (setq ar (abs ar)
- ai (abs ai))
- (cond ((> ai ar) (setq ar (/ ar ai)) (* ai (sqrt (1+ (* ar ar)))))
- ((> ar ai) (setq ai (/ ai ar)) (* ar (sqrt (1+ (* ai ai)))))
- ((* (sqrt 2.0)
- ar))))
- ;;; This is the algorithm for doing real polynomials. It is algorithm
- ;;; 493 from acm toms vol 1 p 178 (1975) by jenkins. Note that array
- ;;; indexing starts from 0. The names of the arrays have been changed
- ;;; to be the same as for cpoly. The correspondence is: p - pr-sl, qp
- ;;; - qpr-sl, k - hr-sl, qk - qhr-sl, svk - shr-sl, temp - shi-sl.
- ;;; the roots *are* put in pr-sl and pi-sl. The variable *si* appears
- ;;; not to be used here
- (defun rpoly-sl (degree)
- (let ((*logbas* (log 2.0))
- (*infin* most-positive-flonum)
- (*are* flonum-epsilon)
- (*mre* 0.0)
- (xx (sqrt 0.5)) ;; sqrt(0.5)
- (yy 0.0)
- (cosr (cos (float (* 94/180 pi))))
- (sinr (sin (float (* 94/180 pi))))
- (aa 0.0)
- (cc 0.0)
- (bb 0.0)
- (bnd 0.0)
- (*sr* 0.0)
- (*u* 0.0)
- (*v* 0.0)
- (t1 0.0)
- (*szr* 0.0) (*szi* 0.0)
- (*lzr* 0.0) (*lzi* 0.0)
- (*nz* 0)
- (*n* 0)
- (*polysc* 0)
- (*polysc1* 0)
- (zerok 0)
- (conv1 t))
- (setq *mre* *are* yy (- xx))
- (do ((i degree (1- i))) ((not (zerop (aref *pr-sl* i))) (setq *nn* i *n* (1- i))))
- (setq degree *nn*)
- (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i))))
- (if (> *nn* 0) (scale-sl))
- ;; Loop to find all roots
- (do nil
- ((< *nn* 3)
- (cond ((= *nn* 2)
- ;; Solve the final quadratic polynomial
- (quad-sl (aref *pr-sl* 0.) (aref *pr-sl* 1) (aref *pr-sl* 2))
- (cond ((and $polyfactor (not (zerop *szi*)))
- (setf (aref *pr-sl* 2) (/ (aref *pr-sl* 2) (aref *pr-sl* 0)))
- (setf (aref *pr-sl* 1) (/ (aref *pr-sl* 1) (aref *pr-sl* 0)))
- (setf (aref *pi-sl* 2) 1.0))
- (t (setf (aref *pr-sl* 2) *szr*)
- (setf (aref *pi-sl* 2) *szi*)
- (setf (aref *pr-sl* 1) *lzr*)
- (setf (aref *pi-sl* 1) *lzi*))))
- (t
- ;; Solve the final linear polynomial
- (setf (aref *pr-sl* 1) (- (/ (aref *pr-sl* 1) (aref *pr-sl* 0))))))
- (setq *nn* 0))
- ;; Calculate a lower bound on the modulus of the zeros.
- (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i))))
- (setq bnd (cauchy-sl))
- ;; Compute derivative of polynomial. Save result in *hr-sl*.
- (do ((i 1 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i) (/ (* (float (- *n* i)) (aref *pr-sl* i)) (float *n*))))
- (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
- (setq aa (aref *pr-sl* *nn*) bb (aref *pr-sl* *n*) zerok (zerop (aref *hr-sl* *n*)))
- ;; Do 5 steps with no shift.
- (do ((jj 1 (1+ jj)))
- ((> jj 5.))
- (setq cc (aref *hr-sl* *n*))
- (cond (zerok (do ((j *n* (1- j)))
- ((< j 1))
- (setf (aref *hr-sl* j) (aref *hr-sl* (1- j))))
- (setf (aref *hr-sl* 0) 0.0)
- (setq zerok (zerop (aref *hr-sl* *n*))))
- (t (setq t1 (- (/ aa cc)))
- (do ((j *n* (1- j)))
- ((< j 1))
- (setf (aref *hr-sl* j) (+ (* t1 (aref *hr-sl* (1- j))) (aref *pr-sl* j))))
- (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
- (setq zerok (not (> (abs (aref *hr-sl* *n*))
- (* (abs bb) *are* 10.0)))))))
- (do ((i 0 (1+ i))) ((> i *n*)) (setf (aref *shi-sl* i) (aref *hr-sl* i)))
- (do ((cnt 1 (1+ cnt)))
- ((> cnt 20.) (setq conv1 nil))
- (setq xx (prog1
- (- (* cosr xx) (* sinr yy))
- (setq yy (+ (* sinr xx) (* cosr yy))))
- *sr* (* bnd xx)
- *u* (* -2.0 *sr*)
- *v* bnd)
- (fxshfr-sl (* 20 cnt))
- (cond ((> *nz* 0)
- (setf (aref *pr-sl* *nn*) *szr*)
- (setf (aref *pi-sl* *nn*) *szi*)
- (cond ((= *nz* 2)
- (setf (aref *pr-sl* *n*) *lzr*)
- (setf (aref *pi-sl* *n*) *lzi*)
- (cond ((and $polyfactor (not (zerop *szi*)))
- (setf (aref *pr-sl* *nn*) *v*)
- (setf (aref *pr-sl* *n*) *u*)
- (setf (aref *pi-sl* *nn*) 1.0)))))
- (setq *nn* (- *nn* *nz*) *n* (1- *nn*))
- (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *pr-sl* i) (aref *qpr-sl* i)))
- (return nil)))
- (do ((i 0 (1+ i))) ((> i *n*)) (setf (aref *hr-sl* i) (aref *shi-sl* i))))
- (or conv1 (return nil)))
- (cond ($polyfactor
- (do ((i degree (1- i)))
- ((= i *nn*))
- (cond ((zerop (aref *pi-sl* i))
- (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*)))
- (t (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) (* 2 *polysc1*)))
- (setq i (1- i))
- (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*))))))
- (t (do ((i (1+ *nn*) (1+ i)))
- ((> i degree))
- (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) *polysc1*))
- (setf (aref *pi-sl* i) (scale-float (aref *pi-sl* i) *polysc1*)))))
- (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
- ((> i *nn*))
- (setf (aref *pr-sl* i) (scale-float (aref *pr-sl* i) j)))))
- (defun fxshfr-sl (l2)
- (let ((*my-type* 0)
- (*a* 0.0) (*b* 0.0) (*c* 0.0) (*d* 0.0) (*e* 0.0) (*f* 0.0) (*g* 0.0) (*h* 0.0)
- (*a1* 0.0) (*a3* 0.0) (*a7* 0.0))
- (declare (special *my-type*))
- (setq *nz* 0)
- (quadsd-sl)
- (calcsc-sl)
- (do ((j 1 (1+ j))
- (betav 0.25)
- (betas 0.25)
- (oss *sr*)
- (ovv *v*)
- (tvv) (tss) (ss) (vv) (tv) (ts) (ots) (otv)
- (*ui*) (*vi*) (*s*) (svv) (svu) (iflag) (vpass) (spass) (vtry) (stry))
- ((> j l2))
- (nextk-sl)
- (calcsc-sl)
- (newest-sl)
- (setq vv *vi*
- ss 0.0)
- (or (zerop (aref *hr-sl* *n*))
- (setq ss (- (/ (aref *pr-sl* *nn*) (aref *hr-sl* *n*)))))
- (setq tv 1.0 ts 1.0)
- (cond ((not (or (= j 1) (= *my-type* 3)))
- (or (zerop vv) (setq tv (abs (/ (- vv ovv) vv))))
- (or (zerop ss) (setq ts (abs (/ (- ss oss) ss))))
- (setq tvv 1.0)
- (and (< tv otv) (setq tvv (* tv otv)))
- (setq tss 1.0)
- (and (< ts ots) (setq tss (* ts ots)))
- (setq vpass (< tvv betav) spass (< tss betas))
- (cond ((or spass vpass)
- (setq svu *u* svv *v*)
- (do ((i 0 (1+ i)))
- ((> i *n*)) (setf (aref *shr-sl* i)
- (aref *hr-sl* i)))
- (setq *s* ss vtry nil stry nil)
- (and (do ((*bool* (not (and spass (or (not vpass) (< tss tvv)))) t)
- (l50 nil nil))
- (nil)
- (cond (*bool* (quadit-sl)
- (and (> *nz* 0) (return t))
- (setq vtry t
- betav (* 0.25 betav))
- (cond ((or stry (not spass))
- (setq l50 t))
- (t (do ((i 0 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (aref *shr-sl* i)))))))
- (cond ((not l50)
- (setq iflag (realit-sl))
- (and (> *nz* 0) (return t))
- (setq stry t betas (* 0.25 betas))
- (cond ((zerop iflag) (setq l50 t))
- (t (setq *ui* (- (+ *s* *s*))
- *vi* (* *s* *s*))))))
- (cond (l50 (setq *u* svu *v* svv)
- (do ((i 0 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (aref *shr-sl* i)))
- (and (or (not vpass) vtry)
- (return nil)))))
- (return nil))
- (quadsd-sl)
- (calcsc-sl)))))
- (setq ovv vv
- oss ss
- otv tv
- ots ts))))
- (defun quadit-sl nil
- (setq *nz* 0 *u* *ui* *v* *vi*)
- (do ((tried) (j 0) (ee) (mp) (relstp) (omp) (ms))
- (nil)
- (quad-sl 1.0 *u* *v*)
- (and (> (abs (- (abs *szr*) (abs *lzr*))) (* 1e-2 (abs *lzr*)))
- (return nil))
- (quadsd-sl)
- (setf mp (+ (abs (- *a* (* *szr* *b*))) (abs (* *szi* *b*))))
- (setf ms (cmod-sl *szr* *szi*))
- (setf ee (errev-sl ms mp))
- (cond ((not (> mp (* 2e1 ee))) (setq *nz* 2)
- (return nil)))
- (setq j (1+ j))
- (and (> j 20) (return nil))
- (cond ((not (or (< j 2) (> relstp 1e-2) (< mp omp) tried))
- (and (< relstp *are*) (setq relstp *are*))
- (setq relstp (sqrt relstp)
- *u* (- *u* (* *u* relstp))
- *v* (+ *v* (* *v* relstp)))
- (quadsd-sl)
- (do ((i 1 (1+ i)))
- ((> i 5)) (calcsc-sl) (nextk-sl))
- (setq tried t j 0)))
- (setq omp mp)
- (calcsc-sl)
- (nextk-sl)
- (calcsc-sl)
- (newest-sl)
- (and (zerop *vi*) (return nil))
- (setq relstp (abs (/ (- *vi* *v*) *vi*)) *u* *ui* *v* *vi*)))
- (defun realit-sl ()
- (setq *nz* 0)
- (do ((j 0) (pv) (ee) (ms) (mp) (kv) (t1) (omp))
- (nil)
- (setq pv (aref *pr-sl* 0))
- (setf (aref *qpr-sl* 0) pv)
- (do ((i 1 (1+ i)))
- ((> i *nn*))
- (setq pv (+ (* pv *s*) (aref *pr-sl* i)))
- (setf (aref *qpr-sl* i) pv))
- (setq mp (abs pv)
- ms (abs *s*)
- ee (* (/ *mre* (+ *are* *mre*)) (abs (aref *qpr-sl* 0))))
- (do ((i 1 (1+ i)))
- ((> i *nn*)) (setq ee (+ (* ee ms) (abs (aref *qpr-sl* i)))))
- (cond ((not (> mp (* 2e1 (- (* (+ *are* *mre*) ee) (* *mre* mp)))))
- (setq *nz* 1 *szr* *s* *szi* 0.0)
- (return 0)))
- (setq j (1+ j))
- (and (> j 10) (return 0))
- (cond ((not (or (< j 2)
- (> (abs t1) (* 1e-3 (abs (- *s* t1))))
- (not (> mp omp))))
- (return 1)))
- (setq omp mp kv (aref *hr-sl* 0))
- (setf (aref *qhr-sl* 0) kv)
- (do ((i 1 (1+ i)))
- ((> i *n*))
- (setq kv (+ (* kv *s*) (aref *hr-sl* i)))
- (setf (aref *qhr-sl* i) kv))
- (cond ((> (abs kv) (* (abs (aref *hr-sl* *n*)) 1e1 *are*))
- (setq t1 (- (/ pv kv)))
- (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
- (do ((i 1 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (+ (* t1 (aref *qhr-sl* (1- i))) (aref *qpr-sl* i)))))
- (t (setf (aref *hr-sl* 0) 0.0)
- (do ((i 1 (1+ i)))
- ((> i *n*)) (setf (aref *hr-sl* i) (aref *qhr-sl* (1- i))))))
- (setq kv (aref *hr-sl* 0))
- (do ((i 1 (1+ i)))
- ((> i *n*)) (setq kv (+ (* kv *s*) (aref *hr-sl* i))))
- (setq t1 0.0)
- (and (> (abs kv) (* (abs (aref *hr-sl* *n*)) 10.0 *are*))
- (setq t1 (- (/ pv kv))))
- (setq *s* (+ *s* t1))))
- (defun calcsc-sl ()
- (declare (special *my-type*))
- (setq *d* (aref *hr-sl* 0))
- (setf (aref *qhr-sl* 0) *d*)
- (setq *c* (- (aref *hr-sl* 1) (* *u* *d*)))
- (setf (aref *qhr-sl* 1) *c*)
- (do ((i 2 (1+ i))
- (c0))
- ((> i *n*))
- (setq c0 (- (aref *hr-sl* i) (* *u* *c*) (* *v* *d*)))
- (setf (aref *qhr-sl* i) c0)
- (setq *d* *c* *c* c0))
- (cond ((not (or (> (abs *c*) (* (abs (aref *hr-sl* *n*)) 1e2 *are*))
- (> (abs *d*) (* (abs (aref *hr-sl* (1- *n*))) 1e2 *are*))))
- (setq *my-type* 3))
- ((not (< (abs *d*) (abs *c*)))
- (setq *my-type* 2
- *e* (/ *a* *d*)
- *f* (/ *c* *d*)
- *g* (* *u* *b*)
- *h* (* *v* *b*)
- *a3* (+ (* (+ *a* *g*) *e*) (* *h* (/ *b* *d*)))
- *a1* (- (* *b* *f*) *a*)
- *a7* (+ (* (+ *f* *u*) *a*) *h*)))
- (t (setq *my-type* 1
- *e* (/ *a* *c*)
- *f* (/ *d* *c*)
- *g* (* *u* *e*)
- *h* (* *v* *b*)
- *a3* (+ (* *a* *e*) (* (+ (/ *h* *c*) *g*) *b*))
- *a1* (- *b* (* *a* (/ *d* *c*)))
- *a7* (+ *a* (* *g* *d*) (* *h* *f*)))))
- nil)
- (defun nextk-sl ()
- (declare (special *my-type*))
- (cond ((= *my-type* 3)
- (setf (aref *hr-sl* 0) 0.0)
- (setf (aref *hr-sl* 1) 0.0)
- (do ((i 2 (1+ i)))
- ((> i *n*)) (setf (aref *hr-sl* i) (aref *qhr-sl* (- i 2)))))
- ((> (abs *a1*) (* (abs (if (= *my-type* 1) *b* *a*)) 1e1 *are*))
- (setq *a7* (/ *a7* *a1*) *a3* (/ *a3* *a1*))
- (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
- (setf (aref *hr-sl* 1) (- (aref *qpr-sl* 1) (* *a7* (aref *qpr-sl* 0))))
- (do ((i 2 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (+ (* *a3* (aref *qhr-sl* (- i 2)))
- (- (* *a7* (aref *qpr-sl* (1- i))))
- (aref *qpr-sl* i)))))
- (t (setf (aref *hr-sl* 0) 0.0)
- (setf (aref *hr-sl* 1) (- (* *a7* (aref *qpr-sl* 0))))
- (do ((i 2 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (- (* *a3* (aref *qhr-sl* (- i 2)))
- (* *a7* (aref *qpr-sl* (1- i))))))))
- nil)
- (defun newest-sl ()
- (declare (special *my-type*))
- (let ((a4 0.0) (a5 0.0)
- (b1 0.0) (b2 0.0)
- (c1 0.0) (c2 0.0) (c3 0.0) (c4 0.0))
- (cond ((= *my-type* 3)
- (setq *ui* 0.0 *vi* 0.0))
- (t
- (if (= *my-type* 2)
- (setq a4 (+ (* (+ *a* *g*) *f*) *h*)
- a5 (+ (* (+ *f* *u*) *c*) (* *v* *d*)))
- (setq a4 (+ *a* (* *u* *b*) (* *h* *f*))
- a5 (+ *c* (* (+ *u* (* *v* *f*)) *d*))))
- (setq b1 (- (/ (aref *hr-sl* *n*) (aref *pr-sl* *nn*)))
- b2 (- (/ (+ (aref *hr-sl* (1- *n*)) (* b1 (aref *pr-sl* *n*))) (aref *pr-sl* *nn*)))
- c1 (* *v* b2 *a1*)
- c2 (* b1 *a7*)
- c3 (* b1 b1 *a3*)
- c4 (- c1 c2 c3)
- c1 (+ a5 (* b1 a4) (- c4)))
- (if (zerop c1)
- (setq *ui* 0.0 *vi* 0.0)
- (setq *ui* (- *u* (/ (+ (* *u* (+ c3 c2))
- (* *v* (+ (* b1 *a1*) (* b2 *a7*))))
- c1))
- *vi* (* *v* (1+ (/ c4 c1)))))))
- nil))
- ;; Divide the polynomial in *pr-sl* by the quadratic x^2 + (*u*)*x +
- ;; (*v*). Place the quotient in *qpr-sl* and the remainder in *a* and
- ;; *b*. I (rtoy) think the remainder polynomial is (*b*)*x + (*a*).
- (defun quadsd-sl ()
- (setq *b* (aref *pr-sl* 0))
- (setf (aref *qpr-sl* 0) *b*)
- (setq *a* (- (aref *pr-sl* 1) (* *u* *b*)))
- (setf (aref *qpr-sl* 1) *a*)
- (do ((i 2 (1+ i))
- (c0))
- ((> i *nn*))
- (setq c0 (- (aref *pr-sl* i) (* *u* *a*) (* *v* *b*)))
- (setf (aref *qpr-sl* i) c0)
- (setq *b* *a*
- *a* c0)))
- ;; Compute the zeros of the quadratic a0*z^2+b1*z+c0. The larger zero
- ;; is returned in *szr* and *szi*. The smaller zero is in *lzr* and
- ;; *lzi*.
- (defun quad-sl (a0 b1 c0)
- (setq *szr* 0.0 *szi* 0.0 *lzr* 0.0 *lzi* 0.0)
- (let ((b0 0.0)
- (l0 0.0)
- (*e* 0.0))
- ;; Handle the degenerate cases of a0 = 0 or c0 = 0 first.
- (cond ((zerop a0) (or (zerop b1) (setq *szr* (- (/ c0 b1)))))
- ((zerop c0) (setq *lzr* (- (/ b1 a0))))
- (t
- ;; Quadratic formula.
- (setq b0 (/ b1 2.0))
- (cond ((< (abs b0) (abs c0))
- (setq *e* a0)
- (and (< c0 0.0) (setq *e* (- a0)))
- (setq *e* (- (* b0 (/ b0 (abs c0))) *e*)
- l0 (* (sqrt (abs *e*)) (sqrt (abs c0)))))
- (t (setq *e* (- 1.0 (* (/ a0 b0) (/ c0 b0)))
- l0 (* (sqrt (abs *e*)) (abs b0)))))
- (cond ((< *e* 0.0)
- (setq *szr* (- (/ b0 a0))
- *lzr* *szr*
- *szi* (abs (/ l0 a0))
- *lzi* (- *szi*)))
- (t (or (< b0 0.0) (setq l0 (- l0)))
- (setq *lzr* (/ (- l0 b0) a0))
- (or (zerop *lzr*) (setq *szr* (/ c0 *lzr* a0)))))))
- nil))
- ;; This is a very straightforward conversion of $allroots to use
- ;; bfloats instead of floats.
- (defun bf-cpoly-err (expr)
- (merror (intl:gettext "bfallroots: expected a polynomial; found ~M") expr))
- (defun fpzerop (x)
- (equal '(0 0) x))
- ;; (ar+%i*ai)/(br+%i*bi) -> cr+%i*ci.
- (defun bf-cdivid-sl (ar ai br bi)
- (cond ((and (fpzerop br)
- (fpzerop bi))
- ;; Division by zero. Should we do something else besides set
- ;; both parts to be "infinity"?
- (setq *cr* (setq *ci* *infin*)))
- ((fpgreaterp (fpabs bi) (fpabs br))
- (let ((r1 (fpquotient br bi)))
- (setq bi (fpplus bi (fptimes* br r1))
- br (fpplus ai (fptimes* ar r1))
- *cr* (fpquotient br bi)
- br (fpdifference (fptimes* ai r1) ar)
- *ci* (fpquotient br bi))))
- (t
- (let ((r1 (fpquotient bi br)))
- (setq bi (fpplus br (fptimes* bi r1))
- br (fpplus ar (fptimes* ai r1))
- *cr* (fpquotient br bi)
- br (fpdifference ai (fptimes* ar r1))
- *ci* (fpquotient br bi))))))
- (defun fpsqrt (x)
- (fproot (bcons x) 2))
- (defun bf-cmod-sl (ar ai)
- (let ((ar (fpabs ar))
- (ai (fpabs ai)))
- (cond ((fpgreaterp ai ar)
- (setq ar (fpquotient ar ai))
- (fptimes* ai
- (fpsqrt (fpplus (fpone) (fptimes* ar ar)))))
- ((fpgreaterp ar ai)
- (setq ai (fpquotient ai ar))
- (fptimes* ar (fpsqrt (fpplus (fpone) (fptimes* ai ai)))))
- ((fptimes* (fpsqrt (intofp 2))
- ar)))))
- (defun bf-calct-sl nil
- (do ((i 1 (1+ i))
- (tt)
- (hvr (setf (aref *qhr-sl* 0) (aref *hr-sl* 0)))
- (hvi (setf (aref *qhi-sl* 0) (aref *hi-sl* 0))))
- ((> i *n*)
- (setq *bool*
- (not (fpgreaterp (bf-cmod-sl hvr hvi)
- (fptimes* (intofp 10)
- (fptimes* *are*
- (bf-cmod-sl (aref *hr-sl* *n*)
- (aref *hi-sl* *n*)))))))
- (cond ((not *bool*)
- (bf-cdivid-sl (fpminus *pvr*) (fpminus *pvi*) hvr hvi)
- (setq *tr* *cr*
- *ti* *ci*))
- (t (setq *tr* (intofp 0) *ti* (intofp 0))))
- nil)
- (setq tt (fpdifference (fpplus (aref *hr-sl* i)
- (fptimes* hvr *sr*))
- (fptimes* hvi *si*)))
- (setf (aref *qhi-sl* i)
- (setq hvi (fpplus (aref *hi-sl* i)
- (fpplus (fptimes* hvr *si*)
- (fptimes* hvi *sr*)))))
- (setf (aref *qhr-sl* i) (setq hvr tt))))
- (defun bf-nexth-sl ()
- (cond (*bool*
- (do ((j 1 (1+ j)))
- ((> j *n*))
- (setf (aref *hr-sl* j) (aref *qhr-sl* (1- j)))
- (setf (aref *hi-sl* j) (aref *qhi-sl* (1- j))))
- (setf (aref *hr-sl* 0) (intofp 0))
- (setf (aref *hi-sl* 0) (intofp 0)))
- (t
- (do ((j 1. (1+ j))
- (t1)
- (t2))
- ((> j *n*))
- (setq t1 (aref *qhr-sl* (1- j))
- t2 (aref *qhi-sl* (1- j)))
- (setf (aref *hr-sl* j)
- (fpdifference (fpplus (aref *qpr-sl* j)
- (fptimes* t1 *tr*))
- (fptimes* t2 *ti*)))
- (setf (aref *hi-sl* j)
- (fpplus (aref *qpi-sl* j)
- (fpplus (fptimes* t1 *ti*)
- (fptimes* t2 *tr*)))))
- (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
- (setf (aref *hi-sl* 0) (aref *qpi-sl* 0))))
- nil)
- (defun bf-polyev-sl ()
- (setq *pvr* (setf (aref *qpr-sl* 0) (aref *pr-sl* 0))
- *pvi* (setf (aref *qpi-sl* 0) (aref *pi-sl* 0)))
- (do ((i 1 (1+ i))
- (tt))
- ((> i *nn*))
- (setq tt (fpdifference (fpplus (aref *pr-sl* i) (fptimes* *pvr* *sr*))
- (fptimes* *pvi* *si*)))
- (setf (aref *qpi-sl* i)
- (setq *pvi* (fpplus (aref *pi-sl* i)
- (fpplus (fptimes* *pvr* *si*)
- (fptimes* *pvi* *sr*)))))
- (setf (aref *qpr-sl* i)
- (setq *pvr* tt))))
- (defun bf-errev-sl (ms mp)
- (fpdifference
- (fptimes* (do ((j 0 (1+ j))
- (e (fpquotient (fptimes* (bf-cmod-sl (aref *qpr-sl* 0) (aref *qpi-sl* 0))
- *mre*)
- (fpplus *are* *mre*))))
- ((> j *nn*) e)
- (setq e (fpplus (bf-cmod-sl (aref *qpr-sl* j) (aref *qpi-sl* j))
- (fptimes* e ms))))
- (fpplus *are* *mre*))
- (fptimes* mp *mre*)))
- (defun bf-cauchy-sl ()
- (let ((x (fpexp (fpquotient (fpdifference (fplog (aref *shr-sl* *nn*))
- (fplog (aref *shr-sl* 0)))
- (intofp *nn*))))
- (xm (intofp 0)))
- (setf (aref *shr-sl* *nn*) (fpminus (aref *shr-sl* *nn*)))
- (cond ((not (fpzerop (aref *shr-sl* *n*)))
- (setq xm (fpminus (fpquotient (aref *shr-sl* *nn*)
- (aref *shr-sl* *n*))))
- (cond ((fpgreaterp x xm) (setq x xm)))))
- (do ((f))
- (nil)
- (setq xm (fptimes* (intofp 0.1) x)
- f (aref *shr-sl* 0))
- (do ((i 1 (1+ i)))
- ((> i *nn*))
- (setq f (fpplus (aref *shr-sl* i)
- (fptimes* f xm))))
- ;;(cond ((not (< 0.0 f)) (return t)))
- (when (fpgreaterp (intofp 0) f)
- (return t))
- (setq x xm))
- (do ((dx x)
- (df)
- (f))
- ((fpgreaterp (intofp 5e-3)
- (fpabs (fpquotient dx x)))
- x)
- (setq f (aref *shr-sl* 0)
- df f)
- (do ((i 1 (1+ i)))
- ((> i *n*))
- (setq f (fpplus (fptimes* f x)
- (aref *shr-sl* i))
- df (fpplus (fptimes* df x)
- f)))
- (setq f (fpplus (fptimes* f x)
- (aref *shr-sl* *nn*))
- dx (fpquotient f df)
- x (fpdifference x dx)))))
- (defun bf-scale-float (bf scale)
- (destructuring-bind (mantissa exp)
- bf
- (if (zerop mantissa)
- (list mantissa 0)
- (list mantissa
- (+ exp scale)))))
- (defun bf-scale-sl ()
- (do ((i 0 (1+ i))
- (j 0)
- (x (intofp 0))
- (dx (intofp 0)))
- ((> i *nn*)
- (setq x (fpquotient x (intofp (- (1+ *nn*) j)))
- dx (fpquotient (fpdifference (fplog (aref *shr-sl* *nn*))
- (fplog (aref *shr-sl* 0)))
- (intofp *nn*))
- *polysc1* (fpentier (bcons (fpplus (cdr bfhalf)
- (fpquotient dx *logbas*))))
- x (fpplus x (fptimes* (intofp (* *polysc1* *nn*))
- (fptimes* *logbas*
- (cdr bfhalf))))
- *polysc* (fpentier (bcons (fpplus (cdr bfhalf) (fpquotient x *logbas*))))))
- (cond ((equalp (aref *shr-sl* i) (intofp 0))
- (setq j (1+ j)))
- (t
- (setq x (fpplus x (fplog (aref *shr-sl* i)))))))
- (do ((i *nn* (1- i))
- (j (- *polysc*) (+ j *polysc1*)))
- ((< i 0))
- (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) j))
- (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) j)))
- nil)
- (defun bf-noshft-sl (l1)
- (do ((i 0 (1+ i))
- (xni (intofp *nn*) (fpdifference xni (intofp 1)))
- (t1 (fpquotient (fpone) (intofp *nn*))))
- ((> i *n*))
- (setf (aref *hr-sl* i) (fptimes* (aref *pr-sl* i)
- (fptimes* xni t1)))
- (setf (aref *hi-sl* i) (fptimes* (aref *pi-sl* i)
- (fptimes* xni t1))))
- (do ((jj 1 (1+ jj)))
- ((> jj l1))
- (cond ((fpgreaterp (bf-cmod-sl (aref *hr-sl* *n*) (aref *hi-sl* *n*))
- (fptimes* (intofp 10)
- (fptimes* *are*
- (bf-cmod-sl (aref *pr-sl* *n*)
- (aref *pi-sl* *n*)))))
- (bf-cdivid-sl (fpminus (aref *pr-sl* *nn*))
- (fpminus (aref *pi-sl* *nn*))
- (aref *hr-sl* *n*)
- (aref *hi-sl* *n*))
- (setq *tr* *cr*
- *ti* *ci*)
- (do ((j *n* (1- j)) (t1) (t2))
- ((> 1 j))
- (setq t1 (aref *hr-sl* (1- j))
- t2 (aref *hi-sl* (1- j)))
- (setf (aref *hr-sl* j) (fpdifference (fpplus (aref *pr-sl* j)
- (fptimes* t1 *tr*))
- (fptimes* t2 *ti*)))
- (setf (aref *hi-sl* j) (fpplus (aref *pi-sl* j)
- (fpplus (fptimes* t1 *ti*)
- (fptimes* t2 *tr*)))))
- (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
- (setf (aref *hi-sl* 0) (aref *pi-sl* 0)))
- (t (do ((j *n* (1- j)))
- ((> 1 j))
- (setf (aref *hr-sl* j) (aref *hr-sl* (1- j)))
- (setf (aref *hi-sl* j) (aref *hi-sl* (1- j))))
- (setf (aref *hr-sl* 0) (intofp 0))
- (setf (aref *hi-sl* 0) (intofp 0))))))
- (defun bf-vrshft-sl (l3)
- (setq *conv* nil
- *sr* *zr*
- *si* *zi*)
- (do ((i 1 (1+ i))
- (bool1 nil)
- (mp)
- (ms)
- (omp)
- (relstp)
- (tp)
- (r1))
- ((> i l3))
- (bf-polyev-sl)
- (setq mp (bf-cmod-sl *pvr* *pvi*)
- ms (bf-cmod-sl *sr* *si*))
- (cond ((fpgreaterp (fptimes* (intofp 20) (bf-errev-sl ms mp))
- mp)
- (setq *conv* t
- *zr* *sr*
- *zi* *si*)
- (return t)))
- (cond ((= i 1)
- (setq omp mp))
- ((or bool1
- (fpgreaterp omp mp)
- ;;(not (< relstp 0.05))
- (fpgreaterp relstp (cdr bfhalf)))
- (if (fpgreaterp (fptimes* (intofp 0.1) mp)
- omp)
- (return t)
- (setq omp mp)))
- (t
- (setq tp relstp
- bool1 t)
- (when (fpgreaterp *are* relstp)
- (setq tp *are*))
- (setq r1 (fpsqrt tp)
- *sr* (prog1
- (fpdifference (fptimes* (fpplus (fpone) r1)
- *sr*)
- (fptimes* r1 *si*))
- (setq *si* (fpplus (fptimes* (fpplus (fpone) r1)
- *si*)
- (fptimes* r1 *sr*)))))
- (bf-polyev-sl)
- (do ((j 1 (1+ j)))
- ((> j 5))
- (bf-calct-sl)
- (bf-nexth-sl))
- (setq omp *infin*)))
- (bf-calct-sl)
- (bf-nexth-sl)
- (bf-calct-sl)
- (or *bool*
- (setq relstp (fpquotient (bf-cmod-sl *tr* *ti*)
- (bf-cmod-sl *sr* *si*))
- *sr* (fpplus *sr* *tr*)
- *si* (fpplus *si* *ti*)))))
- (defun bf-fxshft-sl (l2)
- (let ((test t)
- (pasd nil)
- (otr (intofp 0))
- (oti (intofp 0))
- (svsr (intofp 0))
- (svsi (intofp 0))
- (*bool* nil)
- (*pvr* (intofp 0))
- (*pvi* (intofp 0)))
- (bf-polyev-sl)
- (setq *conv* nil)
- (bf-calct-sl)
- (do ((j 1 (1+ j)))
- ((> j l2))
- (setq otr *tr*
- oti *ti*)
- (bf-nexth-sl)
- (bf-calct-sl)
- (setq *zr* (fpplus *sr* *tr*)
- *zi* (fpplus *si* *ti*))
- (cond ((and (not *bool*)
- test
- (not (= j l2)))
- (cond ((fpgreaterp (fptimes* (cdr bfhalf) (bf-cmod-sl *zr* *zi*))
- (bf-cmod-sl (fpdifference *tr* otr)
- (fpdifference *ti* oti)))
- (cond (pasd
- (do ((i 0 (1+ i)))
- ((> i *n*))
- (setf (aref *shr-sl* i) (aref *hr-sl* i))
- (setf (aref *shi-sl* i) (aref *hi-sl* i)))
- (setq svsr *sr* svsi *si*)
- (bf-vrshft-sl 10.)
- (when *conv* (return nil))
- (setq test nil)
- (do ((i 0 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i) (aref *shr-sl* i))
- (setf (aref *hi-sl* i) (aref *shi-sl* i)))
- (setq *sr* svsr *si* svsi)
- (bf-polyev-sl)
- (bf-calct-sl))
- ((setq pasd t))))
- ((setq pasd nil))))))
- (or *conv* (bf-vrshft-sl 10))
- nil))
- (defun bf-cpoly-sl (degree)
- (let ( ;; Log of our floating point base. (Do we need this much accuracy?)
- (*logbas* (fplog (intofp 2)))
- ;; "Largest" bfloat. What should we use?
- (*infin* (intofp most-positive-flonum))
- ;; bfloat epsilon. 2^(-fpprec)
- (*are* (bf-scale-float (intofp 2) (- fpprec)))
- (*mre* (intofp 0))
- (xx (fproot bfhalf 2))
- (yy (intofp 0))
- ;; cos(94deg). Probably don't need full bfloat precision here.
- (cosr (intofp -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804))
- ;; sin(94deg). Probably don't need full bfloat precision here.
- (sinr (intofp 0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315))
- (*cr* (intofp 0))
- (*ci* (intofp 0))
- (*sr* (intofp 0))
- (*si* (intofp 0))
- (*tr* (intofp 0))
- (*ti* (intofp 0))
- (*zr* (intofp 0))
- (*zi* (intofp 0))
- (bnd (intofp 0))
- (*n* 0)
- (*polysc* 0)
- (*polysc1* 0)
- (*conv* nil))
- (setq *mre* (fptimes* (intofp 2)
- (fptimes* (fpsqrt (intofp 2)) *are*))
- yy (fpminus xx))
- (do ((i degree (1- i)))
- ((not (and (equalp (intofp 0) (aref *pr-sl* i))
- (equalp (intofp 0) (aref *pi-sl* i))))
- (setq *nn* i
- *n* (1- i))))
- (setq degree *nn*)
- (do ((i 0 (1+ i)))
- ((> i *nn*))
- (setf (aref *shr-sl* i) (bf-cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
- (if (> *nn* 0) (bf-scale-sl))
- (do ()
- ((> 2 *nn*)
- (bf-cdivid-sl (fpminus (aref *pr-sl* 1))
- (fpminus (aref *pi-sl* 1))
- (aref *pr-sl* 0)
- (aref *pi-sl* 0))
- (setf (aref *pr-sl* 1) *cr*)
- (setf (aref *pi-sl* 1) *ci*)
- (setq *nn* 0))
- (do ((i 0 (1+ i)))
- ((> i *nn*))
- (setf (aref *shr-sl* i) (bf-cmod-sl (aref *pr-sl* i) (aref *pi-sl* i))))
- (setq bnd (bf-cauchy-sl))
- (catch 'newroot
- (do ((cnt1 1 (1+ cnt1)))
- ((> cnt1 2))
- (bf-noshft-sl 5)
- (do ((cnt2 1 (1+ cnt2)))
- ((> cnt2 9))
- (setq xx (prog1
- (fpdifference (fptimes* cosr xx)
- (fptimes* sinr yy))
- (setq yy (fpplus (fptimes* sinr xx)
- (fptimes* cosr yy))))
- *sr* (fptimes* bnd xx)
- *si* (fptimes* bnd yy))
- (bf-fxshft-sl (* 10 cnt2))
- (cond (*conv* (setf (aref *pr-sl* *nn*) *zr*)
- (setf (aref *pi-sl* *nn*) *zi*)
- (setq *nn* *n* *n* (1- *n*))
- (do ((i 0 (1+ i)))
- ((> i *nn*))
- (setf (aref *pr-sl* i) (aref *qpr-sl* i))
- (setf (aref *pi-sl* i) (aref *qpi-sl* i)))
- (throw 'newroot t))))))
- (or *conv* (return t)))
- (do ((i (1+ *nn*) (1+ i)))
- ((> i degree))
- (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*))
- (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) *polysc1*)))
- (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
- ((> i *nn*))
- (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) j))
- (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) j)))
- *nn*))
- (defmfun $bfallroots (expr)
- (prog (degree *nn* var res $partswitch
- ($keepfloat t)
- $demoivre
- ($listconstvars t)
- ($algebraic t) complex $ratfac den expr1)
- (setq expr1 (setq expr (meqhk expr)))
- (setq var (delete '$%i (cdr ($listofvars expr)) :test #'eq))
- (or var (setq var (list (gensym))))
- (cond ((not (= (length var) 1))
- (merror (intl:gettext "bfallroots: expected a polynomial in one variable; found variables ~M") `((mlist) ,@var)))
- ((setq var (car var))))
- (setq expr ($rat expr '$%i var)
- res (reverse (car (cdddar expr))))
- (do ((i (- (length res)
- (length (caddar expr)))
- (1- i)))
- ((= i 0))
- (setq res (cdr res)))
- (setq den (cddr expr)
- expr (cadr expr))
- ;; Check denominator is a complex number
- (cond ((numberp den) (setq den (list den 0)))
- ((eq (car den) (cadr res))
- (setq den (cddr den))
- (cond ((numberp (car den))
- (cond ((null (cddr den))
- (setq den (list 0 (car den))))
- ((numberp (caddr den))
- (setq den (list (caddr den) (car den))))
- (t (bf-cpoly-err expr1))))
- (t (bf-cpoly-err expr1))))
- (t (bf-cpoly-err expr1)))
- ;; If the name variable has disappeared, this is caught here
- (setq *nn* 0)
- (cond ((numberp expr)
- (setq expr (list expr 0)))
- ((eq (car expr) (car res))
- (setq *nn* 1))
- ((eq (car expr) (cadr res))
- (setq expr (cddr expr))
- (cond ((numberp (car expr))
- (cond ((null (cddr expr))
- (setq expr (list 0 (car expr))))
- ((numberp (caddr expr))
- (setq expr (list (caddr expr) (car expr))))
- (t
- (bf-cpoly-err expr1))))
- (t (bf-cpoly-err expr1))))
- (t (bf-cpoly-err expr1)))
- (cond ((= *nn* 0)
- (cond ($polyfactor
- (let ((*cr* (intofp 0))
- (*ci* (intofp 0)))
- (bf-cdivid-sl (cdr ($bfloat (car expr)))
- (cdr ($bfloat (cadr expr)))
- (cdr ($bfloat (car den)))
- (cdr ($bfloat (cadr den))))
- (return (add (bcons *cr*)
- (mul '$%i (bcons *ci*))))))
- (t (return (list '(mlist simp)))))))
- (setq degree (cadr expr) *nn* (1+ degree))
- (setq *pr-sl* (make-array *nn* :initial-element (intofp 0)))
- (setq *pi-sl* (make-array *nn* :initial-element (intofp 0)))
- (or (catch 'notpoly
- (errset (do ((expr (cdr expr) (cddr expr)) (l) (%i (cadr res)))
- ((null expr))
- (setq l (- degree (car expr)) res (cadr expr))
- (cond ((numberp res)
- (setf (aref *pr-sl* l) (cdr ($bfloat res))))
- (t
- (or (eq (car res) %i)
- (throw 'notpoly nil))
- (setq res (cddr res))
- (setf (aref *pi-sl* l) (cdr ($bfloat (car res))))
- (setq res (caddr res))
- (and res (setf (aref *pr-sl* l) (cdr ($bfloat res))))
- (setq complex t))))))
- ;; This should catch expressions like sin(x)-x
- (bf-cpoly-err expr1))
- (setq *shr-sl* (make-array *nn* :initial-element (intofp 0)))
- (setq *shi-sl* (make-array *nn* :initial-element (intofp 0)))
- (setq *qpr-sl* (make-array *nn* :initial-element (intofp 0)))
- (setq *hr-sl* (make-array degree :initial-element (intofp 0)))
- (setq *qhr-sl* (make-array degree :initial-element (intofp 0)))
- (setq *qpi-sl* (make-array *nn* :initial-element (intofp 0)))
- (when complex
- (setq *hi-sl* (make-array degree :initial-element (intofp 0)))
- (setq *qhi-sl* (make-array degree :initial-element (intofp 0))))
- (setq *nn* degree)
- (if complex
- (setq res (errset (bf-cpoly-sl degree)))
- (setq res (bf-rpoly-sl degree)))
- (unless res
- (mtell (intl:gettext "bfallroots: unexpected error; treat results with caution.~%")))
- (when (= *nn* degree)
- (merror (intl:gettext "bfallroots: no roots found.")))
- (setq res nil)
- (cond ((not (zerop *nn*))
- (mtell (intl:gettext "bfallroots: only ~S out of ~S roots found.~%") (- degree *nn*) degree)
- (setq expr (bcons (intofp 0)))
- (do ((i 0 (1+ i)))
- ((> i *nn*))
- (setq expr
- (simplify
- (list '(mplus) expr
- (simplify (list '(mtimes)
- (simplify (list '(mplus)
- (simplify (list '(mtimes) '$%i
- (bcons (aref *pi-sl* i))))
- (bcons (aref *pr-sl* i))))
- (simplify (list '(mexpt) var (- *nn* i)))))))))
- (setq res (cons expr res)))
- ($polyfactor
- (setq expr (let ((*cr* (intofp 0))
- (*ci* (intofp 0)))
- (bf-cdivid-sl (aref *pr-sl* 0)
- (aref *pi-sl* 0)
- (cdr ($bfloat (car den)))
- (cdr ($bfloat (cadr den))))
- (add (bcons *cr*) (mul '$%i (bcons *ci*))))
- res (cons expr res))))
- (do ((i degree (1- i)))
- ((= i *nn*))
- ;; zr+%i*zi, where zr and zi parts of the root.
- (setq expr (add (bcons (aref *pr-sl* i))
- (mul '$%i (bcons (aref *pi-sl* i)))))
- (setq res
- (cond ($polyfactor
- (cons (cond ((or complex (fpzerop (aref *pi-sl* i)))
- (add var (mul -1 expr)))
- (t
- (setq i (1- i))
- (simplify (list '(mplus)
- (simplify (list '(mexpt) var 2))
- (simplify (list '(mtimes) var
- (bcons (aref *pr-sl* i))))
- (bcons (aref *pr-sl* (1+ i)))))))
- res))
- (t
- (cons (let ((expr (simplify (list '(mequal) var expr))))
- (if $programmode expr (displine expr)))
- res)))))
- (return (simplify (if $polyfactor
- (cons '(mtimes) res)
- (cons '(mlist) (nreverse res)))))))
- (defun bf-rpoly-sl (degree)
- (let ((*logbas* (fplog (intofp 2)))
- (*infin* (intofp most-positive-flonum))
- (*are* (bf-scale-float (intofp 2) (- fpprec)))
- (*mre* (intofp 0))
- (xx (fproot bfhalf 2))
- (yy (intofp 0))
- ;; cos(94deg)
- (cosr (intofp
- -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804))
- ;; sin(94deg)
- (sinr (intofp
- 0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315))
- (aa (intofp 0))
- (cc (intofp 0))
- (bb (intofp 0))
- (bnd (intofp 0))
- (*sr* (intofp 0))
- (*u* (intofp 0))
- (*v* (intofp 0))
- (t1 (intofp 0))
- (*szr* (intofp 0))
- (*szi* (intofp 0))
- (*lzr* (intofp 0))
- (*lzi* (intofp 0))
- (*nz* 0)
- (*n* 0)
- (*polysc* 0)
- (*polysc1* 0)
- (zerok 0)
- (conv1 t))
- (setq *mre* *are*
- yy (fpminus xx))
- (do ((i degree (1- i)))
- ((not (fpzerop (aref *pr-sl* i)))
- (setq *nn* i *n* (1- i))))
- (setq degree *nn*)
- (do ((i 0 (1+ i)))
- ((> i *nn*))
- (setf (aref *shr-sl* i) (fpabs (aref *pr-sl* i))))
- (if (> *nn* 0) (bf-scale-sl))
- (do nil
- ((< *nn* 3)
- (cond ((= *nn* 2)
- (bf-quad-sl (aref *pr-sl* 0.) (aref *pr-sl* 1) (aref *pr-sl* 2))
- (cond ((and $polyfactor (not (fpzerop *szi*)))
- (setf (aref *pr-sl* 2) (fpquotient (aref *pr-sl* 2)
- (aref *pr-sl* 0)))
- (setf (aref *pr-sl* 1) (fpquotient (aref *pr-sl* 1)
- (aref *pr-sl* 0)))
- (setf (aref *pi-sl* 2) (intofp 1)))
- (t (setf (aref *pr-sl* 2) *szr*)
- (setf (aref *pi-sl* 2) *szi*)
- (setf (aref *pr-sl* 1) *lzr*)
- (setf (aref *pi-sl* 1) *lzi*))))
- (t
- (setf (aref *pr-sl* 1) (fpminus (fpquotient (aref *pr-sl* 1)
- (aref *pr-sl* 0))))))
- (setq *nn* 0))
- (do ((i 0 (1+ i)))
- ((> i *nn*))
- (setf (aref *shr-sl* i) (fpabs (aref *pr-sl* i))))
- (setq bnd (bf-cauchy-sl))
- (do ((i 1 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (fpquotient (fptimes* (intofp (- *n* i))
- (aref *pr-sl* i))
- (intofp *n*))))
- (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
- (setq aa (aref *pr-sl* *nn*)
- bb (aref *pr-sl* *n*)
- zerok (fpzerop (aref *hr-sl* *n*)))
- (do ((jj 1 (1+ jj)))
- ((> jj 5.))
- (setq cc (aref *hr-sl* *n*))
- (cond (zerok (do ((j *n* (1- j)))
- ((< j 1))
- (setf (aref *hr-sl* j) (aref *hr-sl* (1- j))))
- (setf (aref *hr-sl* 0) (intofp 0))
- (setq zerok (fpzerop (aref *hr-sl* *n*))))
- (t
- (setq t1 (fpminus (fpquotient aa cc)))
- (do ((j *n* (1- j)))
- ((< j 1))
- (setf (aref *hr-sl* j)
- (fpplus (fptimes* t1 (aref *hr-sl* (1- j)))
- (aref *pr-sl* j))))
- (setf (aref *hr-sl* 0) (aref *pr-sl* 0))
- (setq zerok (not (fpgreaterp (fpabs (aref *hr-sl* *n*))
- (fptimes* (fpabs bb)
- (fptimes* *are* (intofp 10)))))))))
- (do ((i 0 (1+ i)))
- ((> i *n*))
- (setf (aref *shi-sl* i) (aref *hr-sl* i)))
- (do ((cnt 1 (1+ cnt)))
- ((> cnt 20.)
- (setq conv1 nil))
- (setq xx (prog1
- (fpdifference (fptimes* cosr xx)
- (fptimes* sinr yy))
- (setq yy (fpplus (fptimes* sinr xx)
- (fptimes* cosr yy))))
- *sr* (fptimes* bnd xx)
- *u* (fptimes* (intofp -2) *sr*)
- *v* bnd)
- (bf-fxshfr-sl (* 20 cnt))
- (cond ((> *nz* 0)
- (setf (aref *pr-sl* *nn*) *szr*)
- (setf (aref *pi-sl* *nn*) *szi*)
- (cond ((= *nz* 2)
- (setf (aref *pr-sl* *n*) *lzr*)
- (setf (aref *pi-sl* *n*) *lzi*)
- (cond ((and $polyfactor (not (fpzerop *szi*)))
- (setf (aref *pr-sl* *nn*) *v*)
- (setf (aref *pr-sl* *n*) *u*)
- (setf (aref *pi-sl* *nn*) (intofp 1))))))
- (setq *nn* (- *nn* *nz*) *n* (1- *nn*))
- (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *pr-sl* i) (aref *qpr-sl* i)))
- (return nil)))
- (do ((i 0 (1+ i))) ((> i *n*)) (setf (aref *hr-sl* i) (aref *shi-sl* i))))
- (or conv1 (return nil)))
- (cond ($polyfactor
- (do ((i degree (1- i)))
- ((= i *nn*))
- (cond ((fpzerop (aref *pi-sl* i))
- (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*)))
- (t (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) (* 2 *polysc1*)))
- (setq i (1- i))
- (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*))))))
- (t (do ((i (1+ *nn*) (1+ i)))
- ((> i degree))
- (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) *polysc1*))
- (setf (aref *pi-sl* i) (bf-scale-float (aref *pi-sl* i) *polysc1*)))))
- (do ((i 0 (1+ i)) (j (- *polysc* (* *polysc1* degree)) (+ j *polysc1*)))
- ((> i *nn*))
- (setf (aref *pr-sl* i) (bf-scale-float (aref *pr-sl* i) j)))
- t))
- (defun bf-fxshfr-sl (l2)
- (let ((*my-type* 0)
- (*a* (intofp 0))
- (*b* (intofp 0))
- (*c* (intofp 0))
- (*d* (intofp 0))
- (*e* (intofp 0))
- (*f* (intofp 0))
- (*g* (intofp 0))
- (*h* (intofp 0))
- (*a1* (intofp 0))
- (*a3* (intofp 0))
- (*a7* (intofp 0)))
- (declare (special *my-type*))
- (setq *nz* 0)
- (bf-quadsd-sl)
- (bf-calcsc-sl)
- (do ((j 1 (1+ j))
- (betav (intofp 0.25))
- (betas (intofp 0.25))
- (oss *sr*)
- (ovv *v*)
- (tvv) (tss) (ss) (vv) (tv) (ts) (ots) (otv)
- (*ui*) (*vi*) (*s*) (svv) (svu) (iflag) (vpass) (spass) (vtry) (stry))
- ((> j l2))
- (bf-nextk-sl)
- (bf-calcsc-sl)
- (bf-newest-sl)
- (setq vv *vi*
- ss (intofp 0))
- (or (fpzerop (aref *hr-sl* *n*))
- (setq ss (fpminus (fpquotient (aref *pr-sl* *nn*)
- (aref *hr-sl* *n*)))))
- (setq tv (intofp 1)
- ts (intofp 1))
- (cond ((not (or (= j 1)
- (= *my-type* 3)))
- (or (fpzerop vv)
- (setq tv (fpabs (fpquotient (fpdifference vv ovv)
- vv))))
- (or (fpzerop ss)
- (setq ts (fpabs (fpquotient (fpdifference ss oss)
- ss))))
- (setq tvv (intofp 1))
- (and (not (fpgreaterp tv otv))
- (setq tvv (fptimes* tv otv)))
- (setq tss (intofp 1))
- (and (not (fpgreaterp ts ots))
- (setq tss (fptimes* ts ots)))
- (setq vpass (not (fpgreaterp tvv betav))
- spass (not (fpgreaterp tss betas)))
- (cond ((or spass vpass)
- (setq svu *u* svv *v*)
- (do ((i 0 (1+ i)))
- ((> i *n*))
- (setf (aref *shr-sl* i)
- (aref *hr-sl* i)))
- (setq *s* ss vtry nil stry nil)
- (and (do ((bool (not (and spass (or (not vpass)
- (not (fpgreaterp tss tvv)))))
- t)
- (l50 nil nil))
- (nil)
- (cond (bool
- (bf-quadit-sl)
- (and (> *nz* 0) (return t))
- (setq vtry t
- betav (fptimes* (intofp 0.25) betav))
- (cond ((or stry (not spass))
- (setq l50 t))
- (t (do ((i 0 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (aref *shr-sl* i)))))))
- (cond ((not l50)
- (setq iflag (bf-realit-sl))
- (and (> *nz* 0) (return t))
- (setq stry t
- betas (fptimes* (intofp 0.25) betas))
- (cond ((zerop iflag)
- (setq l50 t))
- (t
- (setq *ui* (fpminus (fpplus *s* *s*))
- *vi* (fptimes* *s* *s*))))))
- (cond (l50
- (setq *u* svu *v* svv)
- (do ((i 0 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (aref *shr-sl* i)))
- (and (or (not vpass) vtry)
- (return nil)))))
- (return nil))
- (bf-quadsd-sl)
- (bf-calcsc-sl)))))
- (setq ovv vv
- oss ss
- otv tv
- ots ts))))
- (defun bf-quadit-sl nil
- (setq *nz* 0 *u* *ui* *v* *vi*)
- (do ((tried)
- (j 0)
- (ee)
- (mp)
- (relstp)
- (omp)
- (ms))
- (nil)
- (bf-quad-sl (intofp 1) *u* *v*)
- (and (fpgreaterp (fpabs (fpdifference (fpabs *szr*)
- (fpabs *lzr*)))
- (fptimes* (intofp 1e-2) (fpabs *lzr*)))
- (return nil))
- (bf-quadsd-sl)
- (setq mp (fpplus (fpabs (fpdifference *a* (fptimes* *szr* *b*)))
- (fpabs (fptimes* *szi* *b*)))
- ms (bf-cmod-sl *szr* *szi*)
- ee (bf-errev-sl ms mp))
- (cond ((not (fpgreaterp mp (fptimes* (intofp 2e1) ee)))
- (setq *nz* 2)
- (return nil)))
- (setq j (1+ j))
- (and (> j 20) (return nil))
- (cond ((not (or (< j 2)
- (fpgreaterp relstp (intofp 1e-2))
- (not (fpgreaterp mp omp))
- tried))
- (and (not (fpgreaterp relstp *are*))
- (setq relstp *are*))
- (setq relstp (fpsqrt relstp)
- *u* (fpdifference *u* (fptimes* *u* relstp))
- *v* (fpplus *v* (fptimes* *v* relstp)))
- (bf-quadsd-sl)
- (do ((i 1 (1+ i)))
- ((> i 5))
- (bf-calcsc-sl)
- (bf-nextk-sl))
- (setq tried t j 0)))
- (setq omp mp)
- (bf-calcsc-sl)
- (bf-nextk-sl)
- (bf-calcsc-sl)
- (bf-newest-sl)
- (and (fpzerop *vi*) (return nil))
- (setq relstp (fpabs (fpquotient (fpdifference *vi* *v*)
- *vi*))
- *u* *ui*
- *v* *vi*)))
- (defun bf-realit-sl ()
- (setq *nz* 0)
- (do ((j 0)
- (pv)
- (ee)
- (ms)
- (mp)
- (kv)
- (t1)
- (omp))
- (nil)
- (setq pv (aref *pr-sl* 0))
- (setf (aref *qpr-sl* 0) pv)
- (do ((i 1 (1+ i)))
- ((> i *nn*))
- (setq pv (fpplus (fptimes* pv *s*)
- (aref *pr-sl* i)))
- (setf (aref *qpr-sl* i) pv))
- (setq mp (fpabs pv)
- ms (fpabs *s*)
- ee (fptimes* (fpquotient *mre* (fpplus *are* *mre*))
- (fpabs (aref *qpr-sl* 0))))
- (do ((i 1 (1+ i)))
- ((> i *nn*))
- (setq ee (fpplus (fptimes* ee ms)
- (fpabs (aref *qpr-sl* i)))))
- (cond ((not (fpgreaterp mp
- (fptimes* (intofp 2e1)
- (fpdifference (fptimes* (fpplus *are* *mre*)
- ee)
- (fptimes* *mre* mp)))))
- (setq *nz* 1 *szr* *s* *szi* (intofp 0))
- (return 0)))
- (setq j (1+ j))
- (and (> j 10) (return 0))
- (cond ((not (or (< j 2)
- (fpgreaterp (fpabs t1)
- (fptimes* (intofp 1e-3) (fpabs (fpdifference *s* t1))))
- (not (fpgreaterp mp omp))))
- (return 1)))
- (setq omp mp kv (aref *hr-sl* 0))
- (setf (aref *qhr-sl* 0) kv)
- (do ((i 1 (1+ i)))
- ((> i *n*))
- (setq kv (fpplus (fptimes* kv *s*)
- (aref *hr-sl* i)))
- (setf (aref *qhr-sl* i) kv))
- (cond ((fpgreaterp (fpabs kv)
- (fptimes* (fpabs (aref *hr-sl* *n*))
- (fptimes* (intofp 1e1) *are*)))
- (setq t1 (fpminus (fpquotient pv kv)))
- (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
- (do ((i 1 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (fpplus (fptimes* t1 (aref *qhr-sl* (1- i)))
- (aref *qpr-sl* i)))))
- (t
- (setf (aref *hr-sl* 0) (intofp 0))
- (do ((i 1 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i) (aref *qhr-sl* (1- i))))))
- (setq kv (aref *hr-sl* 0))
- (do ((i 1 (1+ i)))
- ((> i *n*))
- (setq kv (fpplus (fptimes* kv *s*)
- (aref *hr-sl* i))))
- (setq t1 (intofp 0))
- (and (fpgreaterp (fpabs kv)
- (fptimes* (fpabs (aref *hr-sl* *n*))
- (fptimes* (intofp 10) *are*)))
- (setq t1 (fpminus (fpquotient pv kv))))
- (setq *s* (fpplus *s* t1))))
- (defun bf-calcsc-sl ()
- (declare (special *my-type*))
- (setq *d* (aref *hr-sl* 0))
- (setf (aref *qhr-sl* 0) *d*)
- (setq *c* (fpdifference (aref *hr-sl* 1)
- (fptimes* *u* *d*)))
- (setf (aref *qhr-sl* 1) *c*)
- (do ((i 2 (1+ i))
- (c0))
- ((> i *n*))
- (setq c0 (fpdifference (fpdifference (aref *hr-sl* i)
- (fptimes* *u* *c*))
- (fptimes* *v* *d*)))
- (setf (aref *qhr-sl* i) c0)
- (setq *d* *c* *c* c0))
- (cond ((not (or (fpgreaterp (fpabs *c*)
- (fptimes* (fpabs (aref *hr-sl* *n*))
- (fptimes* (intofp 100) *are*)))
- (fpgreaterp (fpabs *d*)
- (fptimes* (fpabs (aref *hr-sl* (1- *n*)))
- (fptimes* (intofp 100) *are*)))))
- (setq *my-type* 3))
- ((not (not (fpgreaterp (fpabs *d*) (fpabs *c*))))
- (setq *my-type* 2
- *e* (fpquotient *a* *d*)
- *f* (fpquotient *c* *d*)
- *g* (fptimes* *u* *b*)
- *h* (fptimes* *v* *b*)
- *a3* (fpplus (fptimes* (fpplus *a* *g*) *e*)
- (fptimes* *h* (fpquotient *b* *d*)))
- *a1* (fpdifference (fptimes* *b* *f*)
- *a*)
- *a7* (fpplus (fptimes* (fpplus *f* *u*) *a*)
- *h*)))
- (t (setq *my-type* 1
- *e* (fpquotient *a* *c*)
- *f* (fpquotient *d* *c*)
- *g* (fptimes* *u* *e*)
- *h* (fptimes* *v* *b*)
- *a3* (fpplus (fptimes* *a* *e*)
- (fptimes* (fpplus (fpquotient *h* *c*)
- *g*)
- *b*))
- *a1* (fpdifference *b*
- (fptimes* *a* (fpquotient *d* *c*)))
- *a7* (fpplus *a*
- (fpplus (fptimes* *g* *d*)
- (fptimes* *h* *f*))))))
- nil)
- (defun bf-nextk-sl ()
- (declare (special *my-type*))
- (cond ((= *my-type* 3)
- (setf (aref *hr-sl* 0) (intofp 0))
- (setf (aref *hr-sl* 1) (intofp 0))
- (do ((i 2 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i) (aref *qhr-sl* (- i 2)))))
- ((fpgreaterp (fpabs *a1*)
- (fptimes* (fpabs (if (= *my-type* 1) *b* *a*))
- (fptimes* (intofp 1e1) *are*)))
- (setq *a7* (fpquotient *a7* *a1*)
- *a3* (fpquotient *a3* *a1*))
- (setf (aref *hr-sl* 0) (aref *qpr-sl* 0))
- (setf (aref *hr-sl* 1)
- (fpdifference (aref *qpr-sl* 1)
- (fptimes* *a7* (aref *qpr-sl* 0))))
- (do ((i 2 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (fpplus (fptimes* *a3* (aref *qhr-sl* (- i 2)))
- (fpplus (fpminus (fptimes* *a7* (aref *qpr-sl* (1- i))))
- (aref *qpr-sl* i))))))
- (t
- (setf (aref *hr-sl* 0) (intofp 0))
- (setf (aref *hr-sl* 1) (fpminus (fptimes* *a7* (aref *qpr-sl* 0))))
- (do ((i 2 (1+ i)))
- ((> i *n*))
- (setf (aref *hr-sl* i)
- (fpdifference (fptimes* *a3* (aref *qhr-sl* (- i 2)))
- (fptimes* *a7* (aref *qpr-sl* (1- i))))))))
- nil)
- (defun bf-newest-sl ()
- (declare (special *my-type*))
- (let ((a4 (intofp 0))
- (a5 (intofp 0))
- (b1 (intofp 0))
- (b2 (intofp 0))
- (c1 (intofp 0))
- (c2 (intofp 0))
- (c3 (intofp 0))
- (c4 (intofp 0)))
- (cond ((= *my-type* 3)
- (setq *ui* (intofp 0)
- *vi* (intofp 0)))
- (t
- (if (= *my-type* 2)
- (setq a4 (fpplus (fptimes* (fpplus *a* *g*)
- *f*)
- *h*)
- a5 (fpplus (fptimes* (fpplus *f* *u*)
- *c*)
- (fptimes* *v* *d*)))
- (setq a4 (fpplus *a*
- (fpplus (fptimes* *u* *b*)
- (fptimes* *h* *f*)))
- a5 (fpplus *c*
- (fptimes* (fpplus *u* (fptimes* *v* *f*))
- *d*))))
- (setq b1 (fpminus (fpquotient (aref *hr-sl* *n*)
- (aref *pr-sl* *nn*)))
- b2 (fpminus (fpquotient (fpplus (aref *hr-sl* (1- *n*))
- (fptimes* b1 (aref *pr-sl* *n*)))
- (aref *pr-sl* *nn*)))
- c1 (fptimes* (fptimes* *v* b2) *a1*)
- c2 (fptimes* b1 *a7*)
- c3 (fptimes* (fptimes* b1 b1) *a3*)
- c4 (fpdifference (fpdifference c1 c2) c3)
- c1 (fpplus (fpplus a5 (fptimes* b1 a4))
- (fpminus c4)))
- (if (fpzerop c1)
- (setq *ui* (intofp 0)
- *vi* (intofp 0))
- (setq *ui* (fpdifference
- *u*
- (fpquotient (fpplus (fptimes* *u* (fpplus c3 c2))
- (fptimes* *v*
- (fpplus (fptimes* b1 *a1*)
- (fptimes* b2 *a7*))))
- c1))
- *vi* (fptimes* *v*
- (fpplus (fpone) (fpquotient c4 c1)))))))
- nil))
- (defun bf-quadsd-sl ()
- (setq *b* (aref *pr-sl* 0))
- (setf (aref *qpr-sl* 0) *b*)
- (setq *a* (fpdifference (aref *pr-sl* 1)
- (fptimes* *u* *b*)))
- (setf (aref *qpr-sl* 1) *a*)
- (do ((i 2 (1+ i))
- (c0))
- ((> i *nn*))
- (setq c0 (fpdifference (fpdifference (aref *pr-sl* i)
- (fptimes* *u* *a*))
- (fptimes* *v* *b*)))
- (setf (aref *qpr-sl* i) c0)
- (setq *b* *a*
- *a* c0)))
- (defun bf-quad-sl (a0 b1 c0)
- (setq *szr* (intofp 0)
- *szi* (intofp 0)
- *lzr* (intofp 0)
- *lzi* (intofp 0))
- (let ((b0 (intofp 0))
- (l0 (intofp 0))
- (*e* (intofp 0)))
- (cond ((fpzerop a0)
- (or (fpzerop b1)
- (setq *szr* (fpminus (fpquotient c0 b1)))))
- ((fpzerop c0)
- (setq *lzr* (fpminus (fpquotient b1 a0))))
- (t
- (setq b0 (fpquotient b1 (intofp 2)))
- (cond ((not (fpgreaterp (fpabs b0) (fpabs c0)))
- (setq *e* a0)
- (and (not (fpgreaterp c0 (intofp 0)))
- (setq *e* (fpminus a0)))
- (setq *e* (fpdifference (fptimes* b0 (fpquotient b0 (fpabs c0)))
- *e*)
- l0 (fptimes* (fpsqrt (fpabs *e*))
- (fpsqrt (fpabs c0)))))
- (t (setq *e* (fpdifference (intofp 1)
- (fptimes* (fpquotient a0 b0)
- (fpquotient c0 b0)))
- l0 (fptimes* (fpsqrt (fpabs *e*))
- (fpabs b0)))))
- (cond ((not (fpgreaterp *e* (intofp 0)))
- (setq *szr* (fpminus (fpquotient b0 a0))
- *lzr* *szr*
- *szi* (fpabs (fpquotient l0 a0))
- *lzi* (fpminus *szi*)))
- (t (or (not (fpgreaterp b0 (intofp 0)))
- (setq l0 (fpminus l0)))
- (setq *lzr* (fpquotient (fpdifference l0 b0) a0))
- (or (fpzerop *lzr*)
- (setq *szr* (fpquotient (fpquotient c0 *lzr*) a0)))))))
- nil))