src/intpol.lisp

   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.
  30 The variables intpolabs, intpolrel, and intpolerror have been renamed
  31 to find_root_abs, find_root_rel, and find_root_error, respectively.
  32 Perhaps 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 errset, errcatch,
  48            ;; $errormsg and *mdebug* so we can catch the error, but
  49            ;; disable printing of any error messages.
  50            (let ((maxima::errset nil)
  51                  (maxima::errcatch t)
  52                  (maxima::$errormsg nil)
  53                  (maxima::*mdebug* nil))
  54              (or (car (maxima::errset (to s)))
  55                  s))))
  56     (let (;; Don't want to bind $numer to T here.  This causes things
  57           ;; like log(400)^400 to be computed using double-floats
  58           ;; (which overflows), which is not what we want if we're
  59           ;; doing bfloat arithmetic.  Could bind it for
  60           ;; double-floats, but all find_root tests pass without this.
  61           #+(or)
  62           (maxima::$numer t)
  63           (maxima::$%enumer t))
  64       (setq left (convert left)
  65             right (convert right)))
  66     (unless (and (numberp left) (numberp right))
  67       ;; The interval boundaries must have numerical values
  68       (return-from find-root-subr (values nil left right)))
  69     (when (< right left)
  70       ;; Make left the lower and right the upper bound of the interval
  71       (psetq left right right left))
  72     (let ((lin 0) (a left) (b right)
  73           (fa (convert (funcall f (maxima::to left))))
  74           (fb (convert (funcall f (maxima::to right)))) c fc)
  75       (unless (and (numberp fa) (numberp fb))
  76         (return-from find-root-subr (values nil a b)))
  77       (when (<= (abs fa) (to abserr))
  78         ;; If a or b is already small enough, return it as the root
  79         (return-from find-root-subr a))
  80       (when (<= (abs fb) (to abserr))
  81         (return-from find-root-subr b))
  82       (when (plusp (* (float-sign fa) (float-sign fb)))
  83         (if (eq maxima::$find_root_error t)
  84             (maxima::merror (intl:gettext "find_root: function has same sign at endpoints: ~M, ~M")
  85                             `((maxima::mequal) ((maxima::f) ,a) ,fa)
  86                             `((maxima::mequal) ((maxima::f) ,b) ,fb))
  87             (return-from find-root-subr 'maxima::$find_root_error)))
  88       (when (plusp fa)
  89         (psetq fa fb
  90                fb fa
  91                a b
  92                b a))
  93       ;; Use binary search to close in on the root
  94       (loop while (< lin 3) do
  95            (setq c (* 0.5 (+ a b))
  96                  fc (convert (funcall f (maxima::to c))))
  97            (unless (numberp fc)
  98              (return-from find-root-subr (values nil a b)))
  99            (when (interpolate-check a c b fc abserr relerr)
 100              (return-from find-root-subr c))
 101            (if (< (abs (- fc (* 0.5 (+ fa fb)))) (* 0.1 (- fb fa)))
 102                (incf lin)
 103                (setq lin 0))
 104            (if (plusp fc)
 105                (setq fb fc b c)
 106                (setq fa fc a c)))
 107       ;; Now use the regula falsi
 108       (loop
 109          (setq c (if (plusp (+ fb fa))
 110                      (+ a (* (- b a) (/ fa (- fa fb))))
 111                      (+ b (* (- a b) (/ fb (- fb fa)))))
 112                fc (convert (funcall f (maxima::to c))))
 113          (unless (numberp fc)
 114            (return-from find-root-subr (values nil a b)))
 115          (when (interpolate-check a c b fc abserr relerr)
 116            (return-from find-root-subr c))
 117          (if (plusp fc)
 118              (setq fb fc b c)
 119              (setq fa fc a c))))))
 120
 121 (defun interpolate-check (a c b fc abserr relerr)
 122   (not (and (prog1
 123                 (> (abs fc) (to abserr))
 124               (setq fc (max (abs a) (abs b))))
 125             (> (abs (- b c)) (* (to relerr) fc))
 126             (> (abs (- c a)) (* (to relerr) fc)))))
 127
 128 (in-package :maxima)
 129 (defun %find-root (name fun-or-expr args)
 130   ;; Extract the keyword arguments from args, if any.
 131   (let (non-keyword keywords)
 132     (loop for arg in args
 133        do (if (and (listp arg)
 134                    (eq (caar arg) 'mequal))
 135               (push arg keywords)
 136               (push arg non-keyword)))
 137     (setf non-keyword (nreverse non-keyword))
 138     (setf keywords (nreverse keywords))
 139     (when keywords
 140       (setf keywords
 141             (defmfun-keywords name keywords '($abserr $relerr))))
 142     #+(or)
 143     (progn
 144       (format t "keyword args = ~S~%" keywords)
 145       (format t "non-keyword args = ~S~%" non-keyword))
 146     (multiple-value-bind (coerce-float fl)
 147         ;; The name tells us what error values to use, how to coerce the
 148         ;; function, and what function to use to convert to the desired
 149         ;; float type.
 150         (ecase name
 151           ($find_root
 152            (values 'coerce-float-fun '$float))
 153           ($bf_find_root
 154            (values 'coerce-bfloat-fun '$bfloat)))
 155       (case (length non-keyword)
 156         (2
 157          ;; function case: f, lo, hi
 158          (multiple-value-bind (result left right)
 159              (apply 'bigfloat::find-root-subr (funcall coerce-float fun-or-expr)
 160                     (funcall fl (first non-keyword))
 161                     (funcall fl (second non-keyword))
 162                     keywords)
 163            (if (bigfloat:numberp result)
 164                (to result)
 165                (if (eq result '$find_root_error)
 166                    $find_root_error
 167                    `((,name) ,fun-or-expr ,(to left) ,(to right))))))
 168         (3
 169          ;; expr case: expr, var, lo, hi
 170          (multiple-value-bind (result left right)
 171              (apply 'bigfloat::find-root-subr
 172                     (funcall coerce-float (sub ($lhs fun-or-expr) ($rhs fun-or-expr))
 173                              `((mlist) ,(first non-keyword)))
 174                     (funcall fl (second non-keyword))
 175                     (funcall fl (third non-keyword))
 176                     keywords)
 177            (if (bigfloat:numberp result)
 178                (to result)
 179                (if (eq result '$find_root_error)
 180                    $find_root_error
 181                    `((,name) ,fun-or-expr ,(first non-keyword) ,(to left) ,(to right))))))
 182         (t
 183          ;; wrong number of args
 184          (wna-err name))))))
 185
 186 (defmfun $find_root (fun-or-expr &rest args)
 187   (%find-root '$find_root fun-or-expr args))
 188
 189 ;; Like find_root but operates on bfloats and returns a bfloat result.
 190 (defmfun $bf_find_root (fun-or-expr &rest args)
 191   (%find-root '$bf_find_root fun-or-expr args))