1
0
mirror of https://github.com/PDP-10/its.git synced 2026-02-14 03:54:00 +00:00
Files
PDP-10.its/src/cffk/cpoly.64
Eric Swenson 85994ed770 Added files to support building and running Macsyma.
Resolves #284.

Commented out uses of time-origin in maxtul; mcldmp (init) until we
can figure out why it gives arithmetic overflows under the emulators.

Updated the expect script statements in build_macsyma_portion to not
attempt to match expected strings, but simply sleep for some time
since in some cases the matching appears not to work.
2018-03-11 13:10:19 -07:00

902 lines
30 KiB
Common Lisp
Raw Blame History

This file contains invisible Unicode characters
This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
;;;;;;;;;;;;;;;;;;; -*- Mode: Lisp; Package: Macsyma -*- ;;;;;;;;;;;;;;;;;;;
;;; (c) Copyright 1981 Massachusetts Institute of Technology ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(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\ and pi\, rather than
;;;into their own arrays. The algorithm seems to benefit be taking are
;;;mre 0.01 times the published values.
(declare
(*expr displa $listofvars meqhk displine)
(special logbas infin smalno are mre cr ci sr si tr ti zr zi n nn bool conv pvr
pvi $partswitch $keepfloat $demoivre $listconstvars $algebraic acp\
$polyfactor polysc polysc1 $ratfac $programmode)
(flonum logbas infin smalno are mre cr ci sr si tr ti zr zi xx yy cosr sinr bnd
xni t1 t2 otr oti svsr svsi pvr pvi mp ms omp relstp tp hvr hvi e ar
ai br bi x xm f dx df r1 $t hi lo max min acp\)
(fixnum degree n nn j l l1 l2 l3 cnt1 cnt2 jj polysc polysc1))
(declare (notype ($cpoly notype) (noshft\ fixnum) (fxshft\ fixnum)
(vrshft\ fixnum) (calct\) (nexth\) (polyev\)
(cdivid\ flonum flonum flonum flonum) (scale\))
(flonum (errev\ flonum flonum) (cauchy\) (cmod\ flonum flonum))
(fixnum (cpoly\ fixnum))
#-PDP10(flonum (*f flonum flonum) (//f flonum flonum) (_f flonum fixnum))
#-PDP10(*expr *f //f _f)
(*lexpr $rat error))
;; The Lisp Machine needs some way of declaring "functional" arrays.
;; It warns about these symbols being referenced as functions but not defined.
(declare
(array* (flonum pr\ 1. pi\ 1. shr\ 1. shi\ 1. qpr\ 1. qpi\ 1. hr\ 1. hi\
1. qhr\ 1. qhi\ 1.)))
;; Fixed for Twenex systems?
#+PDP10
(and (not (get '*f 'subr))
(mapc '(lambda (x) (putprop x '(arith fasl dsk macsym) 'autoload))
'(*f //f +f -f _f)))
;;; It is harder to underflow on lisp machine, but I suppose someday -BEE
#-PDP10
(progn 'compile
(defmacro *f (a b) `(*$ ,a ,b))
(defmacro //f (a b) `(//$ ,a ,b))
(defmacro +f (a b) `(+$ ,a ,b))
(defmacro -f (a b) `(-$ ,a ,b))
#+Multics(defmacro _f (a b) `(fsc ,a ,b))
)
(defmacro cpoly-large-flonum ()
#-LISPM '(fsc (lsh -1 -1) 0.)
#+LISPM (let ((a (float 0)))
(%p-dpb -1 1013 a) (%p-dpb -1 0007 a)
(%p-dpb-offset -1 0030 a 1) a))
(defmacro cpoly-small-flonum ()
#-LISPM '(fsc (rot 1. -10.) 0.)
#+LISPM (%p-dpb 100 0007 (float 0)))
(defmacro float-precision (pres)
pres ;Ignored on Lisp Machine
#-LISPM `(*$ (fsc (lsh 205. 26.) 0.) ,pres)
#+LISPM (let ((a (float 1)))
(%p-dpb-offset 1 0001 a 1) (- a 1.0)))
;; #+Franz
;; (defun _f (number scale) somebody-write-this)
#+LISPM
(defun _f (number scale)
(let ((ans (+ 0.0 number))
(exp))
(setq exp (+ scale (%p-ldb 1013 ans)))
(cond ((zerop number) 0.0)
((> exp 3777) (ferror nil "_F Overflow -- see MC:CFFK;CPOLY"))
;; Should check zunderflow
((< exp 0) (ferror nil "_F Underflow -- see MC:CFFK;CPOLY"))
(t (%p-dpb exp 1013 ans) ans))))
(setq acp\ 0.2)
(DEFMFUN $allroots (expr)
(prog (degree nn var res $partswitch $keepfloat $demoivre $listconstvars
$algebraic complex $ratfac den)
(setq $keepfloat t $listconstvars t $algebraic t)
(setq expr (meqhk expr) var (delq '$%i (cdr ($listofvars expr))))
(or var (setq var (list (gensym))))
(cond ((not (= (length var) 1.))
(merror "polynomial not univariate: ~M" 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 (error '|not a polynomial|))))
(t (error '|not a polynomial|))))
(t (error '|not a polynomial|)))
;;;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 (error '|not a polynomial|))))
(t (error '|not a polynomial|))))
(t (error '|not a polynomial|)))
(cond ((= nn 0)
(cond ($polyfactor
((lambda (cr ci)
(cdivid\ (float (car expr)) (float (cadr expr))
(float (car den)) (float (cadr den)))
(return (simplify (list '(mplus)
(simplify (list '(mtimes)
'$%i ci))
cr))))
0.0 0.0))
(t (return (list '(mlist simp)))))))
(setq degree (cadr expr) nn (1+ degree))
(array pr\ flonum nn) #+LISPM (fillarray 'pr\ '(0.0))
(array pi\ flonum nn) #+LISPM (fillarray 'pi\ '(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) (store (pr\ l) (float res)))
(t (or (eq (car res) %i) (*throw 'notpoly nil))
(setq res (cddr res))
(store (pi\ l) (float (car res)))
(setq res (caddr res))
(and res (store (pr\ l) (float res)))
(setq complex t))))))
;;;this should catch expressions like sin(x)-x
(progn (*rearray 'pr\)
(*rearray 'pi\)
(error '|not a polynomial|)))
(array shr\ flonum nn) #+LISPM (fillarray 'shr\ '(0.0))
(array shi\ flonum nn) #+LISPM (fillarray 'shi\ '(0.0))
(array qpr\ flonum nn) #+LISPM (fillarray 'qpr\ '(0.0))
(array hr\ flonum degree) #+LISPM (fillarray 'hr\ '(0.0))
(array qhr\ flonum degree) #+LISPM (fillarray 'qhr\ '(0.0))
(cond (complex (array qpi\ flonum nn)
#+LISPM (fillarray 'qpi\ '(0.0))
(array hi\ flonum degree)
#+LISPM (fillarray 'hi\ '(0.0))
(array qhi\ flonum degree)
#+LISPM (fillarray 'qhi\ '(0.0))))
(setq nn degree)
(cond (complex (setq res (errset (cpoly\ degree))))
((setq res (errset (rpoly\ degree)))))
(*rearray 'shr\)
(*rearray 'shi\)
(*rearray 'qpr\)
(*rearray 'hr\)
(*rearray 'qhr\)
(cond (complex (*rearray 'qpi\)
(*rearray 'hi\)
(*rearray 'qhi\)))
(or res
(mtell "~%Unexpected error. Treat results with caution."))
(cond ((= nn degree)
(*rearray 'pr\)
(*rearray 'pi\)
(merror "~%No roots found")))
(setq res nil)
(cond
((not (= nn 0.))
(mtell "~%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
(pi\ i)))
(pr\ i)))
(simplify (list '(mexpt)
var
(- nn i)))))))))
(setq res (cons expr res)))
($polyfactor
(setq expr ((lambda (cr ci)
(cdivid\ (pr\ 0) (pi\ 0)
(float (car den))
(float (cadr den)))
(simplify (list '(mplus)
(simplify (list '(mtimes)
'$%i ci))
cr)))
0.0 0.0)
res (cons expr res))))
(do
((i degree (1- i)))
((= i nn))
(setq expr (simplify (list '(mplus)
(simplify (list '(mtimes)
'$%i
(pi\ i)))
(pr\ i))))
(setq
res
(cond
($polyfactor (cons (cond ((or complex (= (pi\ i) 0.0))
(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
(pr\ i)))
(pr\ (1+ i))))))
res))
((cons ((lambda (expr) (cond ($programmode expr)
(t (displine expr))))
(simplify (list '(mequal) var expr)))
res)))))
(*rearray 'pr\)
(*rearray 'pi\)
(return (simplify (cond ($polyfactor (cons '(mtimes) res))
((cons '(mlist) (nreverse res))))))))
(defun cpoly\ (degree)
((lambda (logbas infin smalno are mre xx yy cosr sinr cr ci sr si tr ti zr zi bnd
n polysc polysc1 conv)
(setq mre (*$ 2.0 (sqrt 2.0) are) yy (-$ xx))
(do ((i degree (1- i)))
((not (and (= (pr\ i) 0.0) (= (pi\ i) 0.0))) (setq nn i n (1- i))))
(setq degree nn)
(do ((i 0. (1+ i)))
((> i nn))
(store (shr\ i) (cmod\ (pr\ i) (pi\ i))))
(scale\)
(do nil
((> 2. nn)
(cdivid\ (-$ (pr\ 1.)) (-$ (pi\ 1.)) (pr\ 0.) (pi\ 0.))
(store (pr\ 1.) cr)
(store (pi\ 1.) ci)
(setq nn 0.))
(do ((i 0. (1+ i)))
((> i nn))
(store (shr\ i) (cmod\ (pr\ i) (pi\ i))))
(setq bnd (cauchy\))
(*catch 'newroot
(do ((cnt1 1. (1+ cnt1)))
((> cnt1 2.))
(noshft\ 5.)
(do ((cnt2 1. (1+ cnt2)))
((> cnt2 9.))
(setq xx (prog2 nil
(-$ (*$ cosr xx) (*$ sinr yy))
(setq yy (+$ (*$ sinr xx)
(*$ cosr yy))))
sr (*$ bnd xx)
si (*$ bnd yy))
(fxshft\ (* 10. cnt2))
(cond (conv (store (pr\ nn) zr)
(store (pi\ nn) zi)
(setq nn n n (1- n))
(do ((i 0. (1+ i)))
((> i nn))
(store (pr\ i) (qpr\ i))
(store (pi\ i) (qpi\ i)))
(*throw 'newroot t))))))
(or conv (return t)))
(do ((i (1+ nn) (1+ i)))
((> i degree))
(store (pr\ i) (_f (pr\ i) polysc1))
(store (pi\ i) (_f (pi\ i) polysc1)))
(do ((i 0. (1+ i)) (j (- polysc (* polysc1 degree)) (+ j polysc1)))
((> i nn))
(store (pr\ i) (_f (pr\ i) j))
(store (pi\ i) (_f (pi\ i) j)))
nn)
(log 2.0) (cpoly-large-flonum)
(cpoly-small-flonum) (float-precision acp\)
0.0 0.70710677 0.0 -0.069756474 0.99756405
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0. 0. 0. nil))
(defun noshft\ (l1)
(do ((i 0. (1+ i)) (xni (float nn) (1-$ xni)) (t1 (//$ (float nn))))
((> i n))
(store (hr\ i) (*$ (pr\ i) xni t1))
(store (hi\ i) (*$ (pi\ i) xni t1)))
(do ((jj 1. (1+ jj)))
((> jj l1))
(cond ((> (cmod\ (hr\ n) (hi\ n)) (*$ 10.0 are (cmod\ (pr\ n) (pi\ n))))
(cdivid\ (-$ (pr\ nn)) (-$ (pi\ nn)) (hr\ n) (hi\ n))
(setq tr cr ti ci)
(do ((j n (1- j)) (t1) (t2))
((> 1. j))
(setq t1 (hr\ (1- j)) t2 (hi\ (1- j)))
(store (hr\ j) (-$ (+$ (pr\ j) (*f t1 tr)) (*f t2 ti)))
(store (hi\ j) (+$ (pi\ j) (*f t1 ti) (*f t2 tr))))
(store (hr\ 0.) (pr\ 0.))
(store (hi\ 0.) (pi\ 0.)))
(t (do ((j n (1- j)))
((> 1. j))
(store (hr\ j) (hr\ (1- j)))
(store (hi\ j) (hi\ (1- j))))
(store (hr\ 0.) 0.0)
(store (hi\ 0.) 0.0)))))
(defun fxshft\ (l2)
((lambda (test pasd otr oti svsr svsi bool pvr pvi)
(polyev\)
(setq conv nil)
(calct\)
(do ((j 1. (1+ j)))
((> j l2))
(setq otr tr oti ti)
(nexth\)
(calct\)
(setq zr (+$ sr tr) zi (+$ si ti))
(cond ((and (not bool) test (not (= j l2)))
(cond ((> (*$ 0.5 (cmod\ zr zi))
(cmod\ (-$ tr otr) (-$ ti oti)))
(cond (pasd (do ((i 0. (1+ i)))
((> i n))
(store (shr\ i) (hr\ i))
(store (shi\ i) (hi\ i)))
(setq svsr sr svsi si)
(vrshft\ 10.)
(and conv (return nil))
(setq test nil)
(do ((i 0. (1+ i)))
((> i n))
(store (hr\ i) (shr\ i))
(store (hi\ i) (shi\ i)))
(setq sr svsr si svsi)
(polyev\)
(calct\))
((setq pasd t))))
((setq pasd nil))))))
(or conv (vrshft\ 10.))
nil)
t nil 0.0 0.0 0.0 0.0 nil 0.0 0.0))
(defun vrshft\ (l3)
(setq conv nil sr zr si zi)
(do ((i 1. (1+ i)) (bool1 nil) (mp) (ms) (omp) (relstp) (tp) (r1))
((> i l3))
(polyev\)
(setq mp (cmod\ pvr pvi) ms (cmod\ sr si))
(cond ((> (*$ 20.0 (errev\ 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)))
(cond ((> (*$ 0.1 mp) omp) (return t)) (t (setq omp mp))))
(t (setq tp relstp bool1 t)
(cond ((> are relstp) (setq tp are)))
(setq r1 (sqrt tp)
sr (prog2 nil
(-$ (*$ (1+$ r1) sr) (*f r1 si))
(setq si (+$ (*$ (1+$ r1) si) (*f r1 sr)))))
(polyev\)
(do ((j 1. (1+ j))) ((> j 5.)) (calct\) (nexth\))
(setq omp infin)))
(calct\)
(nexth\)
(calct\)
(or bool
(setq relstp (//$ (cmod\ tr ti) (cmod\ sr si))
sr (+$ sr tr)
si (+$ si ti)))))
(defun calct\ nil
(do ((i 1. (1+ i))
($t)
(hvr (store (qhr\ 0.) (hr\ 0.)))
(hvi (store (qhi\ 0.) (hi\ 0.))))
((> i n)
(setq bool (not (> (cmod\ hvr hvi) (*$ 10.0 are (cmod\ (hr\ n) (hi\ n))))))
(cond ((not bool) (cdivid\ (-$ pvr) (-$ pvi) hvr hvi) (setq tr cr ti ci))
(t (setq tr 0.0 ti 0.0)))
nil)
(setq $t (-$ (+$ (hr\ i) (*f hvr sr)) (*f hvi si)))
(store (qhi\ i) (setq hvi (+$ (hi\ i) (*f hvr si) (*f hvi sr))))
(store (qhr\ i) (setq hvr $t))))
(defun nexth\ nil
(cond (bool (do ((j 1. (1+ j)))
((> j n))
(store (hr\ j) (qhr\ (1- j)))
(store (hi\ j) (qhi\ (1- j))))
(store (hr\ 0.) 0.0)
(store (hi\ 0.) 0.0))
(t (do ((j 1. (1+ j)) (t1) (t2))
((> j n))
(setq t1 (qhr\ (1- j)) t2 (qhi\ (1- j)))
(store (hr\ j) (-$ (+$ (qpr\ j) (*f t1 tr)) (*f t2 ti)))
(store (hi\ j) (+$ (qpi\ j) (*f t1 ti) (*f t2 tr))))
(store (hr\ 0.) (qpr\ 0.))
(store (hi\ 0.) (qpi\ 0.))))
nil)
(defun polyev\ nil
(setq pvr (store (qpr\ 0.) (pr\ 0.)) pvi (store (qpi\ 0.) (pi\ 0.)))
(do ((i 1. (1+ i)) ($t))
((> i nn))
(setq $t (-$ (+$ (pr\ i) (*f pvr sr)) (*f pvi si)))
(store (qpi\ i) (setq pvi (+$ (pi\ i) (*f pvr si) (*f pvi sr))))
(store (qpr\ i) (setq pvr $t))))
(defun errev\ (ms mp)
(-$ (*$ (do ((j 0. (1+ j))
(e (//$ (*$ (cmod\ (qpr\ 0.) (qpi\ 0.)) mre) (+$ are mre))))
((> j nn) e)
(setq e (+$ (cmod\ (qpr\ j) (qpi\ j)) (*$ e ms))))
(+$ are mre))
(*$ mp mre)))
(defun cauchy\ nil
((lambda (x xm)
(store (shr\ nn) (-$ (shr\ nn)))
(cond ((not (= (shr\ n) 0.0))
(setq xm (-$ (//$ (shr\ nn) (shr\ n))))
(cond ((> x xm) (setq x xm)))))
(do ((f))
(nil)
(setq xm (*$ 0.1 x) f (shr\ 0.))
(do ((i 1. (1+ i))) ((> i nn)) (setq f (+$ (shr\ i) (*f f xm))))
(cond ((not (< 0.0 f)) (return t)))
(setq x xm))
(do ((dx x) (df) (f))
((> 5.0e-3 (abs (//$ dx x))) x)
(setq f (shr\ 0.) df f)
(do ((i 1. (1+ i)))
((> i n))
(setq f (+$ (*$ f x) (shr\ i)) df (+$ (*$ df x) f)))
(setq f (+$ (*$ f x) (shr\ nn)) dx (//$ f df) x (-$ x dx))))
(exp (//$ (-$ (log (shr\ nn)) (log (shr\ 0.))) (float nn)))
0.0))
(defun scale\ nil
(do ((i 0. (1+ i)) (j 0.) (x 0.0) (dx 0.0))
((> i nn)
(setq x (//$ x (float (- (1+ nn) j)))
dx (//$ (-$ (log (shr\ nn)) (log (shr\ 0.))) (float nn))
polysc1 (fix (+$ 0.5 (//$ dx logbas)))
x (+$ x (*$ (float (* polysc1 nn)) logbas 0.5))
polysc (fix (+$ 0.5 (//$ x logbas)))))
(cond ((= (shr\ i) 0.0) (setq j (1+ j)))
(t (setq x (+$ x (log (shr\ i)))))))
(do ((i nn (1- i)) (j (- polysc) (+ j polysc1)))
((< i 0.))
(store (pr\ i) (_f (pr\ i) j))
(store (pi\ i) (_f (pi\ i) j))))
;; (defun scale\ nil
;; ((lambda (hi lo max min x l)
;; (do ((i 0. (1+ i)))
;; ((> i nn))
;; (setq x (shr\ i))
;; (cond ((> x max) (setq max x)))
;; (cond ((and (not (= x 0.0)) (< x min)) (setq min x))))
;; (cond ((or (> lo min) (> max hi))
;; (setq x (//$ lo min))
;; (cond ((> x 1.0)
;; (cond ((> max (//$ infin x))
;; ;;;acm has < here but imsl agrees with me
;; (setq x 1.0))))
;; ((setq x (//$ (*$ (sqrt max) (sqrt min))))))
;; (setq l (fix (+$ 0.5 (//$ (log x) logbas))))
;; (cond ((not (= l 0.))
;; (do ((i 0. (1+ i)))
;; ((> i nn))
;; (store (pr\ i) (_f (pr\ i) l))
;; (store (pi\ i) (_f (pi\ i) l)))))))
;; l)
;; (sqrt infin) (//$ smalno are) 0.0 infin 0.0 0.))
(defun cdivid\ (ar ai br bi)
((lambda (r1) (cond ((and (= br 0.0) (= bi 0.0)) (setq cr (setq ci infin)))
((> (abs bi) (abs br))
(setq r1 (//f br bi)
bi (+$ bi (*f br r1))
br (+$ ai (*f ar r1))
cr (//f br bi)
br (-$ (*f ai r1) ar)
ci (//f br bi)))
((setq r1 (//f bi br)
bi (+$ br (*f bi r1))
br (+$ ar (*f ai r1))
cr (//f br bi)
br (-$ ai (*f ar r1))
ci (//f br bi)))))
0.0)
nil)
(defun cmod\ (ar ai)
(setq ar (abs ar) ai (abs ai))
(cond ((> ai ar) (setq ar (//f ar ai)) (*$ ai (sqrt (1+$ (*f ar ar)))))
((> ar ai) (setq ai (//f ai ar)) (*$ ar (sqrt (1+$ (*f ai ai)))))
((*$ 1.41421357 ar))))
;;*page
;;;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\, qp - qpr\, k - hr\, qk - qhr\, svk - shr\,
;;;temp - shi\. the roots are put in pr\ and pi\.
;;;the variable si appears not to be used here
(declare (special sr u v a b c d a1 a3 a7 e f g h szr szi lzr lzi are mre n nn nz
type ui vi s $polyfactor arp\)
(flonum a a0 a1 a3 a4 a5 a6 a7 aa are b b0 b1 b2 logbas bb betas betav bnd c c0
c1 c2 c3 c4 cc cosr d d0 e ee f g h infin kv lzi lzr mp mre ms omp
oss ots otv ovv pv relstp s sinr smalno sr ss svu svv szi szr t1 ts
tss tv tvv u ui v vi vv xx yy zm arp\)
(fixnum cnt degree i iflag j jj l l2 n nn nz type))
(declare (fixnum (realit\))
(notype (rpoly\ fixnum) (fxshfr\ fixnum) (quadit\) (calcsc\) (nextk\)
(newest\) (quadsd\) (quad\ flonum flonum flonum)))
(setq arp\ 1.0)
(defun rpoly\ (degree)
((lambda (logbas infin smalno are mre xx yy cosr sinr aa cc bb bnd sr u v t1 szr
szi lzr lzi nz n polysc polysc1 zerok conv1)
(setq mre are yy (-$ xx))
(do ((i degree (1- i))) ((not (= (pr\ i) 0.0)) (setq nn i n (1- i))))
(setq degree nn)
(do ((i 0. (1+ i))) ((> i nn)) (store (shr\ i) (abs (pr\ i))))
(scale\)
(do nil
((< nn 3.)
(cond ((= nn 2.)
(quad\ (pr\ 0.) (pr\ 1.) (pr\ 2.))
(cond ((and $polyfactor (not (= szi 0.0)))
(store (pr\ 2.) (//$ (pr\ 2.) (pr\ 0.)))
(store (pr\ 1.) (//$ (pr\ 1.) (pr\ 0.)))
(store (pi\ 2.) 1.0))
(t (store (pr\ 2.) szr)
(store (pi\ 2.) szi)
(store (pr\ 1.) lzr)
(store (pi\ 1.) lzi))))
(t (store (pr\ 1.) (-$ (//$ (pr\ 1.) (pr\ 0.))))))
(setq nn 0.))
(do ((i 0. (1+ i))) ((> i nn)) (store (shr\ i) (abs (pr\ i))))
(setq bnd (cauchy\))
(do ((i 1. (1+ i)))
((> i n))
(store (hr\ i) (//$ (*$ (float (- n i)) (pr\ i)) (float n))))
(store (hr\ 0.) (pr\ 0.))
(setq aa (pr\ nn) bb (pr\ n) zerok (= (hr\ n) 0.0))
(do ((jj 1. (1+ jj)))
((> jj 5.))
(setq cc (hr\ n))
(cond (zerok (do ((j n (1- j)))
((< j 1.))
(store (hr\ j) (hr\ (1- j))))
(store (hr\ 0.) 0.0)
(setq zerok (= (hr\ n) 0.0)))
(t (setq t1 (-$ (//$ aa cc)))
(do ((j n (1- j)))
((< j 1.))
(store (hr\ j) (+$ (*$ t1 (hr\ (1- j))) (pr\ j))))
(store (hr\ 0.) (pr\ 0.))
(setq zerok (not (> (abs (hr\ n))
(*$ (abs bb) are 10.0)))))))
(do ((i 0. (1+ i))) ((> i n)) (store (shi\ i) (hr\ i)))
(do ((cnt 1. (1+ cnt)))
((> cnt 20.) (setq conv1 nil))
(setq xx (prog2 nil
(-$ (*$ cosr xx) (*$ sinr yy))
(setq yy (+$ (*$ sinr xx) (*$ cosr yy))))
sr (*$ bnd xx)
u (*$ -2.0 sr)
v bnd)
(fxshfr\ (* 20. cnt))
(cond ((> nz 0.)
(store (pr\ nn) szr)
(store (pi\ nn) szi)
(cond ((= nz 2.)
(store (pr\ n) lzr)
(store (pi\ n) lzi)
(cond ((and $polyfactor (not (= szi 0.0)))
(store (pr\ nn) v)
(store (pr\ n) u)
(store (pi\ nn) 1.0)))))
(setq nn (- nn nz) n (1- nn))
(do ((i 0. (1+ i))) ((> i nn)) (store (pr\ i) (qpr\ i)))
(return nil)))
(do ((i 0. (1+ i))) ((> i n)) (store (hr\ i) (shi\ i))))
(or conv1 (return nil)))
(cond ($polyfactor
(do ((i degree (1- i)))
((= i nn))
(cond ((= (pi\ i) 0.0)
(store (pr\ i) (_f (pr\ i) polysc1)))
(t (store (pr\ i) (_f (pr\ i) (* 2. polysc1)))
(setq i (1- i))
(store (pr\ i) (_f (pr\ i) polysc1))))))
(t (do ((i (1+ nn) (1+ i)))
((> i degree))
(store (pr\ i) (_f (pr\ i) polysc1))
(store (pi\ i) (_f (pi\ i) polysc1)))))
(do ((i 0. (1+ i)) (j (- polysc (* polysc1 degree)) (+ j polysc1)))
((> i nn))
(store (pr\ i) (_f (pr\ i) j))))
(log 2.0) (cpoly-large-flonum)
(cpoly-small-flonum) (float-precision arp\)
0.0 0.70710677 0.0 -0.069756474 0.99756405
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0. 0. 0. 0. 0. t))
(defun fxshfr\ (l2)
((lambda (type a b c d e f g h a1 a3 a7)
(setq nz 0.)
(quadsd\)
(calcsc\)
(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\)
(calcsc\)
(newest\)
(setq vv vi ss 0.0)
(or (= (hr\ n) 0.0) (setq ss (-$ (//$ (pr\ nn) (hr\ n)))))
(setq tv 1.0 ts 1.0)
(cond ((not (or (= j 1.) (= type 3.)))
(or (= vv 0.0) (setq tv (abs (//$ (-$ vv ovv) vv))))
(or (= ss 0.0) (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)) (store (shr\ i) (hr\ 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\)
(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))
(store (hr\ i)
(shr\ i)))))))
(cond ((not l50)
(setq iflag (realit\))
(and (> nz 0.) (return t))
(setq stry t betas (*$ 0.25 betas))
(cond ((= iflag 0.) (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))
(store (hr\ i) (shr\ i)))
(and (or (not vpass) vtry)
(return nil)))))
(return nil))
(quadsd\)
(calcsc\)))))
(setq ovv vv oss ss otv tv ots ts)))
0. 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0))
(defun quadit\ nil
(setq nz 0. u ui v vi)
(do ((tried) (j 0.) (ee) (zm) (t1) (mp) (relstp) (omp))
(nil)
(quad\ 1.0 u v)
(and (> (abs (-$ (abs szr) (abs lzr))) (*$ 0.01 (abs lzr))) (return nil))
(quadsd\)
(setq mp (+$ (abs (-$ a (*$ szr b))) (abs (*$ szi b)))
zm (sqrt (abs v))
ee (*$ 2.0 (abs (qpr\ 0.)))
t1 (-$ (*$ szr b)))
(do ((i 1. (1+ n))) ((> i n)) (setq ee (+$ (*$ ee zm) (abs (qpr\ i)))))
(setq ee (+$ (*$ ee zm) (abs (+$ a t1)))
ee (-$ (*$ (+$ (*$ 5.0 mre) (*$ 4.0 are)) ee)
(*$ (+$ (*$ 5.0 mre) (*$ 2.0 are))
(+$ (abs (+$ a t1)) (*$ (abs b) zm)))
(*$ -2.0 are (abs t1))))
(cond ((not (> mp (*$ 20.0 ee))) (setq nz 2.) (return nil)))
(setq j (1+ j))
(and (> j 20.) (return nil))
(cond ((not (or (< j 2.) (> relstp 0.01) (< mp omp) tried))
(and (< relstp are) (setq relstp are))
(setq relstp (sqrt relstp)
u (-$ u (*$ u relstp))
v (+$ v (*$ v relstp)))
(quadsd\)
(do ((i 1. (1+ i))) ((> i 5.)) (calcsc\) (nextk\))
(setq tried t j 0.)))
(setq omp mp)
(calcsc\)
(nextk\)
(calcsc\)
(newest\)
(and (= vi 0.0) (return nil))
(setq relstp (abs (//$ (-$ vi v) vi)) u ui v vi)))
(defun realit\ nil
(setq nz 0.)
(do ((j 0.) (pv) (ee) (ms) (mp) (kv) (t1) (omp))
(nil)
(setq pv (pr\ 0.))
(store (qpr\ 0.) pv)
(do ((i 1. (1+ i)))
((> i nn))
(setq pv (+$ (*$ pv s) (pr\ i)))
(store (qpr\ i) pv))
(setq mp (abs pv) ms (abs s) ee (*$ (//$ mre (+$ are mre)) (abs (qpr\ 0.))))
(do ((i 1. (1+ i))) ((> i nn)) (setq ee (+$ (*$ ee ms) (abs (qpr\ i)))))
(cond ((not (> mp (*$ 20.0 (-$ (*$ (+$ 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) (*$ 1.0e-3 (abs (-$ s t1))))
(not (> mp omp))))
(return 1.)))
(setq omp mp kv (hr\ 0.))
(store (qhr\ 0.) kv)
(do ((i 1. (1+ i)))
((> i n))
(setq kv (+$ (*$ kv s) (hr\ i)))
(store (qhr\ i) kv))
(cond ((> (abs kv) (*$ (abs (hr\ n)) 10.0 are))
(setq t1 (-$ (//$ pv kv)))
(store (hr\ 0.) (qpr\ 0.))
(do ((i 1. (1+ i)))
((> i n))
(store (hr\ i) (+$ (*$ t1 (qhr\ (1- i))) (qpr\ i)))))
(t (store (hr\ 0.) 0.0)
(do ((i 1. (1+ i))) ((> i n)) (store (hr\ i) (qhr\ (1- i))))))
(setq kv (hr\ 0.))
(do ((i 1. (1+ i))) ((> i n)) (setq kv (+$ (*$ kv s) (hr\ i))))
(setq t1 0.0)
(and (> (abs kv) (*$ (abs (hr\ n)) 10.0 are)) (setq t1 (-$ (//$ pv kv))))
(setq s (+$ s t1))))
(defun calcsc\ nil
(setq d (hr\ 0.))
(store (qhr\ 0.) d)
(setq c (-$ (hr\ 1.) (*$ u d)))
(store (qhr\ 1.) c)
(do ((i 2. (1+ i)) (c0))
((> i n))
(setq c0 (-$ (hr\ i) (*$ u c) (*$ v d)))
(store (qhr\ i) c0)
(setq d c c c0))
(cond ((not (or (> (abs c) (*$ (abs (hr\ n)) 100.0 are))
(> (abs d) (*$ (abs (hr\ (1- n))) 100.0 are))))
(setq type 3.))
((not (< (abs d) (abs c)))
(setq 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 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\ nil
(cond ((= type 3.)
(store (hr\ 0.) 0.0)
(store (hr\ 1.) 0.0)
(do ((i 2. (1+ i))) ((> i n)) (store (hr\ i) (qhr\ (- i 2.)))))
((> (abs a1) (*$ (abs (cond ((= type 1.) b) (a))) 10.0 are))
(setq a7 (//$ a7 a1) a3 (//$ a3 a1))
(store (hr\ 0.) (qpr\ 0.))
(store (hr\ 1.) (-$ (qpr\ 1.) (*$ a7 (qpr\ 0.))))
(do ((i 2. (1+ i)))
((> i n))
(store (hr\ i)
(+$ (*$ a3 (qhr\ (- i 2.)))
(-$ (*$ a7 (qpr\ (1- i))))
(qpr\ i)))))
(t (store (hr\ 0.) 0.0)
(store (hr\ 1.) (-$ (*$ a7 (qpr\ 0.))))
(do ((i 2. (1+ i)))
((> i n))
(store (hr\ i)
(-$ (*$ a3 (qhr\ (- i 2.))) (*$ a7 (qpr\ (1- i))))))))
nil)
(defun newest\ nil
((lambda (a4 a5 b1 b2 c1 c2 c3 c4)
(cond ((= type 3.) (setq ui 0.0 vi 0.0))
(t (cond ((= type 2.)
(setq a4 (+$ (*$ (+$ a g) f) h)
a5 (+$ (*$ (+$ f u) c) (*$ v d))))
(t (setq a4 (+$ a (*$ u b) (*$ h f))
a5 (+$ c (*$ (+$ u (*$ v f)) d)))))
(setq b1 (-$ (//$ (hr\ n) (pr\ nn)))
b2 (-$ (//$ (+$ (hr\ (1- n)) (*$ b1 (pr\ n))) (pr\ nn)))
c1 (*$ v b2 a1)
c2 (*$ b1 a7)
c3 (*$ b1 b1 a3)
c4 (-$ c1 c2 c3)
c1 (+$ a5 (*$ b1 a4) (-$ c4)))
(cond ((= c1 0.0) (setq ui 0.0 vi 0.0))
(t (setq ui (-$ u
(//$ (+$ (*$ u (+$ c3 c2))
(*$ v
(+$ (*$ b1 a1)
(*$ b2 a7))))
c1))
vi (*$ v (1+$ (//$ c4 c1))))))))
nil)
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0))
(defun quadsd\ nil
(setq b (pr\ 0.))
(store (qpr\ 0.) b)
(setq a (-$ (pr\ 1.) (*$ u b)))
(store (qpr\ 1.) a)
(do ((i 2. (1+ i)) (c0))
((> i nn))
(setq c0 (-$ (pr\ i) (*$ u a) (*$ v b)))
(store (qpr\ i) c0)
(setq b a a c0)))
(defun quad\ (a0 b1 c0)
(setq szr 0.0 szi 0.0 lzr 0.0 lzi 0.0)
((lambda (b0 d0 e)
(cond ((= a0 0.0) (or (= b1 0.0) (setq szr (-$ (//$ c0 b1)))))
((= c0 0.0) (setq lzr (-$ (//$ b1 a0))))
(t (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)
d0 (*$ (sqrt (abs e)) (sqrt (abs c0)))))
(t (setq e (-$ 1.0 (*$ (//$ a0 b0) (//$ c0 b0)))
d0 (*$ (sqrt (abs e)) (abs b0)))))
(cond ((< e 0.0)
(setq szr (-$ (//$ b0 a0))
lzr szr
szi (abs (//$ d0 a0))
lzi (-$ szi)))
(t (or (< b0 0.0) (setq d0 (-$ d0)))
(setq lzr (//$ (-$ d0 b0) a0))
(or (= lzr 0.0) (setq szr (//$ c0 lzr a0)))))))
nil)
0.0 0.0 0.0))
(declare (unspecial logbas infin smalno are mre cr ci sr si tr ti zr zi
n nn bool conv pvr pvi acp\ polysc polysc1 sr u v a
b c d a1 a3 a7 e f g h szr szi lzr lzi are mre n nn nz
type ui vi s arp\))