1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;; 2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3;;; The data in this file contains enhancments. ;;;;; 4;;; ;;;;; 5;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;; 6;;; All rights reserved ;;;;; 7;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 8;;; (c) Copyright 1981 Massachusetts Institute of Technology ;;; 9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13;;; Interpolation routine by CFFK. 14;;; Bigfloat version by rtoy. 15 16(macsyma-module intpol) 17 18(load-macsyma-macros transm) 19 20(defmvar $find_root_abs 0.0 21 "Desired absolute error in the root found by find_root") 22(defmvar $find_root_rel 0.0 23 "Desired relative error in the root found by find_root") 24(defmvar $find_root_error t 25 "If true, find_root and bf_find_root prints an error message. 26 Otherwise the value of find_root_error is returned.") 27 28(defmspec $interpolate (form) 29 (format t "NOTE: The interpolate function has been renamed to find_root. 30The variables intpolabs, intpolrel, and intpolerror have been renamed 31to find_root_abs, find_root_rel, and find_root_error, respectively. 32Perhaps you meant to enter `~a'.~%" 33 (print-invert-case (implode (mstring `(($find_root) ,@(cdr form)))))) 34 '$done) 35 36(in-package :bigfloat) 37 38;; Define FIND-ROOT-SUBR and INTERPOLATE-CHECK in the BIGFLOAT package 39;; so we don't have to write BIGFLOAT::foo for all of the arithmetic 40;; operations. 41 42(defun find-root-subr (f left right 43 &key (abserr maxima::$find_root_abs) 44 (relerr maxima::$find_root_rel)) 45 (flet ((convert (s) 46 ;; Try to convert to a BIGFLOAT type. If that fails, just 47 ;; return the argument. Set the flags errcatch and erromsg 48 ;; so we can catch the error, but disable printing of any 49 ;; error messages. 50 (let ((maxima::errcatch t) 51 (maxima::$errormsg nil)) 52 (or (car (maxima::errset (to s))) 53 s)))) 54 (let (;; Don't want to bind $numer to T here. This causes things 55 ;; like log(400)^400 to be computed using double-floats 56 ;; (which overflows), which is not what we want if we're 57 ;; doing bfloat arithmetic. Could bind it for 58 ;; double-floats, but all find_root tests pass without this. 59 #+(or) 60 (maxima::$numer t) 61 (maxima::$%enumer t)) 62 (setq left (convert left) 63 right (convert right))) 64 (unless (and (numberp left) (numberp right)) 65 ;; The interval boundaries must have numerical values 66 (return-from find-root-subr (values nil left right))) 67 (when (< right left) 68 ;; Make left the lower and right the upper bound of the interval 69 (psetq left right right left)) 70 (let ((lin 0) (a left) (b right) 71 (fa (convert (funcall f (maxima::to left)))) 72 (fb (convert (funcall f (maxima::to right)))) c fc) 73 (unless (and (numberp fa) (numberp fb)) 74 (return-from find-root-subr (values nil a b))) 75 (when (<= (abs fa) (to abserr)) 76 ;; If a or b is already small enough, return it as the root 77 (return-from find-root-subr a)) 78 (when (<= (abs fb) (to abserr)) 79 (return-from find-root-subr b)) 80 (when (plusp (* (float-sign fa) (float-sign fb))) 81 (if (eq maxima::$find_root_error t) 82 (maxima::merror (intl:gettext "find_root: function has same sign at endpoints: ~M, ~M") 83 `((maxima::mequal) ((maxima::f) ,a) ,fa) 84 `((maxima::mequal) ((maxima::f) ,b) ,fb)) 85 (return-from find-root-subr 'maxima::$find_root_error))) 86 (when (plusp fa) 87 (psetq fa fb 88 fb fa 89 a b 90 b a)) 91 ;; Use binary search to close in on the root 92 (loop while (< lin 3) do 93 (setq c (* 0.5 (+ a b)) 94 fc (convert (funcall f (maxima::to c)))) 95 (unless (numberp fc) 96 (return-from find-root-subr (values nil a b))) 97 (when (interpolate-check a c b fc abserr relerr) 98 (return-from find-root-subr c)) 99 (if (< (abs (- fc (* 0.5 (+ fa fb)))) (* 0.1 (- fb fa))) 100 (incf lin) 101 (setq lin 0)) 102 (if (plusp fc) 103 (setq fb fc b c) 104 (setq fa fc a c))) 105 ;; Now use the regula falsi 106 (loop 107 (setq c (if (plusp (+ fb fa)) 108 (+ a (* (- b a) (/ fa (- fa fb)))) 109 (+ b (* (- a b) (/ fb (- fb fa))))) 110 fc (convert (funcall f (maxima::to c)))) 111 (unless (numberp fc) 112 (return-from find-root-subr (values nil a b))) 113 (when (interpolate-check a c b fc abserr relerr) 114 (return-from find-root-subr c)) 115 (if (plusp fc) 116 (setq fb fc b c) 117 (setq fa fc a c)))))) 118 119(defun interpolate-check (a c b fc abserr relerr) 120 (not (and (prog1 121 (> (abs fc) (to abserr)) 122 (setq fc (max (abs a) (abs b)))) 123 (> (abs (- b c)) (* (to relerr) fc)) 124 (> (abs (- c a)) (* (to relerr) fc))))) 125 126(in-package :maxima) 127(defun %find-root (name fun-or-expr args) 128 ;; Extract the keyword arguments from args, if any. 129 (let (non-keyword keywords) 130 (loop for arg in args 131 do (if (and (listp arg) 132 (eq (caar arg) 'mequal)) 133 (push arg keywords) 134 (push arg non-keyword))) 135 (setf non-keyword (nreverse non-keyword)) 136 (setf keywords (nreverse keywords)) 137 (when keywords 138 (setf keywords (lispify-maxima-keyword-options keywords '($abserr $relerr)))) 139 #+(or) 140 (progn 141 (format t "keyword args = ~S~%" keywords) 142 (format t "non-keyword args = ~S~%" non-keyword)) 143 (multiple-value-bind (coerce-float fl) 144 ;; The name tells us what error values to use, how to coerce the 145 ;; function, and what function to use to convert to the desired 146 ;; float type. 147 (ecase name 148 ($find_root 149 (values 'coerce-float-fun '$float)) 150 ($bf_find_root 151 (values 'coerce-bfloat-fun '$bfloat))) 152 (case (length non-keyword) 153 (2 154 ;; function case: f, lo, hi 155 (multiple-value-bind (result left right) 156 (apply 'bigfloat::find-root-subr (funcall coerce-float fun-or-expr) 157 (funcall fl (first non-keyword)) 158 (funcall fl (second non-keyword)) 159 keywords) 160 (if (bigfloat:numberp result) 161 (to result) 162 (if (eq result '$find_root_error) 163 $find_root_error 164 `((,name) ,fun-or-expr ,(to left) ,(to right)))))) 165 (3 166 ;; expr case: expr, var, lo, hi 167 (multiple-value-bind (result left right) 168 (apply 'bigfloat::find-root-subr 169 (funcall coerce-float (sub ($lhs fun-or-expr) ($rhs fun-or-expr)) 170 `((mlist) ,(first non-keyword))) 171 (funcall fl (second non-keyword)) 172 (funcall fl (third non-keyword)) 173 keywords) 174 (if (bigfloat:numberp result) 175 (to result) 176 (if (eq result '$find_root_error) 177 $find_root_error 178 `((,name) ,fun-or-expr ,(first non-keyword) ,(to left) ,(to right)))))) 179 (t 180 ;; wrong number of args 181 (wna-err name)))))) 182 183(defmfun $find_root (fun-or-expr &rest args) 184 (%find-root '$find_root fun-or-expr args)) 185 186;; Like find_root but operates on bfloats and returns a bfloat result. 187(defmfun $bf_find_root (fun-or-expr &rest args) 188 (%find-root '$bf_find_root fun-or-expr args)) 189