From: Andreas E. <ar...@us...> - 2009-02-28 15:44:16
|
Update of /cvsroot/maxima/maxima/src In directory 23jxhf1.ch3.sourceforge.com:/tmp/cvs-serv11170/src Modified Files: cpoly.lisp float.lisp laplac.lisp rpart.lisp Log Message: changed lambda binding to let binding for better readability and clarity. Index: cpoly.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/src/cpoly.lisp,v retrieving revision 1.18 retrieving revision 1.19 diff -u -d -r1.18 -r1.19 --- cpoly.lisp 17 Feb 2009 15:39:36 -0000 1.18 +++ cpoly.lisp 28 Feb 2009 15:44:01 -0000 1.19 @@ -33,7 +33,7 @@ (declare-top (special *u* *v* *a* *b* *c* *d* *a1* *a3* *a7* *e* *f* *g* *h* *szr* *szi* *lzr* *lzi* *nz* *ui* *vi* *s*)) - + #+clisp (eval-when (:compile-toplevel :load-toplevel :execute) (ext:without-package-lock () @@ -81,13 +81,12 @@ (t (cpoly-err expr1))) (cond ((= *nn* 0) (cond ($polyfactor - ((lambda (*cr* *ci*) - (cdivid-sl (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)) + (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)) @@ -97,7 +96,7 @@ ((null expr)) (setq l (- degree (car expr)) res (cadr expr)) (cond ((numberp res) (setf (aref *pr-sl* l) (float res))) - (t + (t (or (eq (car res) %i) (throw 'notpoly nil)) (setq res (cddr res)) @@ -140,12 +139,11 @@ (simplify (list '(mexpt) var (- *nn* i))))))))) (setq res (cons expr res))) ($polyfactor - (setq expr ((lambda (*cr* *ci*) - (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*))) - 0.0 0.0) + (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*)) @@ -162,85 +160,88 @@ (aref *pr-sl* i))) (aref *pr-sl* (1+ i)))))) res)) - ((cons ((lambda (expr) (if $programmode expr (displine expr))) - (simplify (list '(mequal) var expr))) + ((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))))))) + (cons '(mlist) (nreverse res))))))) (defun cpoly-err (expr) (merror (intl:gettext "allroots: expected a polynomial; found ~M") expr)) -(defun cpoly-sl (degree) - ((lambda (*logbas* *infin* *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 (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*) - (log 2.0) - most-positive-flonum - flonum-epsilon - 0.0 - ;; sqrt(0.5) - 0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207863675069231154561485124624180279253686063220607485499679157066113329637527963778999752505763910302857350547799858029851372672984310073642587093204445993047761646152421543571607254198813018139976257039948436267 - 0.0 - ;; cos(94deg) - -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804 - ;; sin(94deg) - 0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315 - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0 0 nil)) +(defun cpoly-sl (degree) + (let ((*logbas* (log 2.0)) + (*infin* most-positive-flonum) + (*are* flonum-epsilon) + (*mre* 0.0) + (xx 0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207863675069231154561485124624180279253686063220607485499679157066113329637527963778999752505763910302857350547799858029851372672984310073642587093204445993047761646152421543571607254198813018139976257039948436267) ;; sqrt(0.5) + (yy 0.0) + (cosr -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804) ;; cos(94deg) + (sinr 0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315) ;; sin(94deg) + (*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) +(defun noshft-sl (l1) (do ((i 0 (1+ i)) (xni (float *nn*) (1- xni)) (t1 (/ (float *nn*)))) @@ -265,44 +266,48 @@ (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))))) + (setf (aref *hi-sl* 0) 0.0))))) -(defun fxshft-sl (l2) - ((lambda (test pasd otr oti svsr svsi *bool* *pvr* *pvi*) - (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) - t nil 0.0 0.0 0.0 0.0 nil 0.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) +(defun vrshft-sl (l3) (setq *conv* nil *sr* *zr* *si* *zi*) (do ((i 1 (1+ i)) (bool1 nil) @@ -320,7 +325,7 @@ (setq omp mp))) (t (setq tp relstp bool1 t) (when (> *are* relstp) (setq tp *are*)) - (setq r1 (sqrt tp) + (setq r1 (sqrt tp) *sr* (prog1 (- (* (1+ r1) *sr*) (* r1 *si*)) (setq *si* (+ (* (1+ r1) *si*) (* r1 *sr*))))) @@ -331,11 +336,11 @@ (nexth-sl) (calct-sl) (or *bool* - (setq relstp (/ (cmod-sl *tr* *ti*) (cmod-sl *sr* *si*)) - *sr* (+ *sr* *tr*) - *si* (+ *si* *ti*))))) + (setq relstp (/ (cmod-sl *tr* *ti*) (cmod-sl *sr* *si*)) + *sr* (+ *sr* *tr*) + *si* (+ *si* *ti*))))) -(defun calct-sl nil +(defun calct-sl nil (do ((i 1 (1+ i)) ($t) (hvr (setf (aref *qhr-sl* 0) (aref *hr-sl* 0))) @@ -347,7 +352,7 @@ 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)))) + (setf (aref *qhr-sl* i) (setq hvr $t)))) (defun nexth-sl () (cond (*bool* (do ((j 1 (1+ j))) @@ -363,7 +368,7 @@ (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) + nil) (defun polyev-sl () (setq *pvr* (setf (aref *qpr-sl* 0) (aref *pr-sl* 0)) @@ -372,37 +377,36 @@ ((> 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)))) + (setf (aref *qpr-sl* i) (setq *pvr* $t)))) -(defun errev-sl (ms mp) +(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*))) + (* mp *mre*))) (defun cauchy-sl () - ((lambda (x xm) - (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)))) - (exp (/ (- (log (aref *shr-sl* *nn*)) (log (aref *shr-sl* 0))) (float *nn*))) - 0.0)) + (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))))) (defun scale-sl () (do ((i 0 (1+ i)) (j 0) (x 0.0) (dx 0.0)) @@ -419,26 +423,26 @@ (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) - ((lambda (r1) - (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))))) - 0.0) - nil) +(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) +(defun cmod-sl (ar ai) (setq ar (abs ar) ai (abs ai)) (cond ((> ai ar) (setq ar (/ ar ai)) (* ai (sqrt (1+ (* ar ar))))) @@ -454,176 +458,187 @@ ;;; the roots *are* put in pr-sl and pi-sl. The variable *si* appears ;;; not to be used here -(defun rpoly-sl (degree) - ((lambda (*logbas* *infin* *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 (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)) - (do nil - ((< *nn* 3) - (cond ((= *nn* 2) - (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 (setf (aref *pr-sl* 1) (- (/ (aref *pr-sl* 1) (aref *pr-sl* 0)))))) - (setq *nn* 0)) - (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i)))) - (setq bnd (cauchy-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 ((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)))) - (log 2.0) - most-positive-flonum - flonum-epsilon - 0.0 - ;; sqrt(0.5) - 0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207863675069231154561485124624180279253686063220607485499679157066113329637527963778999752505763910302857350547799858029851372672984310073642587093204445993047761646152421543571607254198813018139976257039948436267 - 0.0 - ;; cos(94deg) - -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804 - ;; sin(94deg) - 0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315 - 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 rpoly-sl (degree) + (let ((*logbas* (log 2.0)) + (*infin* most-positive-flonum) + (*are* flonum-epsilon) + (*mre* 0.0) + (xx 0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207863675069231154561485124624180279253686063220607485499679157066113329637527963778999752505763910302857350547799858029851372672984310073642587093204445993047761646152421543571607254198813018139976257039948436267) ;; sqrt(0.5) + (yy 0.0) + (cosr -0.0697564737441253007759588351941433286009032016527965250436172961370711270667891229125378568280742923028942076107741717160209821557740512756197740925891665208235244345674420755726285778495732000059330205461129612198466216775458241726113210999152981126990497403794217445425671287263223529689424188857433131142804) ;; cos(94deg) + (sinr 0.9975640502598242476131626806442550263693776603842215362259956088982191814110818430852892116754760376426967121558233963175758546629687044962793968705262246733087781690124900795021134880736278349857522534853084644420433826380899280074903993378273609249428279246573946968632240992430211366514177713203298481315) ;; sin(94deg) + (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)) + (do nil + ((< *nn* 3) + (cond ((= *nn* 2) + (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 (setf (aref *pr-sl* 1) (- (/ (aref *pr-sl* 1) (aref *pr-sl* 0)))))) + (setq *nn* 0)) + (do ((i 0 (1+ i))) ((> i *nn*)) (setf (aref *shr-sl* i) (abs (aref *pr-sl* i)))) + (setq bnd (cauchy-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 ((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) - ((lambda (*my-type* *a* *b* *c* *d* *e* *f* *g* *h* *a1* *a3* *a7*) - (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) +(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) + (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))) - 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0)) + (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 +(defun quadit-sl nil (setq *nz* 0 *u* *ui* *v* *vi*) (do ((tried) (j 0) (ee) (zm) (t1) (mp) (relstp) (omp)) (nil) @@ -631,13 +646,13 @@ (and (> (abs (- (abs *szr*) (abs *lzr*))) (* 1e-2 (abs *lzr*))) (return nil)) (quadsd-sl) - (setq mp (+ (abs (- *a* (* *szr* *b*))) (abs (* *szi* *b*))) - zm (sqrt (abs *v*)) + (setq mp (+ (abs (- *a* (* *szr* *b*))) (abs (* *szi* *b*))) + zm (sqrt (abs *v*)) ee (* 2.0 (abs (aref *qpr-sl* 0))) t1 (- (* *szr* *b*))) (do ((i 1 (1+ *n*))) ((> i *n*)) (setq ee (+ (* ee zm) (abs (aref *qpr-sl* i))))) - (setq ee (+ (* ee zm) (abs (+ *a* t1))) + (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))) @@ -648,8 +663,8 @@ (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)) + (setq relstp (sqrt relstp) + *u* (- *u* (* *u* relstp)) *v* (+ *v* (* *v* relstp))) (quadsd-sl) (do ((i 1 (1+ i))) @@ -727,21 +742,21 @@ (> (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*) + (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*))) + *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) @@ -773,28 +788,31 @@ (defun newest-sl () (declare (special *my-type*)) - ((lambda (a4 a5 b1 b2 c1 c2 c3 c4) - (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) - 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0)) + (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)) (defun quadsd-sl () (setq *b* (aref *pr-sl* 0)) @@ -809,29 +827,30 @@ (setq *b* *a* *a* c0))) -(defun quad-sl (a0 b1 c0) +(defun quad-sl (a0 b1 c0) (setq *szr* 0.0 *szi* 0.0 *lzr* 0.0 *lzi* 0.0) - ((lambda (b0 l0 *e*) - (cond ((zerop a0) (or (zerop b1) (setq *szr* (- (/ c0 b1))))) - ((zerop c0) (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*) - 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) - 0.0 0.0 0.0)) + (let ((b0 0.0) + (l0 0.0) + (*e* 0.0)) + (cond ((zerop a0) (or (zerop b1) (setq *szr* (- (/ c0 b1))))) + ((zerop c0) (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*) + 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. @@ -1415,9 +1434,8 @@ (bcons (aref *pr-sl* (1+ i))))))) res)) (t - (cons ((lambda (expr) - (if $programmode expr (displine expr))) - (simplify (list '(mequal) var expr))) + (cons (let ((expr (simplify (list '(mequal) var expr)))) + (if $programmode expr (displine expr))) res))))) (return (simplify (if $polyfactor (cons '(mtimes) res) Index: float.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/src/float.lisp,v retrieving revision 1.51 retrieving revision 1.52 diff -u -d -r1.51 -r1.52 --- float.lisp 27 Feb 2009 02:04:24 -0000 1.51 +++ float.lisp 28 Feb 2009 15:44:01 -0000 1.52 @@ -163,20 +163,20 @@ (list '|0| '|.| '|0| '|b| '|0|)) (t ;; L IS ALWAYS POSITIVE FP NUMBER (let ((extradigs (floor (1+ (quotient (integer-length (caddr l)) #.(/ (log 10.0) (log 2.0)))))) - (*m 1) (*cancelled 0)) + (*m 1) + (*cancelled 0)) (setq l - ((lambda (*decfp fpprec of l expon) - (setq expon (- (cadr l) of)) - (setq l (if (minusp expon) - (fpquotient (intofp (car l)) (fpintexpt 2 (- expon) of)) - (fptimes* (intofp (car l)) (fpintexpt 2 expon of)))) - (incf fpprec (- extradigs)) - (list (fpround (car l)) (+ (- extradigs) *m (cadr l)))) - t - (+ extradigs (decimalsin (- (caddar l) 2))) - (caddar l) - (cdr l) - nil))) + (let ((*decfp t) + (fpprec (+ extradigs (decimalsin (- (caddar l) 2)))) + (of (caddar l)) + (l (cdr l)) + (expon nil)) + (setq expon (- (cadr l) of)) + (setq l (if (minusp expon) + (fpquotient (intofp (car l)) (fpintexpt 2 (- expon) of)) + (fptimes* (intofp (car l)) (fpintexpt 2 expon of)))) + (incf fpprec (- extradigs)) + (list (fpround (car l)) (+ (- extradigs) *m (cadr l)))))) (let ((*print-base* 10.) *print-radix* (l1 nil)) @@ -231,49 +231,52 @@ (defun bigfloat2rat (x) (setq x (bigfloatp x)) - ((lambda ($float2bf exp y sign) - (setq exp (cond ((minusp (cadr x)) - (setq sign t y (fpration1 (cons (car x) (fpabs (cdr x))))) - (rplaca y (* -1 (car y)))) - (t (fpration1 x)))) - (cond ($ratprint (princ "`rat' replaced ") - (when sign (princ "-")) - (princ (maknam (fpformat (cons (car x) (fpabs (cdr x)))))) - (princ " by ") - (princ (car exp)) - (write-char #\/) - (princ (cdr exp)) - (princ " = ") - (setq x ($bfloat (list '(rat simp) (car exp) (cdr exp)))) - (when sign (princ "-")) - (princ (maknam (fpformat (cons (car x) (fpabs (cdr x)))))) - (terpri))) - exp) - t nil nil nil)) + (let (($float2bf t) + (exp nil) + (y nil) + (sign nil)) + (setq exp (cond ((minusp (cadr x)) + (setq sign t + y (fpration1 (cons (car x) (fpabs (cdr x))))) + (rplaca y (* -1 (car y)))) + (t (fpration1 x)))) + (when $ratprint + (princ "`rat' replaced ") + (when sign (princ "-")) + (princ (maknam (fpformat (cons (car x) (fpabs (cdr x)))))) + (princ " by ") + (princ (car exp)) + (write-char #\/) + (princ (cdr exp)) + (princ " = ") + (setq x ($bfloat (list '(rat simp) (car exp) (cdr exp)))) + (when sign (princ "-")) + (princ (maknam (fpformat (cons (car x) (fpabs (cdr x)))))) + (terpri)) + exp)) (defun fpration1 (x) - ((lambda (fprateps) - (or (and (equal x bigfloatzero) (cons 0 1)) - (prog (y a) - (return (do ((xx x (setq y (invertbigfloat - (bcons (fpdifference (cdr xx) (cdr ($bfloat a))))))) - (num (setq a (fpentier x)) - (+ (* (setq a (fpentier y)) num) onum)) - (den 1 (+ (* a den) oden)) - (onum 1 num) - (oden 0 den)) - ((and (not (zerop den)) - (not (fpgreaterp - (fpabs (fpquotient - (fpdifference (cdr x) - (fpquotient (cdr ($bfloat num)) - (cdr ($bfloat den)))) - (cdr x))) - fprateps))) - (cons num den))))))) - (cdr ($bfloat (if $bftorat - (list '(rat simp) 1 (exptrl 2 (1- fpprec))) - $ratepsilon))))) + (let ((fprateps (cdr ($bfloat (if $bftorat + (list '(rat simp) 1 (exptrl 2 (1- fpprec))) + $ratepsilon))))) + (or (and (equal x bigfloatzero) (cons 0 1)) + (prog (y a) + (return (do ((xx x (setq y (invertbigfloat + (bcons (fpdifference (cdr xx) (cdr ($bfloat a))))))) + (num (setq a (fpentier x)) + (+ (* (setq a (fpentier y)) num) onum)) + (den 1 (+ (* a den) oden)) + (onum 1 num) + (oden 0 den)) + ((and (not (zerop den)) + (not (fpgreaterp + (fpabs (fpquotient + (fpdifference (cdr x) + (fpquotient (cdr ($bfloat num)) + (cdr ($bfloat den)))) + (cdr x))) + fprateps))) + (cons num den)))))))) (defun float-nan-p (x) (and (floatp x) (not (= x x)))) @@ -629,16 +632,18 @@ ;; r = floor(x). (defun fpexp (x) (prog (r s) - (if (not (signp ge (car x))) - (return (fpquotient (fpone) (fpexp (fpabs x))))) + (unless (signp ge (car x)) + (return (fpquotient (fpone) (fpexp (fpabs x))))) (setq r (fpintpart x)) - (return (cond ((< r 2) (fpexp1 x)) - (t (setq s (fpexp1 (fpdifference x (intofp r)))) - (fptimes* s - (cdr (bigfloatp - ((lambda (fpprec r) (bcons (fpexpt (fpe) r))) ; patch for full precision %E - (+ fpprec (integer-length r) -1) - r))))))))) + (return (cond ((< r 2) + (fpexp1 x)) + (t + (setq s (fpexp1 (fpdifference x (intofp r)))) + (fptimes* s + (cdr (bigfloatp + (let ((fpprec (+ fpprec (integer-length r) -1)) + (r r)) + (bcons (fpexpt (fpe) r))))))))))) ; patch for full precision %E ;; exp(x) for small x, using Taylor series. (defun fpexp1 (x) @@ -722,26 +727,26 @@ (defun fpe1 nil (bcons (list (fpround (compe (+ fpprec 12))) (+ -12 *m)))) ;; -;; compe is the bignum part of the bfloat(%e) computation +;; compe is the bignum part of the bfloat(%e) computation ;; (compe N)/(2.0^N) is an approximation to E ;; The algorithm is based on the series ;; ;; %e = sum( 1/i! ,i,0,inf ) -;; +;; ;; but sums up 1001 terms to one. -;; +;; (defun compe (prec) (let (s h (n 1) d (k 1001)) (setq h (ash 1 prec)) (setq s h) (do ((i k (+ i k))) - ((zerop h)) + ((zerop h)) (setq d (do ((j 1 (1+ j)) (p i)) - ((> j (1- k)) (* p n)) - (setq p (* p (- i j)))) ) + ((> j (1- k)) (* p n)) + (setq p (* p (- i j)))) ) (setq n (do ((j (- k 2) (1- j)) (p 1)) - ((< j 0) p) - (setq p (1+ (* p (- i j))))) ) + ((< j 0) p) + (setq p (1+ (* p (- i j))))) ) (setq h (*quo (* h n) d)) (setq s (+ s h))) s)) @@ -753,33 +758,33 @@ ;; fppi1 is the bigfloat part of the bfloat(%pi) computation ;; (defun fppi1 nil - (bcons - (fpquotient + (bcons + (fpquotient (fprt18231_) (list (fpround (comppi (+ fpprec 12))) (+ -12 *m)) ))) ;; -;; comppi is the bignum part of the bfloat(%pi) computation +;; comppi is the bignum part of the bfloat(%pi) computation ;; (comppi N)/(2.0^N) is an approximation to 640320^(3/2)/12 * 1/PI ;; ;; Chudnovsky & Chudnovsky (1987): ;; -;; 640320^(3/2) / (12 * %pi) = +;; 640320^(3/2) / (12 * %pi) = ;; ;; sum( (-1)^i*(6*i)!*(545140134*i+13591409) / (i!^3*(3*i)!*640320^(3*i)) ,i,0,inf ) -;; +;; (defun comppi (prec) - (let (s h n d) + (let (s h n d) (setq s (ash 13591409 prec)) (setq h (neg (*quo (ash 67047785160 prec) 262537412640768000))) (setq s (+ s h)) (do ((i 2 (1+ i))) - ((zerop h)) + ((zerop h)) (setq n (* 12 (- (* 6 i) 5) (- (* 6 i) 4) (- (* 2 i) 1) (- (* 6 i) 1) (+ (* i 545140134) 13591409) )) (setq d (* (- (* 3 i) 2) (expt i 3) (- (* i 545140134) 531548725) 262537412640768000)) (setq h (neg (*quo (* h n) d))) (setq s (+ s h))) s )) -;; +;; ;; fprt18231_ computes sqrt(640320^3/12^2) ;; = sqrt(1823176476672000) = 42698670.666333... ;; @@ -796,11 +801,11 @@ ;; quadratic Heron algorithm: x[0] = ----, , sqrt(a) = ------ ;; d[0] d[i+1] = 2*n[i]*d[i] d[inf] #+gcl -(defun fprt18231_ () +(defun fprt18231_ () (let ((a 1823176476672000) - (n 42698670666) - (d 1000) - h ) + (n 42698670666) + (d 1000) + h ) (do ((prec 32 (* 2 prec))) ((> prec fpprec)) (setq h n) @@ -999,8 +1004,9 @@ ((and (< (cadr p) 0) (ratnump n)) ($bfloat ($expand (list '(mtimes) - ($bfloat ((lambda ($domain $m1pbranch) (power -1 n)) - '$complex t)) + ($bfloat (let (($domain '$complex) + ($m1pbranch t)) + (power -1 n))) (exptbigfloat (bcons (fpminus (cdr p))) n))))) ((and (< (cadr p) 0) (not (integerp n))) (cond ((or (equal n 0.5) (equal n bfhalf)) @@ -1024,16 +1030,13 @@ (cond ((not ($bfloatp n)) (list '(mexpt) p n)) (t - ((lambda (extrabits) - (setq p - ((lambda (fpprec) - (fpexp (fptimes* (cdr (bigfloatp n)) (fplog (cdr (bigfloatp p)))))) - (+ extrabits fpprec))) - (setq p (list (fpround (car p)) - (+ (- extrabits) *m (cadr p)))) - (bcons p)) - (max 1 (+ (caddr n) (integer-length (caddr p)))))))) - ; The number of extra bits required + (let ((extrabits (max 1 (+ (caddr n) (integer-length (caddr p)))))) + (setq p + (let ((fpprec (+ extrabits fpprec))) + (fpexp (fptimes* (cdr (bigfloatp n)) (fplog (cdr (bigfloatp p))))))) + (setq p (list (fpround (car p)) (+ (- extrabits) *m (cadr p)))) + (bcons p))))) + ;; The number of extra bits required ((< n 0) (invertbigfloat (exptbigfloat p (- n)))) (t (bcons (fpexpt (cdr p) n))))) @@ -1142,35 +1145,35 @@ (return (cdr (bigfloatp - ((lambda (fpprec xt *cancelled oldprec) - (prog (x) - loop (setq x (cdr (bigfloatp xt))) - (setq piby2 (fpquotient (fppi) (intofp 2))) - (setq r (fpintpart (fpquotient x piby2))) - (setq x (fpplus x (fptimes* (intofp (- r)) piby2))) - (setq k *cancelled) - (fpplus x (fpminus piby2)) - (setq *cancelled (max k *cancelled)) - (when *fpsincheck* - (print `(*canc= ,*cancelled fpprec= ,fpprec oldprec= ,oldprec))) - (cond ((not (> oldprec (- fpprec *cancelled))) - (setq r (rem r 4)) - (setq res - (cond (fl (cond ((= r 0) (fpsin1 x)) - ((= r 1) (fpcos1 x)) - ((= r 2) (fpminus (fpsin1 x))) - ((= r 3) (fpminus (fpcos1 x))))) - (t (cond ((= r 0) (fpcos1 x)) - ((= r 1) (fpminus (fpsin1 x))) - ((= r 2) (fpminus (fpcos1 x))) - ((= r 3) (fpsin1 x)))))) - (return (bcons (if sign res (fpminus res))))) - (t (incf fpprec *cancelled) - (go loop))))) - (max fpprec (+ fpprec (cadr x))) - (bcons x) - 0 - fpprec)))))) + (let ((fpprec (max fpprec (+ fpprec (cadr x)))) + (xt (bcons x)) + (*cancelled 0) + (oldprec fpprec)) + (prog (x) + loop (setq x (cdr (bigfloatp xt))) + (setq piby2 (fpquotient (fppi) (intofp 2))) + (setq r (fpintpart (fpquotient x piby2))) + (setq x (fpplus x (fptimes* (intofp (- r)) piby2))) + (setq k *cancelled) + (fpplus x (fpminus piby2)) + (setq *cancelled (max k *cancelled)) + (when *fpsincheck* + (print `(*canc= ,*cancelled fpprec= ,fpprec oldprec= ,oldprec))) + (cond ((not (> oldprec (- fpprec *cancelled))) + (setq r (rem r 4)) + (setq res + (cond (fl (cond ((= r 0) (fpsin1 x)) + ((= r 1) (fpcos1 x)) + ((= r 2) (fpminus (fpsin1 x))) + ((= r 3) (fpminus (fpcos1 x))))) + (t (cond ((= r 0) (fpcos1 x)) + ((= r 1) (fpminus (fpsin1 x))) + ((= r 2) (fpminus (fpcos1 x))) + ((= r 3) (fpsin1 x)))))) + (return (bcons (if sign res (fpminus res))))) + (t + (incf fpprec *cancelled) + (go loop)))))))))) (defun fpcos1 (x) (fpsincos1 x nil)) @@ -1210,25 +1213,24 @@ (* (car f) (expt 2 (- m))))))) (defun logbigfloat (a) - ((lambda (minus) - (setq a ((lambda (fpprec) - (cond (($bfloatp (car a)) - (setq a ($bfloat (car a))) - (cond ((zerop (cadr a)) (merror "log(0.0b0) has been generated")) - ((minusp (cadr a)) - (setq minus t) (fplog (list (- (cadr a)) (caddr a)))) - (t (fplog (cdr a))))) - (t (list '(%log) (car a))))) - (+ 2 fpprec))) - (when (numberp (car a)) - (setq a (if (zerop (car a)) - (list 0 0) - (list (fpround (car a)) (+ -2 *m (cadr a))))) - (setq a (bcons a))) - (if minus - (add a (mul '$%i ($bfloat '$%pi))) - a)) - nil)) + (let ((minus nil)) + (setq a (let ((fpprec (+ 2 fpprec))) + (cond (($bfloatp (car a)) + (setq a ($bfloat (car a))) + (cond ((zerop (cadr a)) (merror "log(0.0b0) has been generated")) + ((minusp (cadr a)) + (setq minus t) (fplog (list (- (cadr a)) (caddr a)))) + (t (fplog (cdr a))))) + (t + (list '(%log) (car a)))))) + (when (numberp (car a)) + (setq a (if (zerop (car a)) + (list 0 0) + (list (fpround (car a)) (+ -2 *m (cadr a))))) + (setq a (bcons a))) + (if minus + (add a (mul '$%i ($bfloat '$%pi))) + a))) ;;; Computes the log of a bigfloat number. ;;; Index: laplac.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/src/laplac.lisp,v retrieving revision 1.10 retrieving revision 1.11 diff -u -d -r1.10 -r1.11 --- laplac.lisp 12 Jan 2009 21:34:08 -0000 1.10 +++ laplac.lisp 28 Feb 2009 15:44:01 -0000 1.11 @@ -134,24 +134,24 @@ (simptimes (list '(mtimes) -1 (sdiff (laptimes (cdr fun)) parm)) 1 t)) - (t ((lambda (op) - (cond ((eq op 'mexpt) - (lapexpt (car fun) (cdr fun))) - ((eq op 'mplus) - (laplus ($multthru (fixuprest (cdr fun)) (car fun)))) - ((eq op '%sin) - (lapsin (car fun) (cdr fun) nil)) - ((eq op '%cos) - (lapsin (car fun) (cdr fun) t)) - ((eq op '%sinh) - (lapsinh (car fun) (cdr fun) nil)) - ((eq op '%cosh) - (lapsinh (car fun) (cdr fun) t)) - ((eq op '$delta) - (lapdelta (car fun) (cdr fun))) - - (t (lapshift (car fun) (cdr fun))))) - (caaar fun))))) + (t + (let ((op (caaar fun))) + (cond ((eq op 'mexpt) + (lapexpt (car fun) (cdr fun))) + ((eq op 'mplus) + (laplus ($multthru (fixuprest (cdr fun)) (car fun)))) + ((eq op '%sin) + (lapsin (car fun) (cdr fun) nil)) + ((eq op '%cos) + (lapsin (car fun) (cdr fun) t)) + ((eq op '%sinh) + (lapsinh (car fun) (cdr fun) nil)) + ((eq op '%cosh) + (lapsinh (car fun) (cdr fun) t)) + ((eq op '$delta) + (lapdelta (car fun) (cdr fun))) + (t + (lapshift (car fun) (cdr fun)))))))) (defun lapexpt (fun rest) ;;;HANDLES %E**(A*T+B)*REST(T), %E**(A*T**2+B*T+C), @@ -295,16 +295,12 @@ rest)))) (t (lapshift fun rest)))))) +;;;INTEGRAL FROM A TO INFINITY OF F(X) (defun mydefint (f x a) - ;;;INTEGRAL FROM A TO INFINITY OF F(X) - ((lambda (tryint) (cond (tryint (car tryint)) - (t (list '(%integrate simp) - f - x - a - '$inf)))) - (and (not ($unknown f)) - (errset ($defint f x a '$inf))))) + (let ((tryint (and (not ($unknown f)) (errset ($defint f x a '$inf))))) + (if tryint + (car tryint) + (list '(%integrate simp) f x a '$inf)))) ;;;CREATES UNIQUE NAMES FOR VARIABLE OF INTEGRATION (defun createname (head tail) @@ -313,7 +309,8 @@ ;;;REDUCES LAPLACE(F(T)/T**N,T,S) CASE TO LAPLACE(F(T)/T**(N-1),T,S) CASE (defun hackit (exponent rest) (cond ((equal exponent -1) - ((lambda (parm) (laptimes rest)) (createname parm 1))) + (let ((parm (createname parm 1))) + (laptimes rest))) (t (mydefint (hackit (1+ exponent) rest) (createname parm (- -1 exponent)) (createname parm (- exponent)))))) @@ -334,80 +331,72 @@ '(laplace)) (cdr fun)))))))) +;;;COMPUTES %E**(W*B*%I)*F(S-W*A*%I) WHERE W=-1 IF SIGN IS T ELSE W=1 (defun mostpart (f parm sign a b) - ;;;COMPUTES %E**(W*B*%I)*F(S-W*A*%I) WHERE W=-1 IF SIGN IS T ELSE W=1 - ((lambda (substinfun) - (cond ((zerop1 b) substinfun) - (t (list '(mtimes) - (exponentiate (afixsign (list '(mtimes) - b - '$%i) - (null sign))) - substinfun)))) - ($at f - (list '(mequal simp) - parm - (list '(mplus simp) - parm - (afixsign (list '(mtimes) a '$%i) sign)))))) + (let ((substinfun ($at f + (list '(mequal simp) + parm + (list '(mplus simp) parm (afixsign (list '(mtimes) a '$%i) sign)))))) + (if (zerop1 b) + substinfun + (list '(mtimes) + (exponentiate (afixsign (list '(mtimes) b '$%i) (null sign))) + substinfun)))) ;;;IF WHICHSIGN IS NIL THEN SIN TRANSFORM ELSE COS TRANSFORM (defun compose (fun parm whichsign a b) - ((lambda (result) - ($ratsimp (simptimes (cons '(mtimes) - (cond (whichsign result) - (t (cons '$%i result)))) - 1 nil))) - (list '((rat) 1 2) - (list '(mplus) - (mostpart fun parm t a b) - (afixsign (mostpart fun parm nil a b) - whichsign))))) + (let ((result (list '((rat) 1 2) + (list '(mplus) + (mostpart fun parm t a b) + (afixsign (mostpart fun parm nil a b) + whichsign))))) + ($ratsimp (simptimes (cons '(mtimes) + (if whichsign + result + (cons '$%i result))) + 1 nil)))) ;;;FUN IS OF THE FORM SIN(A*T+B)*REST(T) OR COS (defun lapsin (fun rest trigswitch) - ((lambda (ab) - (cond - (ab - (cond - (rest (compose (laptimes rest) - parm - trigswitch - (car ab) - (cdr ab))) - (t (simptimes - (list - '(mtimes) - (cond - ((zerop1 (cdr ab)) - (cond (trigswitch parm) (t (car ab)))) - (t (cond (trigswitch (list '(mplus) - (list '(mtimes) - parm - (list '(%cos) - (cdr ab))) - (list '(mtimes) - -1 - (car ab) - (list '(%sin) - (cdr ab))))) - (t (list '(mplus) - (list '(mtimes) - parm - (list '(%sin) - (cdr ab))) - (list '(mtimes) - (car ab) - (list '(%cos) - (cdr ab)))))))) - (list '(mexpt) - (list '(mplus) - (list '(mexpt) parm 2) - (list '(mexpt) (car ab) 2)) - -1)) - 1 nil)))) - (t (lapshift fun rest)))) - (islinear (cadr fun) var))) + (let ((ab (islinear (cadr fun) var))) + (cond (ab + (cond (rest + (compose (laptimes rest) + parm + trigswitch + (car ab) + (cdr ab))) + (t + (simptimes + (list '(mtimes) + (cond ((zerop1 (cdr ab)) + (if trigswitch parm (car ab))) + (t + (cond (trigswitch + (list '(mplus) + (list '(mtimes) + parm + (list '(%cos) (cdr ab))) + (list '(mtimes) + -1 + (car ab) + (list '(%sin) (cdr ab))))) + (t + (list '(mplus) + (list '(mtimes) + parm + (list '(%sin) (cdr ab))) + (list '(mtimes) + (car ab) + (list '(%cos) (cdr ab)))))))) + (list '(mexpt) + (list '(mplus) + (list '(mexpt) parm 2) + (list '(mexpt) (car ab) 2)) + -1)) + 1 nil)))) + (t + (lapshift fun rest))))) ;;;FUN IS OF THE FORM SINH(A*T+B)*REST(T) OR IS COSH (defun lapsinh (fun rest switch) @@ -436,74 +425,64 @@ ;;;FUN IS OF THE FORM LOG(A*T) (defun laplog (fun) - ((lambda (ab) - (cond ((and ab (zerop1 (cdr ab))) - (simptimes (list '(mtimes) - (list '(mplus) - (subfunmake '$psi '(0) (ncons 1)) - (list '(%log) (car ab)) - (list '(mtimes) - -1 - (list '(%log) parm))) - (list '(mexpt) parm -1)) - 1 nil)) - (t (lapdefint fun)))) - (islinear (cadr fun) var))) + (let ((ab (islinear (cadr fun) var))) + (cond ((and ab (zerop1 (cdr ab))) + (simptimes (list '(mtimes) + (list '(mplus) + (subfunmake '$psi '(0) (ncons 1)) + (list '(%log) (car ab)) + (list '(mtimes) -1 (list '(%log) parm))) + (list '(mexpt) parm -1)) + 1 nil)) + (t + (lapdefint fun))))) (defun raiseup (fbase exponent) - (cond ((equal exponent 1) fbase) - (t (list '(mexpt) fbase exponent)))) + (if (equal exponent 1) + fbase + (list '(mexpt) fbase exponent))) +;;TAKES TRANSFORM OF DELTA(A*T+B)*F(T) (defun lapdelta (fun rest) - ;;TAKES TRANSFORM OF DELTA(A*T+B)*F(T) - ((lambda (ab sign recipa) - (cond - (ab - (setq recipa (power (car ab) -1) ab (div (cdr ab) (car ab))) - (setq sign (asksign ab) recipa (simplifya (list '(mabs) recipa) nil)) - (simplifya (cond ((eq sign '$positive) 0) - ((eq sign '$zero) - (list '(mtimes) - (maxima-substitute 0 var (fixuprest rest)) - recipa)) - (t (list '(mtimes) - (maxima-substitute (neg ab) - var - (fixuprest rest)) - (list '(mexpt) - '$%e - (cons '(mtimes) - (cons parm (ncons ab)))) - recipa))) - nil)) - (t (lapshift fun rest)))) - (islinear (cadr fun) var) nil nil)) + (let ((ab (islinear (cadr fun) var)) + (sign nil) + (recipa nil)) + (cond (ab + (setq recipa (power (car ab) -1) ab (div (cdr ab) (car ab))) + (setq sign (asksign ab) recipa (simplifya (list '(mabs) recipa) nil)) + (simplifya (cond ((eq sign '$positive) + 0) + ((eq sign '$zero) + (list '(mtimes) + (maxima-substitute 0 var (fixuprest rest)) + recipa)) + (t + (list '(mtimes) + (maxima-substitute (neg ab) var (fixuprest rest)) + (list '(mexpt) '$%e (cons '(mtimes) (cons parm (ncons ab)))) + recipa))) + nil)) + (t + (lapshift fun rest))))) (defun laperf (fun) - ((lambda (ab) - (cond ((and ab (equal (cdr ab) 0)) - (simptimes (list '(mtimes) - (div* (exponentiate (div* (list '(mexpt) - parm - 2) - (list '(mtimes) - 4 - (list '(mexpt) - (car ab) - 2)))) - parm) - (list '(mplus) - 1 - (list '(mtimes) - -1 - (list '(%erf) - (div* parm - (list '(mtimes) - 2 - (car ab)))) - ))) 1 nil)) - (t (lapdefint fun)))) - (islinear (cadr fun) var))) + (let ((ab (islinear (cadr fun) var))) + (cond ((and ab (equal (cdr ab) 0)) + (simptimes (list '(mtimes) + (div* (exponentiate (div* (list '(mexpt) parm 2) + (list '(mtimes) + 4 + (list '(mexpt) (car ab) 2)))) + parm) + (list '(mplus) + 1 + (list '(mtimes) + -1 + (list '(%erf) (div* parm (list '(mtimes) 2 (car ab))))))) + 1 + nil)) + (t + (lapdefint fun))))) (defun lapdefint (fun) (prog (tryint mult) Index: rpart.lisp =================================================================== RCS file: /cvsroot/maxima/maxima/src/rpart.lisp,v retrieving revision 1.16 retrieving revision 1.17 diff -u -d -r1.16 -r1.17 --- rpart.lisp 24 Feb 2009 21:55:51 -0000 1.16 +++ rpart.lisp 28 Feb 2009 15:44:01 -0000 1.17 @@ -24,10 +24,10 @@ $logarc rischp $keepfloat complexsign)) (defmvar implicit-real nil "If t RPART assumes radicals and logs - of real quantities are real and doesn't ask sign questions") + of real quantities are real and doesn't ask sign questions") ... [truncated message content] |