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