Mercurial > emacs
view lisp/calc/calcalg2.el @ 88018:4e0d27f39575
*** empty log message ***
author | Thien-Thi Nguyen <ttn@gnuvola.org> |
---|---|
date | Sun, 27 Jan 2008 11:06:45 +0000 |
parents | 107ccd98fa12 |
children | 606f2d163a64 1e3a407766b9 |
line wrap: on
line source
;;; calcalg2.el --- more algebraic functions for Calc ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004, ;; 2005, 2006, 2007, 2008 Free Software Foundation, Inc. ;; Author: David Gillespie <daveg@synaptics.com> ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com> ;; This file is part of GNU Emacs. ;; GNU Emacs is free software; you can redistribute it and/or modify ;; it under the terms of the GNU General Public License as published by ;; the Free Software Foundation; either version 3, or (at your option) ;; any later version. ;; GNU Emacs is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;; GNU General Public License for more details. ;; You should have received a copy of the GNU General Public License ;; along with GNU Emacs; see the file COPYING. If not, write to the ;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, ;; Boston, MA 02110-1301, USA. ;;; Commentary: ;;; Code: ;; This file is autoloaded from calc-ext.el. (require 'calc-ext) (require 'calc-macs) (defun calc-derivative (var num) (interactive "sDifferentiate with respect to: \np") (calc-slow-wrapper (when (< num 0) (error "Order of derivative must be positive")) (let ((func (if (calc-is-hyperbolic) 'calcFunc-tderiv 'calcFunc-deriv)) n expr) (if (or (equal var "") (equal var "$")) (setq n 2 expr (calc-top-n 2) var (calc-top-n 1)) (setq var (math-read-expr var)) (when (eq (car-safe var) 'error) (error "Bad format in expression: %s" (nth 1 var))) (setq n 1 expr (calc-top-n 1))) (while (>= (setq num (1- num)) 0) (setq expr (list func expr var))) (calc-enter-result n "derv" expr)))) (defun calc-integral (var &optional arg) (interactive "sIntegration variable: \nP") (if arg (calc-tabular-command 'calcFunc-integ "Integration" "intg" nil var nil nil) (calc-slow-wrapper (if (or (equal var "") (equal var "$")) (calc-enter-result 2 "intg" (list 'calcFunc-integ (calc-top-n 2) (calc-top-n 1))) (let ((var (math-read-expr var))) (if (eq (car-safe var) 'error) (error "Bad format in expression: %s" (nth 1 var))) (calc-enter-result 1 "intg" (list 'calcFunc-integ (calc-top-n 1) var))))))) (defun calc-num-integral (&optional varname lowname highname) (interactive "sIntegration variable: ") (calc-tabular-command 'calcFunc-ninteg "Integration" "nint" nil varname lowname highname)) (defun calc-summation (arg &optional varname lowname highname) (interactive "P\nsSummation variable: ") (calc-tabular-command 'calcFunc-sum "Summation" "sum" arg varname lowname highname)) (defun calc-alt-summation (arg &optional varname lowname highname) (interactive "P\nsSummation variable: ") (calc-tabular-command 'calcFunc-asum "Summation" "asum" arg varname lowname highname)) (defun calc-product (arg &optional varname lowname highname) (interactive "P\nsIndex variable: ") (calc-tabular-command 'calcFunc-prod "Index" "prod" arg varname lowname highname)) (defun calc-tabulate (arg &optional varname lowname highname) (interactive "P\nsIndex variable: ") (calc-tabular-command 'calcFunc-table "Index" "tabl" arg varname lowname highname)) (defun calc-tabular-command (func prompt prefix arg varname lowname highname) (calc-slow-wrapper (let (var (low nil) (high nil) (step nil) stepname stepnum (num 1) expr) (if (consp arg) (setq stepnum 1) (setq stepnum 0)) (if (or (equal varname "") (equal varname "$") (null varname)) (setq high (calc-top-n (+ stepnum 1)) low (calc-top-n (+ stepnum 2)) var (calc-top-n (+ stepnum 3)) num (+ stepnum 4)) (setq var (if (stringp varname) (math-read-expr varname) varname)) (if (eq (car-safe var) 'error) (error "Bad format in expression: %s" (nth 1 var))) (or lowname (setq lowname (read-string (concat prompt " variable: " varname ", from: ")))) (if (or (equal lowname "") (equal lowname "$")) (setq high (calc-top-n (+ stepnum 1)) low (calc-top-n (+ stepnum 2)) num (+ stepnum 3)) (setq low (if (stringp lowname) (math-read-expr lowname) lowname)) (if (eq (car-safe low) 'error) (error "Bad format in expression: %s" (nth 1 low))) (or highname (setq highname (read-string (concat prompt " variable: " varname ", from: " lowname ", to: ")))) (if (or (equal highname "") (equal highname "$")) (setq high (calc-top-n (+ stepnum 1)) num (+ stepnum 2)) (setq high (if (stringp highname) (math-read-expr highname) highname)) (if (eq (car-safe high) 'error) (error "Bad format in expression: %s" (nth 1 high))) (if (consp arg) (progn (setq stepname (read-string (concat prompt " variable: " varname ", from: " lowname ", to: " highname ", step: "))) (if (or (equal stepname "") (equal stepname "$")) (setq step (calc-top-n 1) num 2) (setq step (math-read-expr stepname)) (if (eq (car-safe step) 'error) (error "Bad format in expression: %s" (nth 1 step))))))))) (or step (if (consp arg) (setq step (calc-top-n 1)) (if arg (setq step (prefix-numeric-value arg))))) (setq expr (calc-top-n num)) (calc-enter-result num prefix (append (list func expr var low high) (and step (list step))))))) (defun calc-solve-for (var) (interactive "sVariable(s) to solve for: ") (calc-slow-wrapper (let ((func (if (calc-is-inverse) (if (calc-is-hyperbolic) 'calcFunc-ffinv 'calcFunc-finv) (if (calc-is-hyperbolic) 'calcFunc-fsolve 'calcFunc-solve)))) (if (or (equal var "") (equal var "$")) (calc-enter-result 2 "solv" (list func (calc-top-n 2) (calc-top-n 1))) (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var) (not (string-match "\\[" var))) (math-read-expr (concat "[" var "]")) (math-read-expr var)))) (if (eq (car-safe var) 'error) (error "Bad format in expression: %s" (nth 1 var))) (calc-enter-result 1 "solv" (list func (calc-top-n 1) var))))))) (defun calc-poly-roots (var) (interactive "sVariable to solve for: ") (calc-slow-wrapper (if (or (equal var "") (equal var "$")) (calc-enter-result 2 "prts" (list 'calcFunc-roots (calc-top-n 2) (calc-top-n 1))) (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var) (not (string-match "\\[" var))) (math-read-expr (concat "[" var "]")) (math-read-expr var)))) (if (eq (car-safe var) 'error) (error "Bad format in expression: %s" (nth 1 var))) (calc-enter-result 1 "prts" (list 'calcFunc-roots (calc-top-n 1) var)))))) (defun calc-taylor (var nterms) (interactive "sTaylor expansion variable: \nNNumber of terms: ") (calc-slow-wrapper (let ((var (math-read-expr var))) (if (eq (car-safe var) 'error) (error "Bad format in expression: %s" (nth 1 var))) (calc-enter-result 1 "tylr" (list 'calcFunc-taylor (calc-top-n 1) var (prefix-numeric-value nterms)))))) ;; The following are global variables used by math-derivative and some ;; related functions (defvar math-deriv-var) (defvar math-deriv-total) (defvar math-deriv-symb) (defvar math-decls-cache) (defvar math-decls-all) (defun math-derivative (expr) (cond ((equal expr math-deriv-var) 1) ((or (Math-scalarp expr) (eq (car expr) 'sdev) (and (eq (car expr) 'var) (or (not math-deriv-total) (math-const-var expr) (progn (math-setup-declarations) (memq 'const (nth 1 (or (assq (nth 2 expr) math-decls-cache) math-decls-all))))))) 0) ((eq (car expr) '+) (math-add (math-derivative (nth 1 expr)) (math-derivative (nth 2 expr)))) ((eq (car expr) '-) (math-sub (math-derivative (nth 1 expr)) (math-derivative (nth 2 expr)))) ((memq (car expr) '(calcFunc-eq calcFunc-neq calcFunc-lt calcFunc-gt calcFunc-leq calcFunc-geq)) (list (car expr) (math-derivative (nth 1 expr)) (math-derivative (nth 2 expr)))) ((eq (car expr) 'neg) (math-neg (math-derivative (nth 1 expr)))) ((eq (car expr) '*) (math-add (math-mul (nth 2 expr) (math-derivative (nth 1 expr))) (math-mul (nth 1 expr) (math-derivative (nth 2 expr))))) ((eq (car expr) '/) (math-sub (math-div (math-derivative (nth 1 expr)) (nth 2 expr)) (math-div (math-mul (nth 1 expr) (math-derivative (nth 2 expr))) (math-sqr (nth 2 expr))))) ((eq (car expr) '^) (let ((du (math-derivative (nth 1 expr))) (dv (math-derivative (nth 2 expr)))) (or (Math-zerop du) (setq du (math-mul (nth 2 expr) (math-mul (math-normalize (list '^ (nth 1 expr) (math-add (nth 2 expr) -1))) du)))) (or (Math-zerop dv) (setq dv (math-mul (math-normalize (list 'calcFunc-ln (nth 1 expr))) (math-mul expr dv)))) (math-add du dv))) ((eq (car expr) '%) (math-derivative (nth 1 expr))) ; a reasonable definition ((eq (car expr) 'vec) (math-map-vec 'math-derivative expr)) ((and (memq (car expr) '(calcFunc-conj calcFunc-re calcFunc-im)) (= (length expr) 2)) (list (car expr) (math-derivative (nth 1 expr)))) ((and (memq (car expr) '(calcFunc-subscr calcFunc-mrow calcFunc-mcol)) (= (length expr) 3)) (let ((d (math-derivative (nth 1 expr)))) (if (math-numberp d) 0 ; assume x and x_1 are independent vars (list (car expr) d (nth 2 expr))))) (t (or (and (symbolp (car expr)) (if (= (length expr) 2) (let ((handler (get (car expr) 'math-derivative))) (and handler (let ((deriv (math-derivative (nth 1 expr)))) (if (Math-zerop deriv) deriv (math-mul (funcall handler (nth 1 expr)) deriv))))) (let ((handler (get (car expr) 'math-derivative-n))) (and handler (funcall handler expr))))) (and (not (eq math-deriv-symb 'pre-expand)) (let ((exp (math-expand-formula expr))) (and exp (or (let ((math-deriv-symb 'pre-expand)) (catch 'math-deriv (math-derivative expr))) (math-derivative exp))))) (if (or (Math-objvecp expr) (eq (car expr) 'var) (not (symbolp (car expr)))) (if math-deriv-symb (throw 'math-deriv nil) (list (if math-deriv-total 'calcFunc-tderiv 'calcFunc-deriv) expr math-deriv-var)) (let ((accum 0) (arg expr) (n 1) derv) (while (setq arg (cdr arg)) (or (Math-zerop (setq derv (math-derivative (car arg)))) (let ((func (intern (concat (symbol-name (car expr)) "'" (if (> n 1) (int-to-string n) "")))) (prop (cond ((= (length expr) 2) 'math-derivative-1) ((= (length expr) 3) 'math-derivative-2) ((= (length expr) 4) 'math-derivative-3) ((= (length expr) 5) 'math-derivative-4) ((= (length expr) 6) 'math-derivative-5)))) (setq accum (math-add accum (math-mul derv (let ((handler (get func prop))) (or (and prop handler (apply handler (cdr expr))) (if (and math-deriv-symb (not (get func 'calc-user-defn))) (throw 'math-deriv nil) (cons func (cdr expr)))))))))) (setq n (1+ n))) accum)))))) (defun calcFunc-deriv (expr math-deriv-var &optional deriv-value math-deriv-symb) (let* ((math-deriv-total nil) (res (catch 'math-deriv (math-derivative expr)))) (or (eq (car-safe res) 'calcFunc-deriv) (null res) (setq res (math-normalize res))) (and res (if deriv-value (math-expr-subst res math-deriv-var deriv-value) res)))) (defun calcFunc-tderiv (expr math-deriv-var &optional deriv-value math-deriv-symb) (math-setup-declarations) (let* ((math-deriv-total t) (res (catch 'math-deriv (math-derivative expr)))) (or (eq (car-safe res) 'calcFunc-tderiv) (null res) (setq res (math-normalize res))) (and res (if deriv-value (math-expr-subst res math-deriv-var deriv-value) res)))) (put 'calcFunc-inv\' 'math-derivative-1 (function (lambda (u) (math-neg (math-div 1 (math-sqr u)))))) (put 'calcFunc-sqrt\' 'math-derivative-1 (function (lambda (u) (math-div 1 (math-mul 2 (list 'calcFunc-sqrt u)))))) (put 'calcFunc-deg\' 'math-derivative-1 (function (lambda (u) (math-div-float '(float 18 1) (math-pi))))) (put 'calcFunc-rad\' 'math-derivative-1 (function (lambda (u) (math-pi-over-180)))) (put 'calcFunc-ln\' 'math-derivative-1 (function (lambda (u) (math-div 1 u)))) (put 'calcFunc-log10\' 'math-derivative-1 (function (lambda (u) (math-div (math-div 1 (math-normalize '(calcFunc-ln 10))) u)))) (put 'calcFunc-lnp1\' 'math-derivative-1 (function (lambda (u) (math-div 1 (math-add u 1))))) (put 'calcFunc-log\' 'math-derivative-2 (function (lambda (x b) (and (not (Math-zerop b)) (let ((lnv (math-normalize (list 'calcFunc-ln b)))) (math-div 1 (math-mul lnv x))))))) (put 'calcFunc-log\'2 'math-derivative-2 (function (lambda (x b) (let ((lnv (list 'calcFunc-ln b))) (math-neg (math-div (list 'calcFunc-log x b) (math-mul lnv b))))))) (put 'calcFunc-exp\' 'math-derivative-1 (function (lambda (u) (math-normalize (list 'calcFunc-exp u))))) (put 'calcFunc-expm1\' 'math-derivative-1 (function (lambda (u) (math-normalize (list 'calcFunc-expm1 u))))) (put 'calcFunc-sin\' 'math-derivative-1 (function (lambda (u) (math-to-radians-2 (math-normalize (list 'calcFunc-cos u)))))) (put 'calcFunc-cos\' 'math-derivative-1 (function (lambda (u) (math-neg (math-to-radians-2 (math-normalize (list 'calcFunc-sin u))))))) (put 'calcFunc-tan\' 'math-derivative-1 (function (lambda (u) (math-to-radians-2 (math-sqr (math-normalize (list 'calcFunc-sec u))))))) (put 'calcFunc-sec\' 'math-derivative-1 (function (lambda (u) (math-to-radians-2 (math-mul (math-normalize (list 'calcFunc-sec u)) (math-normalize (list 'calcFunc-tan u))))))) (put 'calcFunc-csc\' 'math-derivative-1 (function (lambda (u) (math-neg (math-to-radians-2 (math-mul (math-normalize (list 'calcFunc-csc u)) (math-normalize (list 'calcFunc-cot u)))))))) (put 'calcFunc-cot\' 'math-derivative-1 (function (lambda (u) (math-neg (math-to-radians-2 (math-sqr (math-normalize (list 'calcFunc-csc u)))))))) (put 'calcFunc-arcsin\' 'math-derivative-1 (function (lambda (u) (math-from-radians-2 (math-div 1 (math-normalize (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))))) (put 'calcFunc-arccos\' 'math-derivative-1 (function (lambda (u) (math-from-radians-2 (math-div -1 (math-normalize (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))))) (put 'calcFunc-arctan\' 'math-derivative-1 (function (lambda (u) (math-from-radians-2 (math-div 1 (math-add 1 (math-sqr u))))))) (put 'calcFunc-sinh\' 'math-derivative-1 (function (lambda (u) (math-normalize (list 'calcFunc-cosh u))))) (put 'calcFunc-cosh\' 'math-derivative-1 (function (lambda (u) (math-normalize (list 'calcFunc-sinh u))))) (put 'calcFunc-tanh\' 'math-derivative-1 (function (lambda (u) (math-sqr (math-normalize (list 'calcFunc-sech u)))))) (put 'calcFunc-sech\' 'math-derivative-1 (function (lambda (u) (math-neg (math-mul (math-normalize (list 'calcFunc-sech u)) (math-normalize (list 'calcFunc-tanh u))))))) (put 'calcFunc-csch\' 'math-derivative-1 (function (lambda (u) (math-neg (math-mul (math-normalize (list 'calcFunc-csch u)) (math-normalize (list 'calcFunc-coth u))))))) (put 'calcFunc-coth\' 'math-derivative-1 (function (lambda (u) (math-neg (math-sqr (math-normalize (list 'calcFunc-csch u))))))) (put 'calcFunc-arcsinh\' 'math-derivative-1 (function (lambda (u) (math-div 1 (math-normalize (list 'calcFunc-sqrt (math-add (math-sqr u) 1))))))) (put 'calcFunc-arccosh\' 'math-derivative-1 (function (lambda (u) (math-div 1 (math-normalize (list 'calcFunc-sqrt (math-add (math-sqr u) -1))))))) (put 'calcFunc-arctanh\' 'math-derivative-1 (function (lambda (u) (math-div 1 (math-sub 1 (math-sqr u)))))) (put 'calcFunc-bern\'2 'math-derivative-2 (function (lambda (n x) (math-mul n (list 'calcFunc-bern (math-add n -1) x))))) (put 'calcFunc-euler\'2 'math-derivative-2 (function (lambda (n x) (math-mul n (list 'calcFunc-euler (math-add n -1) x))))) (put 'calcFunc-gammag\'2 'math-derivative-2 (function (lambda (a x) (math-deriv-gamma a x 1)))) (put 'calcFunc-gammaG\'2 'math-derivative-2 (function (lambda (a x) (math-deriv-gamma a x -1)))) (put 'calcFunc-gammaP\'2 'math-derivative-2 (function (lambda (a x) (math-deriv-gamma a x (math-div 1 (math-normalize (list 'calcFunc-gamma a))))))) (put 'calcFunc-gammaQ\'2 'math-derivative-2 (function (lambda (a x) (math-deriv-gamma a x (math-div -1 (math-normalize (list 'calcFunc-gamma a))))))) (defun math-deriv-gamma (a x scale) (math-mul scale (math-mul (math-pow x (math-add a -1)) (list 'calcFunc-exp (math-neg x))))) (put 'calcFunc-betaB\' 'math-derivative-3 (function (lambda (x a b) (math-deriv-beta x a b 1)))) (put 'calcFunc-betaI\' 'math-derivative-3 (function (lambda (x a b) (math-deriv-beta x a b (math-div 1 (list 'calcFunc-beta a b)))))) (defun math-deriv-beta (x a b scale) (math-mul (math-mul (math-pow x (math-add a -1)) (math-pow (math-sub 1 x) (math-add b -1))) scale)) (put 'calcFunc-erf\' 'math-derivative-1 (function (lambda (x) (math-div 2 (math-mul (list 'calcFunc-exp (math-sqr x)) (if calc-symbolic-mode '(calcFunc-sqrt (var pi var-pi)) (math-sqrt-pi))))))) (put 'calcFunc-erfc\' 'math-derivative-1 (function (lambda (x) (math-div -2 (math-mul (list 'calcFunc-exp (math-sqr x)) (if calc-symbolic-mode '(calcFunc-sqrt (var pi var-pi)) (math-sqrt-pi))))))) (put 'calcFunc-besJ\'2 'math-derivative-2 (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besJ (math-add v -1) z) (list 'calcFunc-besJ (math-add v 1) z)) 2)))) (put 'calcFunc-besY\'2 'math-derivative-2 (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besY (math-add v -1) z) (list 'calcFunc-besY (math-add v 1) z)) 2)))) (put 'calcFunc-sum 'math-derivative-n (function (lambda (expr) (if (math-expr-contains (cons 'vec (cdr (cdr expr))) math-deriv-var) (throw 'math-deriv nil) (cons 'calcFunc-sum (cons (math-derivative (nth 1 expr)) (cdr (cdr expr)))))))) (put 'calcFunc-prod 'math-derivative-n (function (lambda (expr) (if (math-expr-contains (cons 'vec (cdr (cdr expr))) math-deriv-var) (throw 'math-deriv nil) (math-mul expr (cons 'calcFunc-sum (cons (math-div (math-derivative (nth 1 expr)) (nth 1 expr)) (cdr (cdr expr))))))))) (put 'calcFunc-integ 'math-derivative-n (function (lambda (expr) (if (= (length expr) 3) (if (equal (nth 2 expr) math-deriv-var) (nth 1 expr) (math-normalize (list 'calcFunc-integ (math-derivative (nth 1 expr)) (nth 2 expr)))) (if (= (length expr) 5) (let ((lower (math-expr-subst (nth 1 expr) (nth 2 expr) (nth 3 expr))) (upper (math-expr-subst (nth 1 expr) (nth 2 expr) (nth 4 expr)))) (math-add (math-sub (math-mul upper (math-derivative (nth 4 expr))) (math-mul lower (math-derivative (nth 3 expr)))) (if (equal (nth 2 expr) math-deriv-var) 0 (math-normalize (list 'calcFunc-integ (math-derivative (nth 1 expr)) (nth 2 expr) (nth 3 expr) (nth 4 expr))))))))))) (put 'calcFunc-if 'math-derivative-n (function (lambda (expr) (and (= (length expr) 4) (list 'calcFunc-if (nth 1 expr) (math-derivative (nth 2 expr)) (math-derivative (nth 3 expr))))))) (put 'calcFunc-subscr 'math-derivative-n (function (lambda (expr) (and (= (length expr) 3) (list 'calcFunc-subscr (nth 1 expr) (math-derivative (nth 2 expr))))))) (defvar math-integ-var '(var X ---)) (defvar math-integ-var-2 '(var Y ---)) (defvar math-integ-vars (list 'f math-integ-var math-integ-var-2)) (defvar math-integ-var-list (list math-integ-var)) (defvar math-integ-var-list-list (list math-integ-var-list)) ;; math-integ-depth is a local variable for math-try-integral, but is used ;; by math-integral and math-tracing-integral ;; which are called (directly or indirectly) by math-try-integral. (defvar math-integ-depth) ;; math-integ-level is a local variable for math-try-integral, but is used ;; by math-integral, math-do-integral, math-tracing-integral, ;; math-sub-integration, math-integrate-by-parts and ;; math-integrate-by-substitution, which are called (directly or ;; indirectly) by math-try-integral. (defvar math-integ-level) ;; math-integral-limit is a local variable for calcFunc-integ, but is ;; used by math-tracing-integral, math-sub-integration and ;; math-try-integration. (defvar math-integral-limit) (defmacro math-tracing-integral (&rest parts) (list 'and 'trace-buffer (list 'save-excursion '(set-buffer trace-buffer) '(goto-char (point-max)) (list 'and '(bolp) '(insert (make-string (- math-integral-limit math-integ-level) 32) (format "%2d " math-integ-depth) (make-string math-integ-level 32))) ;;(list 'condition-case 'err (cons 'insert parts) ;; '(error (insert (prin1-to-string err)))) '(sit-for 0)))) ;;; The following wrapper caches results and avoids infinite recursion. ;;; Each cache entry is: ( A B ) Integral of A is B; ;;; ( A N ) Integral of A failed at level N; ;;; ( A busy ) Currently working on integral of A; ;;; ( A parts ) Currently working, integ-by-parts; ;;; ( A parts2 ) Currently working, integ-by-parts; ;;; ( A cancelled ) Ignore this cache entry; ;;; ( A [B] ) Same result as for math-cur-record = B. ;; math-cur-record is a local variable for math-try-integral, but is used ;; by math-integral, math-replace-integral-parts and math-integrate-by-parts ;; which are called (directly or indirectly) by math-try-integral, as well as ;; by calc-dump-integral-cache (defvar math-cur-record) ;; math-enable-subst and math-any-substs are local variables for ;; calcFunc-integ, but are used by math-integral and math-try-integral. (defvar math-enable-subst) (defvar math-any-substs) ;; math-integ-msg is a local variable for math-try-integral, but is ;; used (both locally and non-locally) by math-integral. (defvar math-integ-msg) (defvar math-integral-cache nil) (defvar math-integral-cache-state nil) (defun math-integral (expr &optional simplify same-as-above) (let* ((simp math-cur-record) (math-cur-record (assoc expr math-integral-cache)) (math-integ-depth (1+ math-integ-depth)) (val 'cancelled)) (math-tracing-integral "Integrating " (math-format-value expr 1000) "...\n") (and math-cur-record (progn (math-tracing-integral "Found " (math-format-value (nth 1 math-cur-record) 1000)) (and (consp (nth 1 math-cur-record)) (math-replace-integral-parts math-cur-record)) (math-tracing-integral " => " (math-format-value (nth 1 math-cur-record) 1000) "\n"))) (or (and math-cur-record (not (eq (nth 1 math-cur-record) 'cancelled)) (or (not (integerp (nth 1 math-cur-record))) (>= (nth 1 math-cur-record) math-integ-level))) (and (math-integral-contains-parts expr) (progn (setq val nil) t)) (unwind-protect (progn (let (math-integ-msg) (if (eq calc-display-working-message 'lots) (progn (calc-set-command-flag 'clear-message) (setq math-integ-msg (format "Working... Integrating %s" (math-format-flat-expr expr 0))) (message "%s" math-integ-msg))) (if math-cur-record (setcar (cdr math-cur-record) (if same-as-above (vector simp) 'busy)) (setq math-cur-record (list expr (if same-as-above (vector simp) 'busy)) math-integral-cache (cons math-cur-record math-integral-cache))) (if (eq simplify 'yes) (progn (math-tracing-integral "Simplifying...") (setq simp (math-simplify expr)) (setq val (if (equal simp expr) (progn (math-tracing-integral " no change\n") (math-do-integral expr)) (math-tracing-integral " simplified\n") (math-integral simp 'no t)))) (or (setq val (math-do-integral expr)) (eq simplify 'no) (let ((simp (math-simplify expr))) (or (equal simp expr) (progn (math-tracing-integral "Trying again after " "simplification...\n") (setq val (math-integral simp 'no t)))))))) (if (eq calc-display-working-message 'lots) (message "%s" math-integ-msg))) (setcar (cdr math-cur-record) (or val (if (or math-enable-subst (not math-any-substs)) math-integ-level 'cancelled))))) (setq val math-cur-record) (while (vectorp (nth 1 val)) (setq val (aref (nth 1 val) 0))) (setq val (if (memq (nth 1 val) '(parts parts2)) (progn (setcar (cdr val) 'parts2) (list 'var 'PARTS val)) (and (consp (nth 1 val)) (nth 1 val)))) (math-tracing-integral "Integral of " (math-format-value expr 1000) " is " (math-format-value val 1000) "\n") val)) (defun math-integral-contains-parts (expr) (if (Math-primp expr) (and (eq (car-safe expr) 'var) (eq (nth 1 expr) 'PARTS) (listp (nth 2 expr))) (while (and (setq expr (cdr expr)) (not (math-integral-contains-parts (car expr))))) expr)) (defun math-replace-integral-parts (expr) (or (Math-primp expr) (while (setq expr (cdr expr)) (and (consp (car expr)) (if (eq (car (car expr)) 'var) (and (eq (nth 1 (car expr)) 'PARTS) (consp (nth 2 (car expr))) (if (listp (nth 1 (nth 2 (car expr)))) (progn (setcar expr (nth 1 (nth 2 (car expr)))) (math-replace-integral-parts (cons 'foo expr))) (setcar (cdr math-cur-record) 'cancelled))) (math-replace-integral-parts (car expr))))))) (defvar math-linear-subst-tried t "Non-nil means that a linear substitution has been tried.") ;; The variable math-has-rules is a local variable for math-try-integral, ;; but is used by math-do-integral, which is called (non-directly) by ;; math-try-integral. (defvar math-has-rules) ;; math-old-integ is a local variable for math-do-integral, but is ;; used by math-sub-integration. (defvar math-old-integ) ;; The variables math-t1, math-t2 and math-t3 are local to ;; math-do-integral, math-try-solve-for and math-decompose-poly, but ;; are used by functions they call (directly or indirectly); ;; math-do-integral calls math-do-integral-methods; ;; math-try-solve-for calls math-try-solve-prod, ;; math-solve-find-root-term and math-solve-find-root-in-prod; ;; math-decompose-poly calls math-solve-poly-funny-powers and ;; math-solve-crunch-poly. (defvar math-t1) (defvar math-t2) (defvar math-t3) (defun math-do-integral (expr) (let ((math-linear-subst-tried nil) math-t1 math-t2) (or (cond ((not (math-expr-contains expr math-integ-var)) (math-mul expr math-integ-var)) ((equal expr math-integ-var) (math-div (math-sqr expr) 2)) ((eq (car expr) '+) (and (setq math-t1 (math-integral (nth 1 expr))) (setq math-t2 (math-integral (nth 2 expr))) (math-add math-t1 math-t2))) ((eq (car expr) '-) (and (setq math-t1 (math-integral (nth 1 expr))) (setq math-t2 (math-integral (nth 2 expr))) (math-sub math-t1 math-t2))) ((eq (car expr) 'neg) (and (setq math-t1 (math-integral (nth 1 expr))) (math-neg math-t1))) ((eq (car expr) '*) (cond ((not (math-expr-contains (nth 1 expr) math-integ-var)) (and (setq math-t1 (math-integral (nth 2 expr))) (math-mul (nth 1 expr) math-t1))) ((not (math-expr-contains (nth 2 expr) math-integ-var)) (and (setq math-t1 (math-integral (nth 1 expr))) (math-mul math-t1 (nth 2 expr)))) ((memq (car-safe (nth 1 expr)) '(+ -)) (math-integral (list (car (nth 1 expr)) (math-mul (nth 1 (nth 1 expr)) (nth 2 expr)) (math-mul (nth 2 (nth 1 expr)) (nth 2 expr))) 'yes t)) ((memq (car-safe (nth 2 expr)) '(+ -)) (math-integral (list (car (nth 2 expr)) (math-mul (nth 1 (nth 2 expr)) (nth 1 expr)) (math-mul (nth 2 (nth 2 expr)) (nth 1 expr))) 'yes t)))) ((eq (car expr) '/) (cond ((and (not (math-expr-contains (nth 1 expr) math-integ-var)) (not (math-equal-int (nth 1 expr) 1))) (and (setq math-t1 (math-integral (math-div 1 (nth 2 expr)))) (math-mul (nth 1 expr) math-t1))) ((not (math-expr-contains (nth 2 expr) math-integ-var)) (and (setq math-t1 (math-integral (nth 1 expr))) (math-div math-t1 (nth 2 expr)))) ((and (eq (car-safe (nth 1 expr)) '*) (not (math-expr-contains (nth 1 (nth 1 expr)) math-integ-var))) (and (setq math-t1 (math-integral (math-div (nth 2 (nth 1 expr)) (nth 2 expr)))) (math-mul math-t1 (nth 1 (nth 1 expr))))) ((and (eq (car-safe (nth 1 expr)) '*) (not (math-expr-contains (nth 2 (nth 1 expr)) math-integ-var))) (and (setq math-t1 (math-integral (math-div (nth 1 (nth 1 expr)) (nth 2 expr)))) (math-mul math-t1 (nth 2 (nth 1 expr))))) ((and (eq (car-safe (nth 2 expr)) '*) (not (math-expr-contains (nth 1 (nth 2 expr)) math-integ-var))) (and (setq math-t1 (math-integral (math-div (nth 1 expr) (nth 2 (nth 2 expr))))) (math-div math-t1 (nth 1 (nth 2 expr))))) ((and (eq (car-safe (nth 2 expr)) '*) (not (math-expr-contains (nth 2 (nth 2 expr)) math-integ-var))) (and (setq math-t1 (math-integral (math-div (nth 1 expr) (nth 1 (nth 2 expr))))) (math-div math-t1 (nth 2 (nth 2 expr))))) ((eq (car-safe (nth 2 expr)) 'calcFunc-exp) (math-integral (math-mul (nth 1 expr) (list 'calcFunc-exp (math-neg (nth 1 (nth 2 expr))))))))) ((eq (car expr) '^) (cond ((not (math-expr-contains (nth 1 expr) math-integ-var)) (or (and (setq math-t1 (math-is-polynomial (nth 2 expr) math-integ-var 1)) (math-div expr (math-mul (nth 1 math-t1) (math-normalize (list 'calcFunc-ln (nth 1 expr)))))) (math-integral (list 'calcFunc-exp (math-mul (nth 2 expr) (math-normalize (list 'calcFunc-ln (nth 1 expr))))) 'yes t))) ((not (math-expr-contains (nth 2 expr) math-integ-var)) (if (and (integerp (nth 2 expr)) (< (nth 2 expr) 0)) (math-integral (list '/ 1 (math-pow (nth 1 expr) (- (nth 2 expr)))) nil t) (or (and (setq math-t1 (math-is-polynomial (nth 1 expr) math-integ-var 1)) (setq math-t2 (math-add (nth 2 expr) 1)) (math-div (math-pow (nth 1 expr) math-t2) (math-mul math-t2 (nth 1 math-t1)))) (and (Math-negp (nth 2 expr)) (math-integral (math-div 1 (math-pow (nth 1 expr) (math-neg (nth 2 expr)))) nil t)) nil)))))) ;; Integral of a polynomial. (and (setq math-t1 (math-is-polynomial expr math-integ-var 20)) (let ((accum 0) (n 1)) (while math-t1 (if (setq accum (math-add accum (math-div (math-mul (car math-t1) (math-pow math-integ-var n)) n)) math-t1 (cdr math-t1)) (setq n (1+ n)))) accum)) ;; Try looking it up! (cond ((= (length expr) 2) (and (symbolp (car expr)) (setq math-t1 (get (car expr) 'math-integral)) (progn (while (and math-t1 (not (setq math-t2 (funcall (car math-t1) (nth 1 expr))))) (setq math-t1 (cdr math-t1))) (and math-t2 (math-normalize math-t2))))) ((= (length expr) 3) (and (symbolp (car expr)) (setq math-t1 (get (car expr) 'math-integral-2)) (progn (while (and math-t1 (not (setq math-t2 (funcall (car math-t1) (nth 1 expr) (nth 2 expr))))) (setq math-t1 (cdr math-t1))) (and math-t2 (math-normalize math-t2)))))) ;; Integral of a rational function. (and (math-ratpoly-p expr math-integ-var) (setq math-t1 (calcFunc-apart expr math-integ-var)) (not (equal math-t1 expr)) (math-integral math-t1)) ;; Try user-defined integration rules. (and math-has-rules (let ((math-old-integ (symbol-function 'calcFunc-integ)) (input (list 'calcFunc-integtry expr math-integ-var)) res part) (unwind-protect (progn (fset 'calcFunc-integ 'math-sub-integration) (setq res (math-rewrite input '(var IntegRules var-IntegRules) 1)) (fset 'calcFunc-integ math-old-integ) (and (not (equal res input)) (if (setq part (math-expr-calls res '(calcFunc-integsubst))) (and (memq (length part) '(3 4 5)) (let ((parts (mapcar (function (lambda (x) (math-expr-subst x (nth 2 part) math-integ-var))) (cdr part)))) (math-integrate-by-substitution expr (car parts) t (or (nth 2 parts) (list 'calcFunc-integfailed math-integ-var)) (nth 3 parts)))) (if (not (math-expr-calls res '(calcFunc-integtry calcFunc-integfailed))) res)))) (fset 'calcFunc-integ math-old-integ)))) ;; See if the function is a symbolic derivative. (and (string-match "'" (symbol-name (car expr))) (let ((name (symbol-name (car expr))) (p expr) (n 0) (which nil) (bad nil)) (while (setq n (1+ n) p (cdr p)) (if (equal (car p) math-integ-var) (if which (setq bad t) (setq which n)) (if (math-expr-contains (car p) math-integ-var) (setq bad t)))) (and which (not bad) (let ((prime (if (= which 1) "'" (format "'%d" which)))) (and (string-match (concat prime "\\('['0-9]*\\|$\\)") name) (cons (intern (concat (substring name 0 (match-beginning 0)) (substring name (+ (match-beginning 0) (length prime))))) (cdr expr))))))) ;; Try transformation methods (parts, substitutions). (and (> math-integ-level 0) (math-do-integral-methods expr)) ;; Try expanding the function's definition. (let ((res (math-expand-formula expr))) (and res (math-integral res)))))) (defun math-sub-integration (expr &rest rest) (or (if (or (not rest) (and (< math-integ-level math-integral-limit) (eq (car rest) math-integ-var))) (math-integral expr) (let ((res (apply math-old-integ expr rest))) (and (or (= math-integ-level math-integral-limit) (not (math-expr-calls res 'calcFunc-integ))) res))) (list 'calcFunc-integfailed expr))) ;; math-so-far is a local variable for math-do-integral-methods, but ;; is used by math-integ-try-linear-substitutions and ;; math-integ-try-substitutions. (defvar math-so-far) ;; math-integ-expr is a local variable for math-do-integral-methods, ;; but is used by math-integ-try-linear-substitutions and ;; math-integ-try-substitutions. (defvar math-integ-expr) (defun math-do-integral-methods (math-integ-expr) (let ((math-so-far math-integ-var-list-list) rat-in) ;; Integration by substitution, for various likely sub-expressions. ;; (In first pass, we look only for sub-exprs that are linear in X.) (or (math-integ-try-linear-substitutions math-integ-expr) (math-integ-try-substitutions math-integ-expr) ;; If function has sines and cosines, try tan(x/2) substitution. (and (let ((p (setq rat-in (math-expr-rational-in math-integ-expr)))) (while (and p (memq (car (car p)) '(calcFunc-sin calcFunc-cos calcFunc-tan calcFunc-sec calcFunc-csc calcFunc-cot)) (equal (nth 1 (car p)) math-integ-var)) (setq p (cdr p))) (null p)) (or (and (math-integ-parts-easy math-integ-expr) (math-integ-try-parts math-integ-expr t)) (math-integrate-by-good-substitution math-integ-expr (list 'calcFunc-tan (math-div math-integ-var 2))))) ;; If function has sinh and cosh, try tanh(x/2) substitution. (and (let ((p rat-in)) (while (and p (memq (car (car p)) '(calcFunc-sinh calcFunc-cosh calcFunc-tanh calcFunc-sech calcFunc-csch calcFunc-coth calcFunc-exp)) (equal (nth 1 (car p)) math-integ-var)) (setq p (cdr p))) (null p)) (or (and (math-integ-parts-easy math-integ-expr) (math-integ-try-parts math-integ-expr t)) (math-integrate-by-good-substitution math-integ-expr (list 'calcFunc-tanh (math-div math-integ-var 2))))) ;; If function has square roots, try sin, tan, or sec substitution. (and (let ((p rat-in)) (setq math-t1 nil) (while (and p (or (equal (car p) math-integ-var) (and (eq (car (car p)) 'calcFunc-sqrt) (setq math-t1 (math-is-polynomial (nth 1 (setq math-t2 (car p))) math-integ-var 2))))) (setq p (cdr p))) (and (null p) math-t1)) (if (cdr (cdr math-t1)) (if (math-guess-if-neg (nth 2 math-t1)) (let* ((c (math-sqrt (math-neg (nth 2 math-t1)))) (d (math-div (nth 1 math-t1) (math-mul -2 c))) (a (math-sqrt (math-add (car math-t1) (math-sqr d))))) (math-integrate-by-good-substitution math-integ-expr (list 'calcFunc-arcsin (math-div-thru (math-add (math-mul c math-integ-var) d) a)))) (let* ((c (math-sqrt (nth 2 math-t1))) (d (math-div (nth 1 math-t1) (math-mul 2 c))) (aa (math-sub (car math-t1) (math-sqr d)))) (if (and nil (not (and (eq d 0) (eq c 1)))) (math-integrate-by-good-substitution math-integ-expr (math-add (math-mul c math-integ-var) d)) (if (math-guess-if-neg aa) (math-integrate-by-good-substitution math-integ-expr (list 'calcFunc-arccosh (math-div-thru (math-add (math-mul c math-integ-var) d) (math-sqrt (math-neg aa))))) (math-integrate-by-good-substitution math-integ-expr (list 'calcFunc-arcsinh (math-div-thru (math-add (math-mul c math-integ-var) d) (math-sqrt aa)))))))) (math-integrate-by-good-substitution math-integ-expr math-t2)) ) ;; Try integration by parts. (math-integ-try-parts math-integ-expr) ;; Give up. nil))) (defun math-integ-parts-easy (expr) (cond ((Math-primp expr) t) ((memq (car expr) '(+ - *)) (and (math-integ-parts-easy (nth 1 expr)) (math-integ-parts-easy (nth 2 expr)))) ((eq (car expr) '/) (and (math-integ-parts-easy (nth 1 expr)) (math-atomic-factorp (nth 2 expr)))) ((eq (car expr) '^) (and (natnump (nth 2 expr)) (math-integ-parts-easy (nth 1 expr)))) ((eq (car expr) 'neg) (math-integ-parts-easy (nth 1 expr))) (t t))) ;; math-prev-parts-v is local to calcFunc-integ (as well as ;; math-integrate-by-parts), but is used by math-integ-try-parts. (defvar math-prev-parts-v) ;; math-good-parts is local to calcFunc-integ (as well as ;; math-integ-try-parts), but is used by math-integrate-by-parts. (defvar math-good-parts) (defun math-integ-try-parts (expr &optional math-good-parts) ;; Integration by parts: ;; integ(f(x) g(x),x) = f(x) h(x) - integ(h(x) f'(x),x) ;; where h(x) = integ(g(x),x). (or (let ((exp (calcFunc-expand expr))) (and (not (equal exp expr)) (math-integral exp))) (and (eq (car expr) '*) (let ((first-bad (or (math-polynomial-p (nth 1 expr) math-integ-var) (equal (nth 2 expr) math-prev-parts-v)))) (or (and first-bad ; so try this one first (math-integrate-by-parts (nth 1 expr) (nth 2 expr))) (math-integrate-by-parts (nth 2 expr) (nth 1 expr)) (and (not first-bad) (math-integrate-by-parts (nth 1 expr) (nth 2 expr)))))) (and (eq (car expr) '/) (math-expr-contains (nth 1 expr) math-integ-var) (let ((recip (math-div 1 (nth 2 expr)))) (or (math-integrate-by-parts (nth 1 expr) recip) (math-integrate-by-parts recip (nth 1 expr))))) (and (eq (car expr) '^) (math-integrate-by-parts (math-pow (nth 1 expr) (math-sub (nth 2 expr) 1)) (nth 1 expr))))) (defun math-integrate-by-parts (u vprime) (let ((math-integ-level (if (or math-good-parts (math-polynomial-p u math-integ-var)) math-integ-level (1- math-integ-level))) (math-doing-parts t) v temp) (and (>= math-integ-level 0) (unwind-protect (progn (setcar (cdr math-cur-record) 'parts) (math-tracing-integral "Integrating by parts, u = " (math-format-value u 1000) ", v' = " (math-format-value vprime 1000) "\n") (and (setq v (math-integral vprime)) (setq temp (calcFunc-deriv u math-integ-var nil t)) (setq temp (let ((math-prev-parts-v v)) (math-integral (math-mul v temp) 'yes))) (setq temp (math-sub (math-mul u v) temp)) (if (eq (nth 1 math-cur-record) 'parts) (calcFunc-expand temp) (setq v (list 'var 'PARTS math-cur-record) temp (let (calc-next-why) (math-simplify-extended (math-solve-for (math-sub v temp) 0 v nil))) temp (if (and (eq (car-safe temp) '/) (math-zerop (nth 2 temp))) nil temp))))) (setcar (cdr math-cur-record) 'busy))))) ;;; This tries two different formulations, hoping the algebraic simplifier ;;; will be strong enough to handle at least one. (defun math-integrate-by-substitution (expr u &optional user uinv uinvprime) (and (> math-integ-level 0) (let ((math-integ-level (max (- math-integ-level 2) 0))) (math-integrate-by-good-substitution expr u user uinv uinvprime)))) (defun math-integrate-by-good-substitution (expr u &optional user uinv uinvprime) (let ((math-living-dangerously t) deriv temp) (and (setq uinv (if uinv (math-expr-subst uinv math-integ-var math-integ-var-2) (let (calc-next-why) (math-solve-for u math-integ-var-2 math-integ-var nil)))) (progn (math-tracing-integral "Integrating by substitution, u = " (math-format-value u 1000) "\n") (or (and (setq deriv (calcFunc-deriv u math-integ-var nil (not user))) (setq temp (math-integral (math-expr-subst (math-expr-subst (math-expr-subst (math-div expr deriv) u math-integ-var-2) math-integ-var uinv) math-integ-var-2 math-integ-var) 'yes))) (and (setq deriv (or uinvprime (calcFunc-deriv uinv math-integ-var-2 math-integ-var (not user)))) (setq temp (math-integral (math-mul (math-expr-subst (math-expr-subst (math-expr-subst expr u math-integ-var-2) math-integ-var uinv) math-integ-var-2 math-integ-var) deriv) 'yes))))) (math-simplify-extended (math-expr-subst temp math-integ-var u))))) ;;; Look for substitutions of the form u = a x + b. (defun math-integ-try-linear-substitutions (sub-expr) (setq math-linear-subst-tried t) (and (not (Math-primp sub-expr)) (or (and (not (memq (car sub-expr) '(+ - * / neg))) (not (and (eq (car sub-expr) '^) (integerp (nth 2 sub-expr)))) (math-expr-contains sub-expr math-integ-var) (let ((res nil)) (while (and (setq sub-expr (cdr sub-expr)) (or (not (math-linear-in (car sub-expr) math-integ-var)) (assoc (car sub-expr) math-so-far) (progn (setq math-so-far (cons (list (car sub-expr)) math-so-far)) (not (setq res (math-integrate-by-substitution math-integ-expr (car sub-expr)))))))) res)) (let ((res nil)) (while (and (setq sub-expr (cdr sub-expr)) (not (setq res (math-integ-try-linear-substitutions (car sub-expr)))))) res)))) ;;; Recursively try different substitutions based on various sub-expressions. (defun math-integ-try-substitutions (sub-expr &optional allow-rat) (and (not (Math-primp sub-expr)) (not (assoc sub-expr math-so-far)) (math-expr-contains sub-expr math-integ-var) (or (and (if (and (not (memq (car sub-expr) '(+ - * / neg))) (not (and (eq (car sub-expr) '^) (integerp (nth 2 sub-expr))))) (setq allow-rat t) (prog1 allow-rat (setq allow-rat nil))) (not (eq sub-expr math-integ-expr)) (or (math-integrate-by-substitution math-integ-expr sub-expr) (and (eq (car sub-expr) '^) (integerp (nth 2 sub-expr)) (< (nth 2 sub-expr) 0) (math-integ-try-substitutions (math-pow (nth 1 sub-expr) (- (nth 2 sub-expr))) t)))) (let ((res nil)) (setq math-so-far (cons (list sub-expr) math-so-far)) (while (and (setq sub-expr (cdr sub-expr)) (not (setq res (math-integ-try-substitutions (car sub-expr) allow-rat))))) res)))) ;; The variable math-expr-parts is local to math-expr-rational-in, ;; but is used by math-expr-rational-in-rec (defvar math-expr-parts) (defun math-expr-rational-in (expr) (let ((math-expr-parts nil)) (math-expr-rational-in-rec expr) (mapcar 'car math-expr-parts))) (defun math-expr-rational-in-rec (expr) (cond ((Math-primp expr) (and (equal expr math-integ-var) (not (assoc expr math-expr-parts)) (setq math-expr-parts (cons (list expr) math-expr-parts)))) ((or (memq (car expr) '(+ - * / neg)) (and (eq (car expr) '^) (integerp (nth 2 expr)))) (math-expr-rational-in-rec (nth 1 expr)) (and (nth 2 expr) (math-expr-rational-in-rec (nth 2 expr)))) ((and (eq (car expr) '^) (eq (math-quarter-integer (nth 2 expr)) 2)) (math-expr-rational-in-rec (list 'calcFunc-sqrt (nth 1 expr)))) (t (and (not (assoc expr math-expr-parts)) (math-expr-contains expr math-integ-var) (setq math-expr-parts (cons (list expr) math-expr-parts)))))) (defun math-expr-calls (expr funcs &optional arg-contains) (if (consp expr) (if (or (memq (car expr) funcs) (and (eq (car expr) '^) (eq (car funcs) 'calcFunc-sqrt) (eq (math-quarter-integer (nth 2 expr)) 2))) (and (or (not arg-contains) (math-expr-contains expr arg-contains)) expr) (and (not (Math-primp expr)) (let ((res nil)) (while (and (setq expr (cdr expr)) (not (setq res (math-expr-calls (car expr) funcs arg-contains))))) res))))) (defun math-fix-const-terms (expr except-vars) (cond ((not (math-expr-depends expr except-vars)) 0) ((Math-primp expr) expr) ((eq (car expr) '+) (math-add (math-fix-const-terms (nth 1 expr) except-vars) (math-fix-const-terms (nth 2 expr) except-vars))) ((eq (car expr) '-) (math-sub (math-fix-const-terms (nth 1 expr) except-vars) (math-fix-const-terms (nth 2 expr) except-vars))) (t expr))) ;; Command for debugging the Calculator's symbolic integrator. (defun calc-dump-integral-cache (&optional arg) (interactive "P") (let ((buf (current-buffer))) (unwind-protect (let ((p math-integral-cache) math-cur-record) (display-buffer (get-buffer-create "*Integral Cache*")) (set-buffer (get-buffer "*Integral Cache*")) (erase-buffer) (while p (setq math-cur-record (car p)) (or arg (math-replace-integral-parts math-cur-record)) (insert (math-format-flat-expr (car math-cur-record) 0) " --> " (if (symbolp (nth 1 math-cur-record)) (concat "(" (symbol-name (nth 1 math-cur-record)) ")") (math-format-flat-expr (nth 1 math-cur-record) 0)) "\n") (setq p (cdr p))) (goto-char (point-min))) (set-buffer buf)))) ;; The variable math-max-integral-limit is local to calcFunc-integ, ;; but is used by math-try-integral. (defvar math-max-integral-limit) (defun math-try-integral (expr) (let ((math-integ-level math-integral-limit) (math-integ-depth 0) (math-integ-msg "Working...done") (math-cur-record nil) ; a technicality (math-integrating t) (calc-prefer-frac t) (calc-symbolic-mode t) (math-has-rules (calc-has-rules 'var-IntegRules))) (or (math-integral expr 'yes) (and math-any-substs (setq math-enable-subst t) (math-integral expr 'yes)) (and (> math-max-integral-limit math-integral-limit) (setq math-integral-limit math-max-integral-limit math-integ-level math-integral-limit) (math-integral expr 'yes))))) (defvar var-IntegLimit nil) (defun calcFunc-integ (expr var &optional low high) (cond ;; Do these even if the parts turn out not to be integrable. ((eq (car-safe expr) '+) (math-add (calcFunc-integ (nth 1 expr) var low high) (calcFunc-integ (nth 2 expr) var low high))) ((eq (car-safe expr) '-) (math-sub (calcFunc-integ (nth 1 expr) var low high) (calcFunc-integ (nth 2 expr) var low high))) ((eq (car-safe expr) 'neg) (math-neg (calcFunc-integ (nth 1 expr) var low high))) ((and (eq (car-safe expr) '*) (not (math-expr-contains (nth 1 expr) var))) (math-mul (nth 1 expr) (calcFunc-integ (nth 2 expr) var low high))) ((and (eq (car-safe expr) '*) (not (math-expr-contains (nth 2 expr) var))) (math-mul (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr))) ((and (eq (car-safe expr) '/) (not (math-expr-contains (nth 1 expr) var)) (not (math-equal-int (nth 1 expr) 1))) (math-mul (nth 1 expr) (calcFunc-integ (math-div 1 (nth 2 expr)) var low high))) ((and (eq (car-safe expr) '/) (not (math-expr-contains (nth 2 expr) var))) (math-div (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr))) ((and (eq (car-safe expr) '/) (eq (car-safe (nth 1 expr)) '*) (not (math-expr-contains (nth 1 (nth 1 expr)) var))) (math-mul (nth 1 (nth 1 expr)) (calcFunc-integ (math-div (nth 2 (nth 1 expr)) (nth 2 expr)) var low high))) ((and (eq (car-safe expr) '/) (eq (car-safe (nth 1 expr)) '*) (not (math-expr-contains (nth 2 (nth 1 expr)) var))) (math-mul (nth 2 (nth 1 expr)) (calcFunc-integ (math-div (nth 1 (nth 1 expr)) (nth 2 expr)) var low high))) ((and (eq (car-safe expr) '/) (eq (car-safe (nth 2 expr)) '*) (not (math-expr-contains (nth 1 (nth 2 expr)) var))) (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 2 (nth 2 expr))) var low high) (nth 1 (nth 2 expr)))) ((and (eq (car-safe expr) '/) (eq (car-safe (nth 2 expr)) '*) (not (math-expr-contains (nth 2 (nth 2 expr)) var))) (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 1 (nth 2 expr))) var low high) (nth 2 (nth 2 expr)))) ((eq (car-safe expr) 'vec) (cons 'vec (mapcar (function (lambda (x) (calcFunc-integ x var low high))) (cdr expr)))) (t (let ((state (list calc-angle-mode ;;calc-symbolic-mode ;;calc-prefer-frac calc-internal-prec (calc-var-value 'var-IntegRules) (calc-var-value 'var-IntegSimpRules)))) (or (equal state math-integral-cache-state) (setq math-integral-cache-state state math-integral-cache nil))) (let* ((math-max-integral-limit (or (and (natnump var-IntegLimit) var-IntegLimit) 3)) (math-integral-limit 1) (sexpr (math-expr-subst expr var math-integ-var)) (trace-buffer (get-buffer "*Trace*")) (calc-language (if (eq calc-language 'big) nil calc-language)) (math-any-substs t) (math-enable-subst nil) (math-prev-parts-v nil) (math-doing-parts nil) (math-good-parts nil) (res (if trace-buffer (let ((calcbuf (current-buffer)) (calcwin (selected-window))) (unwind-protect (progn (if (get-buffer-window trace-buffer) (select-window (get-buffer-window trace-buffer))) (set-buffer trace-buffer) (goto-char (point-max)) (or (assq 'scroll-stop (buffer-local-variables)) (progn (make-local-variable 'scroll-step) (setq scroll-step 3))) (insert "\n\n\n") (set-buffer calcbuf) (math-try-integral sexpr)) (select-window calcwin) (set-buffer calcbuf))) (math-try-integral sexpr)))) (if res (progn (if (calc-has-rules 'var-IntegAfterRules) (setq res (math-rewrite res '(var IntegAfterRules var-IntegAfterRules)))) (math-simplify (if (and low high) (math-sub (math-expr-subst res math-integ-var high) (math-expr-subst res math-integ-var low)) (setq res (math-fix-const-terms res math-integ-vars)) (if low (math-expr-subst res math-integ-var low) (math-expr-subst res math-integ-var var))))) (append (list 'calcFunc-integ expr var) (and low (list low)) (and high (list high)))))))) (math-defintegral calcFunc-inv (math-integral (math-div 1 u))) (math-defintegral calcFunc-conj (let ((int (math-integral u))) (and int (list 'calcFunc-conj int)))) (math-defintegral calcFunc-deg (let ((int (math-integral u))) (and int (list 'calcFunc-deg int)))) (math-defintegral calcFunc-rad (let ((int (math-integral u))) (and int (list 'calcFunc-rad int)))) (math-defintegral calcFunc-re (let ((int (math-integral u))) (and int (list 'calcFunc-re int)))) (math-defintegral calcFunc-im (let ((int (math-integral u))) (and int (list 'calcFunc-im int)))) (math-defintegral calcFunc-sqrt (and (equal u math-integ-var) (math-mul '(frac 2 3) (list 'calcFunc-sqrt (math-pow u 3))))) (math-defintegral calcFunc-exp (or (and (equal u math-integ-var) (list 'calcFunc-exp u)) (let ((p (math-is-polynomial u math-integ-var 2))) (and (nth 2 p) (let ((sqa (math-sqrt (math-neg (nth 2 p))))) (math-div (math-mul (math-mul (math-div (list 'calcFunc-sqrt '(var pi var-pi)) sqa) (math-normalize (list 'calcFunc-exp (math-div (math-sub (math-mul (car p) (nth 2 p)) (math-div (math-sqr (nth 1 p)) 4)) (nth 2 p))))) (list 'calcFunc-erf (math-sub (math-mul sqa math-integ-var) (math-div (nth 1 p) (math-mul 2 sqa))))) 2)))))) (math-defintegral calcFunc-ln (or (and (equal u math-integ-var) (math-sub (math-mul u (list 'calcFunc-ln u)) u)) (and (eq (car u) '*) (math-integral (math-add (list 'calcFunc-ln (nth 1 u)) (list 'calcFunc-ln (nth 2 u))))) (and (eq (car u) '/) (math-integral (math-sub (list 'calcFunc-ln (nth 1 u)) (list 'calcFunc-ln (nth 2 u))))) (and (eq (car u) '^) (math-integral (math-mul (nth 2 u) (list 'calcFunc-ln (nth 1 u))))))) (math-defintegral calcFunc-log10 (and (equal u math-integ-var) (math-sub (math-mul u (list 'calcFunc-ln u)) (math-div u (list 'calcFunc-ln 10))))) (math-defintegral-2 calcFunc-log (math-integral (math-div (list 'calcFunc-ln u) (list 'calcFunc-ln v)))) (math-defintegral calcFunc-sin (or (and (equal u math-integ-var) (math-neg (math-from-radians-2 (list 'calcFunc-cos u)))) (and (nth 2 (math-is-polynomial u math-integ-var 2)) (math-integral (math-to-exponentials (list 'calcFunc-sin u)))))) (math-defintegral calcFunc-cos (or (and (equal u math-integ-var) (math-from-radians-2 (list 'calcFunc-sin u))) (and (nth 2 (math-is-polynomial u math-integ-var 2)) (math-integral (math-to-exponentials (list 'calcFunc-cos u)))))) (math-defintegral calcFunc-tan (and (equal u math-integ-var) (math-from-radians-2 (list 'calcFunc-ln (list 'calcFunc-sec u))))) (math-defintegral calcFunc-sec (and (equal u math-integ-var) (math-from-radians-2 (list 'calcFunc-ln (math-add (list 'calcFunc-sec u) (list 'calcFunc-tan u)))))) (math-defintegral calcFunc-csc (and (equal u math-integ-var) (math-from-radians-2 (list 'calcFunc-ln (math-sub (list 'calcFunc-csc u) (list 'calcFunc-cot u)))))) (math-defintegral calcFunc-cot (and (equal u math-integ-var) (math-from-radians-2 (list 'calcFunc-ln (list 'calcFunc-sin u))))) (math-defintegral calcFunc-arcsin (and (equal u math-integ-var) (math-add (math-mul u (list 'calcFunc-arcsin u)) (math-from-radians-2 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))) (math-defintegral calcFunc-arccos (and (equal u math-integ-var) (math-sub (math-mul u (list 'calcFunc-arccos u)) (math-from-radians-2 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))) (math-defintegral calcFunc-arctan (and (equal u math-integ-var) (math-sub (math-mul u (list 'calcFunc-arctan u)) (math-from-radians-2 (math-div (list 'calcFunc-ln (math-add 1 (math-sqr u))) 2))))) (math-defintegral calcFunc-sinh (and (equal u math-integ-var) (list 'calcFunc-cosh u))) (math-defintegral calcFunc-cosh (and (equal u math-integ-var) (list 'calcFunc-sinh u))) (math-defintegral calcFunc-tanh (and (equal u math-integ-var) (list 'calcFunc-ln (list 'calcFunc-cosh u)))) (math-defintegral calcFunc-sech (and (equal u math-integ-var) (list 'calcFunc-arctan (list 'calcFunc-sinh u)))) (math-defintegral calcFunc-csch (and (equal u math-integ-var) (list 'calcFunc-ln (list 'calcFunc-tanh (math-div u 2))))) (math-defintegral calcFunc-coth (and (equal u math-integ-var) (list 'calcFunc-ln (list 'calcFunc-sinh u)))) (math-defintegral calcFunc-arcsinh (and (equal u math-integ-var) (math-sub (math-mul u (list 'calcFunc-arcsinh u)) (list 'calcFunc-sqrt (math-add (math-sqr u) 1))))) (math-defintegral calcFunc-arccosh (and (equal u math-integ-var) (math-sub (math-mul u (list 'calcFunc-arccosh u)) (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))) (math-defintegral calcFunc-arctanh (and (equal u math-integ-var) (math-sub (math-mul u (list 'calcFunc-arctan u)) (math-div (list 'calcFunc-ln (math-add 1 (math-sqr u))) 2)))) ;;; (Ax + B) / (ax^2 + bx + c)^n forms. (math-defintegral-2 / (math-integral-rational-funcs u v)) (defun math-integral-rational-funcs (u v) (let ((pu (math-is-polynomial u math-integ-var 1)) (vpow 1) pv) (and pu (catch 'int-rat (if (and (eq (car-safe v) '^) (natnump (nth 2 v))) (setq vpow (nth 2 v) v (nth 1 v))) (and (setq pv (math-is-polynomial v math-integ-var 2)) (let ((int (math-mul-thru (car pu) (math-integral-q02 (car pv) (nth 1 pv) (nth 2 pv) v vpow)))) (if (cdr pu) (setq int (math-add int (math-mul-thru (nth 1 pu) (math-integral-q12 (car pv) (nth 1 pv) (nth 2 pv) v vpow))))) int)))))) (defun math-integral-q12 (a b c v vpow) (let (q) (cond ((not c) (cond ((= vpow 1) (math-sub (math-div math-integ-var b) (math-mul (math-div a (math-sqr b)) (list 'calcFunc-ln v)))) ((= vpow 2) (math-div (math-add (list 'calcFunc-ln v) (math-div a v)) (math-sqr b))) (t (let ((nm1 (math-sub vpow 1)) (nm2 (math-sub vpow 2))) (math-div (math-sub (math-div a (math-mul nm1 (math-pow v nm1))) (math-div 1 (math-mul nm2 (math-pow v nm2)))) (math-sqr b)))))) ((math-zerop (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b)))) (let ((part (math-div b (math-mul 2 c)))) (math-mul-thru (math-pow c vpow) (math-integral-q12 part 1 nil (math-add math-integ-var part) (* vpow 2))))) ((= vpow 1) (and (math-ratp q) (math-negp q) (let ((calc-symbolic-mode t)) (math-ratp (math-sqrt (math-neg q)))) (throw 'int-rat nil)) ; should have used calcFunc-apart first (math-sub (math-div (list 'calcFunc-ln v) (math-mul 2 c)) (math-mul-thru (math-div b (math-mul 2 c)) (math-integral-q02 a b c v 1)))) (t (let ((n (1- vpow))) (math-sub (math-neg (math-div (math-add (math-mul b math-integ-var) (math-mul 2 a)) (math-mul n (math-mul q (math-pow v n))))) (math-mul-thru (math-div (math-mul b (1- (* 2 n))) (math-mul n q)) (math-integral-q02 a b c v n)))))))) (defun math-integral-q02 (a b c v vpow) (let (q rq part) (cond ((not c) (cond ((= vpow 1) (math-div (list 'calcFunc-ln v) b)) (t (math-div (math-pow v (- 1 vpow)) (math-mul (- 1 vpow) b))))) ((math-zerop (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b)))) (let ((part (math-div b (math-mul 2 c)))) (math-mul-thru (math-pow c vpow) (math-integral-q02 part 1 nil (math-add math-integ-var part) (* vpow 2))))) ((progn (setq part (math-add (math-mul 2 (math-mul c math-integ-var)) b)) (> vpow 1)) (let ((n (1- vpow))) (math-add (math-div part (math-mul n (math-mul q (math-pow v n)))) (math-mul-thru (math-div (math-mul (- (* 4 n) 2) c) (math-mul n q)) (math-integral-q02 a b c v n))))) ((math-guess-if-neg q) (setq rq (list 'calcFunc-sqrt (math-neg q))) ;;(math-div-thru (list 'calcFunc-ln ;; (math-div (math-sub part rq) ;; (math-add part rq))) ;; rq) (math-div (math-mul -2 (list 'calcFunc-arctanh (math-div part rq))) rq)) (t (setq rq (list 'calcFunc-sqrt q)) (math-div (math-mul 2 (math-to-radians-2 (list 'calcFunc-arctan (math-div part rq)))) rq))))) (math-defintegral calcFunc-erf (and (equal u math-integ-var) (math-add (math-mul u (list 'calcFunc-erf u)) (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u)) (list 'calcFunc-sqrt '(var pi var-pi))))))) (math-defintegral calcFunc-erfc (and (equal u math-integ-var) (math-sub (math-mul u (list 'calcFunc-erfc u)) (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u)) (list 'calcFunc-sqrt '(var pi var-pi))))))) (defvar math-tabulate-initial nil) (defvar math-tabulate-function nil) ;; The variables calc-low and calc-high are local to calcFunc-table, ;; but are used by math-scan-for-limits. (defvar calc-low) (defvar calc-high) (defun calcFunc-table (expr var &optional calc-low calc-high step) (or calc-low (setq calc-low '(neg (var inf var-inf)) calc-high '(var inf var-inf))) (or calc-high (setq calc-high calc-low calc-low 1)) (and (or (math-infinitep calc-low) (math-infinitep calc-high)) (not step) (math-scan-for-limits expr)) (and step (math-zerop step) (math-reject-arg step 'nonzerop)) (let ((known (+ (if (Math-objectp calc-low) 1 0) (if (Math-objectp calc-high) 1 0) (if (or (null step) (Math-objectp step)) 1 0))) (count '(var inf var-inf)) vec) (or (= known 2) ; handy optimization (equal calc-high '(var inf var-inf)) (progn (setq count (math-div (math-sub calc-high calc-low) (or step 1))) (or (Math-objectp count) (setq count (math-simplify count))) (if (Math-messy-integerp count) (setq count (math-trunc count))))) (if (Math-negp count) (setq count -1)) (if (integerp count) (let ((var-DUMMY nil) (vec math-tabulate-initial) (math-working-step-2 (1+ count)) (math-working-step 0)) (setq expr (math-evaluate-expr (math-expr-subst expr var '(var DUMMY var-DUMMY)))) (while (>= count 0) (setq math-working-step (1+ math-working-step) var-DUMMY calc-low vec (cond ((eq math-tabulate-function 'calcFunc-sum) (math-add vec (math-evaluate-expr expr))) ((eq math-tabulate-function 'calcFunc-prod) (math-mul vec (math-evaluate-expr expr))) (t (cons (math-evaluate-expr expr) vec))) calc-low (math-add calc-low (or step 1)) count (1- count))) (if math-tabulate-function vec (cons 'vec (nreverse vec)))) (if (Math-integerp count) (calc-record-why 'fixnump calc-high) (if (Math-num-integerp calc-low) (if (Math-num-integerp calc-high) (calc-record-why 'integerp step) (calc-record-why 'integerp calc-high)) (calc-record-why 'integerp calc-low))) (append (list (or math-tabulate-function 'calcFunc-table) expr var) (and (not (and (equal calc-low '(neg (var inf var-inf))) (equal calc-high '(var inf var-inf)))) (list calc-low calc-high)) (and step (list step)))))) (defun math-scan-for-limits (x) (cond ((Math-primp x)) ((and (eq (car x) 'calcFunc-subscr) (Math-vectorp (nth 1 x)) (math-expr-contains (nth 2 x) var)) (let* ((calc-next-why nil) (low-val (math-solve-for (nth 2 x) 1 var nil)) (high-val (math-solve-for (nth 2 x) (1- (length (nth 1 x))) var nil)) temp) (and low-val (math-realp low-val) high-val (math-realp high-val)) (and (Math-lessp high-val low-val) (setq temp low-val low-val high-val high-val temp)) (setq calc-low (math-max calc-low (math-ceiling low-val)) calc-high (math-min calc-high (math-floor high-val))))) (t (while (setq x (cdr x)) (math-scan-for-limits (car x)))))) (defvar math-disable-sums nil) (defun calcFunc-sum (expr var &optional low high step) (if math-disable-sums (math-reject-arg)) (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2))) (math-sum-rec expr var low high step))) (math-disable-sums t)) (math-normalize res))) (defun math-sum-rec (expr var &optional low high step) (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf))) (and low (not high) (setq high low low 1)) (let (t1 t2 val) (setq val (cond ((not (math-expr-contains expr var)) (math-mul expr (math-add (math-div (math-sub high low) (or step 1)) 1))) ((and step (not (math-equal-int step 1))) (if (math-negp step) (math-sum-rec expr var high low (math-neg step)) (let ((lo (math-simplify (math-div low step)))) (if (math-known-num-integerp lo) (math-sum-rec (math-normalize (math-expr-subst expr var (math-mul step var))) var lo (math-simplify (math-div high step))) (math-sum-rec (math-normalize (math-expr-subst expr var (math-add (math-mul step var) low))) var 0 (math-simplify (math-div (math-sub high low) step))))))) ((memq (setq t1 (math-compare low high)) '(0 1)) (if (eq t1 0) (math-expr-subst expr var low) 0)) ((setq t1 (math-is-polynomial expr var 20)) (let ((poly nil) (n 0)) (while t1 (setq poly (math-poly-mix poly 1 (math-sum-integer-power n) (car t1)) n (1+ n) t1 (cdr t1))) (setq n (math-build-polynomial-expr poly high)) (if (= low 1) n (math-sub n (math-build-polynomial-expr poly (math-sub low 1)))))) ((and (memq (car expr) '(+ -)) (setq t1 (math-sum-rec (nth 1 expr) var low high) t2 (math-sum-rec (nth 2 expr) var low high)) (not (and (math-expr-calls t1 '(calcFunc-sum)) (math-expr-calls t2 '(calcFunc-sum))))) (list (car expr) t1 t2)) ((and (eq (car expr) '*) (setq t1 (math-sum-const-factors expr var))) (math-mul (car t1) (math-sum-rec (cdr t1) var low high))) ((and (eq (car expr) '*) (memq (car-safe (nth 1 expr)) '(+ -))) (math-sum-rec (math-add-or-sub (math-mul (nth 1 (nth 1 expr)) (nth 2 expr)) (math-mul (nth 2 (nth 1 expr)) (nth 2 expr)) nil (eq (car (nth 1 expr)) '-)) var low high)) ((and (eq (car expr) '*) (memq (car-safe (nth 2 expr)) '(+ -))) (math-sum-rec (math-add-or-sub (math-mul (nth 1 expr) (nth 1 (nth 2 expr))) (math-mul (nth 1 expr) (nth 2 (nth 2 expr))) nil (eq (car (nth 2 expr)) '-)) var low high)) ((and (eq (car expr) '/) (not (math-primp (nth 1 expr))) (setq t1 (math-sum-const-factors (nth 1 expr) var))) (math-mul (car t1) (math-sum-rec (math-div (cdr t1) (nth 2 expr)) var low high))) ((and (eq (car expr) '/) (setq t1 (math-sum-const-factors (nth 2 expr) var))) (math-div (math-sum-rec (math-div (nth 1 expr) (cdr t1)) var low high) (car t1))) ((eq (car expr) 'neg) (math-neg (math-sum-rec (nth 1 expr) var low high))) ((and (eq (car expr) '^) (not (math-expr-contains (nth 1 expr) var)) (setq t1 (math-is-polynomial (nth 2 expr) var 1))) (let ((x (math-pow (nth 1 expr) (nth 1 t1)))) (math-div (math-mul (math-sub (math-pow x (math-add 1 high)) (math-pow x low)) (math-pow (nth 1 expr) (car t1))) (math-sub x 1)))) ((and (setq t1 (math-to-exponentials expr)) (setq t1 (math-sum-rec t1 var low high)) (not (math-expr-calls t1 '(calcFunc-sum)))) (math-to-exps t1)) ((memq (car expr) '(calcFunc-ln calcFunc-log10)) (list (car expr) (calcFunc-prod (nth 1 expr) var low high))) ((and (eq (car expr) 'calcFunc-log) (= (length expr) 3) (not (math-expr-contains (nth 2 expr) var))) (list 'calcFunc-log (calcFunc-prod (nth 1 expr) var low high) (nth 2 expr))))) (if (equal val '(var nan var-nan)) (setq val nil)) (or val (let* ((math-tabulate-initial 0) (math-tabulate-function 'calcFunc-sum)) (calcFunc-table expr var low high))))) (defun calcFunc-asum (expr var low &optional high step no-mul-flag) (or high (setq high low low 1)) (if (and step (not (math-equal-int step 1))) (if (math-negp step) (math-mul (math-pow -1 low) (calcFunc-asum expr var high low (math-neg step) t)) (let ((lo (math-simplify (math-div low step)))) (if (math-num-integerp lo) (calcFunc-asum (math-normalize (math-expr-subst expr var (math-mul step var))) var lo (math-simplify (math-div high step))) (calcFunc-asum (math-normalize (math-expr-subst expr var (math-add (math-mul step var) low))) var 0 (math-simplify (math-div (math-sub high low) step)))))) (math-mul (if no-mul-flag 1 (math-pow -1 low)) (calcFunc-sum (math-mul (math-pow -1 var) expr) var low high)))) (defun math-sum-const-factors (expr var) (let ((const nil) (not-const nil) (p expr)) (while (eq (car-safe p) '*) (if (math-expr-contains (nth 1 p) var) (setq not-const (cons (nth 1 p) not-const)) (setq const (cons (nth 1 p) const))) (setq p (nth 2 p))) (if (math-expr-contains p var) (setq not-const (cons p not-const)) (setq const (cons p const))) (and const (cons (let ((temp (car const))) (while (setq const (cdr const)) (setq temp (list '* (car const) temp))) temp) (let ((temp (or (car not-const) 1))) (while (setq not-const (cdr not-const)) (setq temp (list '* (car not-const) temp))) temp))))) (defvar math-sum-int-pow-cache (list '(0 1))) ;; Following is from CRC Math Tables, 27th ed, pp. 52-53. (defun math-sum-integer-power (pow) (let ((calc-prefer-frac t) (n (length math-sum-int-pow-cache))) (while (<= n pow) (let* ((new (list 0 0)) (lin new) (pp (cdr (nth (1- n) math-sum-int-pow-cache))) (p 2) (sum 0) q) (while pp (setq q (math-div (car pp) p) new (cons (math-mul q n) new) sum (math-add sum q) p (1+ p) pp (cdr pp))) (setcar lin (math-sub 1 (math-mul n sum))) (setq math-sum-int-pow-cache (nconc math-sum-int-pow-cache (list (nreverse new))) n (1+ n)))) (nth pow math-sum-int-pow-cache))) (defun math-to-exponentials (expr) (and (consp expr) (= (length expr) 2) (let ((x (nth 1 expr)) (pi (if calc-symbolic-mode '(var pi var-pi) (math-pi))) (i (if calc-symbolic-mode '(var i var-i) '(cplx 0 1)))) (cond ((eq (car expr) 'calcFunc-exp) (list '^ '(var e var-e) x)) ((eq (car expr) 'calcFunc-sin) (or (eq calc-angle-mode 'rad) (setq x (list '/ (list '* x pi) 180))) (list '/ (list '- (list '^ '(var e var-e) (list '* x i)) (list '^ '(var e var-e) (list 'neg (list '* x i)))) (list '* 2 i))) ((eq (car expr) 'calcFunc-cos) (or (eq calc-angle-mode 'rad) (setq x (list '/ (list '* x pi) 180))) (list '/ (list '+ (list '^ '(var e var-e) (list '* x i)) (list '^ '(var e var-e) (list 'neg (list '* x i)))) 2)) ((eq (car expr) 'calcFunc-sinh) (list '/ (list '- (list '^ '(var e var-e) x) (list '^ '(var e var-e) (list 'neg x))) 2)) ((eq (car expr) 'calcFunc-cosh) (list '/ (list '+ (list '^ '(var e var-e) x) (list '^ '(var e var-e) (list 'neg x))) 2)) (t nil))))) (defun math-to-exps (expr) (cond (calc-symbolic-mode expr) ((Math-primp expr) (if (equal expr '(var e var-e)) (math-e) expr)) ((and (eq (car expr) '^) (equal (nth 1 expr) '(var e var-e))) (list 'calcFunc-exp (nth 2 expr))) (t (cons (car expr) (mapcar 'math-to-exps (cdr expr)))))) (defvar math-disable-prods nil) (defun calcFunc-prod (expr var &optional low high step) (if math-disable-prods (math-reject-arg)) (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2))) (math-prod-rec expr var low high step))) (math-disable-prods t)) (math-normalize res))) (defun math-prod-rec (expr var &optional low high step) (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf))) (and low (not high) (setq high '(var inf var-inf))) (let (t1 t2 t3 val) (setq val (cond ((not (math-expr-contains expr var)) (math-pow expr (math-add (math-div (math-sub high low) (or step 1)) 1))) ((and step (not (math-equal-int step 1))) (if (math-negp step) (math-prod-rec expr var high low (math-neg step)) (let ((lo (math-simplify (math-div low step)))) (if (math-known-num-integerp lo) (math-prod-rec (math-normalize (math-expr-subst expr var (math-mul step var))) var lo (math-simplify (math-div high step))) (math-prod-rec (math-normalize (math-expr-subst expr var (math-add (math-mul step var) low))) var 0 (math-simplify (math-div (math-sub high low) step))))))) ((and (memq (car expr) '(* /)) (setq t1 (math-prod-rec (nth 1 expr) var low high) t2 (math-prod-rec (nth 2 expr) var low high)) (not (and (math-expr-calls t1 '(calcFunc-prod)) (math-expr-calls t2 '(calcFunc-prod))))) (list (car expr) t1 t2)) ((and (eq (car expr) '^) (not (math-expr-contains (nth 2 expr) var))) (math-pow (math-prod-rec (nth 1 expr) var low high) (nth 2 expr))) ((and (eq (car expr) '^) (not (math-expr-contains (nth 1 expr) var))) (math-pow (nth 1 expr) (calcFunc-sum (nth 2 expr) var low high))) ((eq (car expr) 'sqrt) (math-normalize (list 'calcFunc-sqrt (list 'calcFunc-prod (nth 1 expr) var low high)))) ((eq (car expr) 'neg) (math-mul (math-pow -1 (math-add (math-sub high low) 1)) (math-prod-rec (nth 1 expr) var low high))) ((eq (car expr) 'calcFunc-exp) (list 'calcFunc-exp (calcFunc-sum (nth 1 expr) var low high))) ((and (setq t1 (math-is-polynomial expr var 1)) (setq t2 (cond ((or (and (math-equal-int (nth 1 t1) 1) (setq low (math-simplify (math-add low (car t1))) high (math-simplify (math-add high (car t1))))) (and (math-equal-int (nth 1 t1) -1) (setq t2 low low (math-simplify (math-sub (car t1) high)) high (math-simplify (math-sub (car t1) t2))))) (if (or (math-zerop low) (math-zerop high)) 0 (if (and (or (math-negp low) (math-negp high)) (or (math-num-integerp low) (math-num-integerp high))) (if (math-posp high) 0 (math-mul (math-pow -1 (math-add (math-add low high) 1)) (list '/ (list 'calcFunc-fact (math-neg low)) (list 'calcFunc-fact (math-sub -1 high))))) (list '/ (list 'calcFunc-fact high) (list 'calcFunc-fact (math-sub low 1)))))) ((and (or (and (math-equal-int (nth 1 t1) 2) (setq t2 (math-simplify (math-add (math-mul low 2) (car t1))) t3 (math-simplify (math-add (math-mul high 2) (car t1))))) (and (math-equal-int (nth 1 t1) -2) (setq t2 (math-simplify (math-sub (car t1) (math-mul high 2))) t3 (math-simplify (math-sub (car t1) (math-mul low 2)))))) (or (math-integerp t2) (and (math-messy-integerp t2) (setq t2 (math-trunc t2))) (math-integerp t3) (and (math-messy-integerp t3) (setq t3 (math-trunc t3))))) (if (or (math-zerop t2) (math-zerop t3)) 0 (if (or (math-evenp t2) (math-evenp t3)) (if (or (math-negp t2) (math-negp t3)) (if (math-posp high) 0 (list '/ (list 'calcFunc-dfact (math-neg t2)) (list 'calcFunc-dfact (math-sub -2 t3)))) (list '/ (list 'calcFunc-dfact t3) (list 'calcFunc-dfact (math-sub t2 2)))) (if (math-negp t3) (list '* (list '^ -1 (list '/ (list '- (list '- t2 t3) 2) 2)) (list '/ (list 'calcFunc-dfact (math-neg t2)) (list 'calcFunc-dfact (math-sub -2 t3)))) (if (math-posp t2) (list '/ (list 'calcFunc-dfact t3) (list 'calcFunc-dfact (math-sub t2 2))) nil)))))))) t2))) (if (equal val '(var nan var-nan)) (setq val nil)) (or val (let* ((math-tabulate-initial 1) (math-tabulate-function 'calcFunc-prod)) (calcFunc-table expr var low high))))) (defvar math-solve-ranges nil) (defvar math-solve-sign) ;;; Attempt to reduce math-solve-lhs = math-solve-rhs to ;;; math-solve-var = math-solve-rhs', where math-solve-var appears ;;; in math-solve-lhs but not in math-solve-rhs or math-solve-rhs'; ;;; return math-solve-rhs'. ;;; Uses global values: math-solve-var, math-solve-full. (defvar math-solve-var) (defvar math-solve-full) ;; The variables math-solve-lhs, math-solve-rhs and math-try-solve-sign ;; are local to math-try-solve-for, but are used by math-try-solve-prod. ;; (math-solve-lhs and math-solve-rhs are is also local to ;; math-decompose-poly, but used by math-solve-poly-funny-powers.) (defvar math-solve-lhs) (defvar math-solve-rhs) (defvar math-try-solve-sign) (defun math-try-solve-for (math-solve-lhs math-solve-rhs &optional math-try-solve-sign no-poly) (let (math-t1 math-t2 math-t3) (cond ((equal math-solve-lhs math-solve-var) (setq math-solve-sign math-try-solve-sign) (if (eq math-solve-full 'all) (let ((vec (list 'vec (math-evaluate-expr math-solve-rhs))) newvec var p) (while math-solve-ranges (setq p (car math-solve-ranges) var (car p) newvec (list 'vec)) (while (setq p (cdr p)) (setq newvec (nconc newvec (cdr (math-expr-subst vec var (car p)))))) (setq vec newvec math-solve-ranges (cdr math-solve-ranges))) (math-normalize vec)) math-solve-rhs)) ((Math-primp math-solve-lhs) nil) ((and (eq (car math-solve-lhs) '-) (eq (car-safe (nth 1 math-solve-lhs)) (car-safe (nth 2 math-solve-lhs))) (Math-zerop math-solve-rhs) (= (length (nth 1 math-solve-lhs)) 2) (= (length (nth 2 math-solve-lhs)) 2) (setq math-t1 (get (car (nth 1 math-solve-lhs)) 'math-inverse)) (setq math-t2 (funcall math-t1 '(var SOLVEDUM SOLVEDUM))) (eq (math-expr-contains-count math-t2 '(var SOLVEDUM SOLVEDUM)) 1) (setq math-t3 (math-solve-above-dummy math-t2)) (setq math-t1 (math-try-solve-for (math-sub (nth 1 (nth 1 math-solve-lhs)) (math-expr-subst math-t2 math-t3 (nth 1 (nth 2 math-solve-lhs)))) 0))) math-t1) ((eq (car math-solve-lhs) 'neg) (math-try-solve-for (nth 1 math-solve-lhs) (math-neg math-solve-rhs) (and math-try-solve-sign (- math-try-solve-sign)))) ((and (not (eq math-solve-full 't)) (math-try-solve-prod))) ((and (not no-poly) (setq math-t2 (math-decompose-poly math-solve-lhs math-solve-var 15 math-solve-rhs))) (setq math-t1 (cdr (nth 1 math-t2)) math-t1 (let ((math-solve-ranges math-solve-ranges)) (cond ((= (length math-t1) 5) (apply 'math-solve-quartic (car math-t2) math-t1)) ((= (length math-t1) 4) (apply 'math-solve-cubic (car math-t2) math-t1)) ((= (length math-t1) 3) (apply 'math-solve-quadratic (car math-t2) math-t1)) ((= (length math-t1) 2) (apply 'math-solve-linear (car math-t2) math-try-solve-sign math-t1)) (math-solve-full (math-poly-all-roots (car math-t2) math-t1)) (calc-symbolic-mode nil) (t (math-try-solve-for (car math-t2) (math-poly-any-root (reverse math-t1) 0 t) nil t))))) (if math-t1 (if (eq (nth 2 math-t2) 1) math-t1 (math-solve-prod math-t1 (math-try-solve-for (nth 2 math-t2) 0 nil t))) (calc-record-why "*Unable to find a symbolic solution") nil)) ((and (math-solve-find-root-term math-solve-lhs nil) (eq (math-expr-contains-count math-solve-lhs math-t1) 1)) ; just in case (math-try-solve-for (math-simplify (math-sub (if (or math-t3 (math-evenp math-t2)) (math-pow math-t1 math-t2) (math-neg (math-pow math-t1 math-t2))) (math-expand-power (math-sub (math-normalize (math-expr-subst math-solve-lhs math-t1 0)) math-solve-rhs) math-t2 math-solve-var))) 0)) ((eq (car math-solve-lhs) '+) (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 2 math-solve-lhs) (math-sub math-solve-rhs (nth 1 math-solve-lhs)) math-try-solve-sign)) ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 1 math-solve-lhs) (math-sub math-solve-rhs (nth 2 math-solve-lhs)) math-try-solve-sign)))) ((eq (car math-solve-lhs) 'calcFunc-eq) (math-try-solve-for (math-sub (nth 1 math-solve-lhs) (nth 2 math-solve-lhs)) math-solve-rhs math-try-solve-sign no-poly)) ((eq (car math-solve-lhs) '-) (cond ((or (and (eq (car-safe (nth 1 math-solve-lhs)) 'calcFunc-sin) (eq (car-safe (nth 2 math-solve-lhs)) 'calcFunc-cos)) (and (eq (car-safe (nth 1 math-solve-lhs)) 'calcFunc-cos) (eq (car-safe (nth 2 math-solve-lhs)) 'calcFunc-sin))) (math-try-solve-for (math-sub (nth 1 math-solve-lhs) (list (car (nth 1 math-solve-lhs)) (math-sub (math-quarter-circle t) (nth 1 (nth 2 math-solve-lhs))))) math-solve-rhs)) ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 2 math-solve-lhs) (math-sub (nth 1 math-solve-lhs) math-solve-rhs) (and math-try-solve-sign (- math-try-solve-sign)))) ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 1 math-solve-lhs) (math-add math-solve-rhs (nth 2 math-solve-lhs)) math-try-solve-sign)))) ((and (eq math-solve-full 't) (math-try-solve-prod))) ((and (eq (car math-solve-lhs) '%) (not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var))) (math-try-solve-for (nth 1 math-solve-lhs) (math-add math-solve-rhs (math-solve-get-int (nth 2 math-solve-lhs))))) ((eq (car math-solve-lhs) 'calcFunc-log) (cond ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 1 math-solve-lhs) (math-pow (nth 2 math-solve-lhs) math-solve-rhs))) ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 2 math-solve-lhs) (math-pow (nth 1 math-solve-lhs) (math-div 1 math-solve-rhs)))))) ((and (= (length math-solve-lhs) 2) (symbolp (car math-solve-lhs)) (setq math-t1 (get (car math-solve-lhs) 'math-inverse)) (setq math-t2 (funcall math-t1 math-solve-rhs))) (setq math-t1 (get (car math-solve-lhs) 'math-inverse-sign)) (math-try-solve-for (nth 1 math-solve-lhs) (math-normalize math-t2) (and math-try-solve-sign math-t1 (if (integerp math-t1) (* math-t1 math-try-solve-sign) (funcall math-t1 math-solve-lhs math-try-solve-sign))))) ((and (symbolp (car math-solve-lhs)) (setq math-t1 (get (car math-solve-lhs) 'math-inverse-n)) (setq math-t2 (funcall math-t1 math-solve-lhs math-solve-rhs))) math-t2) ((setq math-t1 (math-expand-formula math-solve-lhs)) (math-try-solve-for math-t1 math-solve-rhs math-try-solve-sign)) (t (calc-record-why "*No inverse known" math-solve-lhs) nil)))) (defun math-try-solve-prod () (cond ((eq (car math-solve-lhs) '*) (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 2 math-solve-lhs) (math-div math-solve-rhs (nth 1 math-solve-lhs)) (math-solve-sign math-try-solve-sign (nth 1 math-solve-lhs)))) ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 1 math-solve-lhs) (math-div math-solve-rhs (nth 2 math-solve-lhs)) (math-solve-sign math-try-solve-sign (nth 2 math-solve-lhs)))) ((Math-zerop math-solve-rhs) (math-solve-prod (let ((math-solve-ranges math-solve-ranges)) (math-try-solve-for (nth 2 math-solve-lhs) 0)) (math-try-solve-for (nth 1 math-solve-lhs) 0))))) ((eq (car math-solve-lhs) '/) (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 2 math-solve-lhs) (math-div (nth 1 math-solve-lhs) math-solve-rhs) (math-solve-sign math-try-solve-sign (nth 1 math-solve-lhs)))) ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 1 math-solve-lhs) (math-mul math-solve-rhs (nth 2 math-solve-lhs)) (math-solve-sign math-try-solve-sign (nth 2 math-solve-lhs)))) ((setq math-t1 (math-try-solve-for (math-sub (nth 1 math-solve-lhs) (math-mul (nth 2 math-solve-lhs) math-solve-rhs)) 0)) math-t1))) ((eq (car math-solve-lhs) '^) (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var)) (math-try-solve-for (nth 2 math-solve-lhs) (math-add (math-normalize (list 'calcFunc-log math-solve-rhs (nth 1 math-solve-lhs))) (math-div (math-mul 2 (math-mul '(var pi var-pi) (math-solve-get-int '(var i var-i)))) (math-normalize (list 'calcFunc-ln (nth 1 math-solve-lhs))))))) ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)) (cond ((and (integerp (nth 2 math-solve-lhs)) (>= (nth 2 math-solve-lhs) 2) (setq math-t1 (math-integer-log2 (nth 2 math-solve-lhs)))) (setq math-t2 math-solve-rhs) (if (and (eq math-solve-full t) (math-known-realp (nth 1 math-solve-lhs))) (progn (while (>= (setq math-t1 (1- math-t1)) 0) (setq math-t2 (list 'calcFunc-sqrt math-t2))) (setq math-t2 (math-solve-get-sign math-t2))) (while (>= (setq math-t1 (1- math-t1)) 0) (setq math-t2 (math-solve-get-sign (math-normalize (list 'calcFunc-sqrt math-t2)))))) (math-try-solve-for (nth 1 math-solve-lhs) (math-normalize math-t2))) ((math-looks-negp (nth 2 math-solve-lhs)) (math-try-solve-for (list '^ (nth 1 math-solve-lhs) (math-neg (nth 2 math-solve-lhs))) (math-div 1 math-solve-rhs))) ((and (eq math-solve-full t) (Math-integerp (nth 2 math-solve-lhs)) (math-known-realp (nth 1 math-solve-lhs))) (setq math-t1 (math-normalize (list 'calcFunc-nroot math-solve-rhs (nth 2 math-solve-lhs)))) (if (math-evenp (nth 2 math-solve-lhs)) (setq math-t1 (math-solve-get-sign math-t1))) (math-try-solve-for (nth 1 math-solve-lhs) math-t1 (and math-try-solve-sign (math-oddp (nth 2 math-solve-lhs)) (math-solve-sign math-try-solve-sign (nth 2 math-solve-lhs))))) (t (math-try-solve-for (nth 1 math-solve-lhs) (math-mul (math-normalize (list 'calcFunc-exp (if (Math-realp (nth 2 math-solve-lhs)) (math-div (math-mul '(var pi var-pi) (math-solve-get-int '(var i var-i) (and (integerp (nth 2 math-solve-lhs)) (math-abs (nth 2 math-solve-lhs))))) (math-div (nth 2 math-solve-lhs) 2)) (math-div (math-mul 2 (math-mul '(var pi var-pi) (math-solve-get-int '(var i var-i) (and (integerp (nth 2 math-solve-lhs)) (math-abs (nth 2 math-solve-lhs)))))) (nth 2 math-solve-lhs))))) (math-normalize (list 'calcFunc-nroot math-solve-rhs (nth 2 math-solve-lhs)))) (and math-try-solve-sign (math-oddp (nth 2 math-solve-lhs)) (math-solve-sign math-try-solve-sign (nth 2 math-solve-lhs))))))))) (t nil))) (defun math-solve-prod (lsoln rsoln) (cond ((null lsoln) rsoln) ((null rsoln) lsoln) ((eq math-solve-full 'all) (cons 'vec (append (cdr lsoln) (cdr rsoln)))) (math-solve-full (list 'calcFunc-if (list 'calcFunc-gt (math-solve-get-sign 1) 0) lsoln rsoln)) (t lsoln))) ;;; This deals with negative, fractional, and symbolic powers of "x". ;; The variable math-solve-b is local to math-decompose-poly, ;; but is used by math-solve-poly-funny-powers. (defvar math-solve-b) (defun math-solve-poly-funny-powers (sub-rhs) ; uses "t1", "t2" (setq math-t1 math-solve-lhs) (let ((pp math-poly-neg-powers) fac) (while pp (setq fac (math-pow (car pp) (or math-poly-mult-powers 1)) math-t1 (math-mul math-t1 fac) math-solve-rhs (math-mul math-solve-rhs fac) pp (cdr pp)))) (if sub-rhs (setq math-t1 (math-sub math-t1 math-solve-rhs))) (let ((math-poly-neg-powers nil)) (setq math-t2 (math-mul (or math-poly-mult-powers 1) (let ((calc-prefer-frac t)) (math-div 1 math-poly-frac-powers))) math-t1 (math-is-polynomial (math-simplify (calcFunc-expand math-t1)) math-solve-b 50)))) ;;; This converts "a x^8 + b x^5 + c x^2" to "(a (x^3)^2 + b (x^3) + c) * x^2". (defun math-solve-crunch-poly (max-degree) ; uses "t1", "t3" (let ((count 0)) (while (and math-t1 (Math-zerop (car math-t1))) (setq math-t1 (cdr math-t1) count (1+ count))) (and math-t1 (let* ((degree (1- (length math-t1))) (scale degree)) (while (and (> scale 1) (= (car math-t3) 1)) (and (= (% degree scale) 0) (let ((p math-t1) (n 0) (new-t1 nil) (okay t)) (while (and p okay) (if (= (% n scale) 0) (setq new-t1 (nconc new-t1 (list (car p)))) (or (Math-zerop (car p)) (setq okay nil))) (setq p (cdr p) n (1+ n))) (if okay (setq math-t3 (cons scale (cdr math-t3)) math-t1 new-t1)))) (setq scale (1- scale))) (setq math-t3 (list (math-mul (car math-t3) math-t2) (math-mul count math-t2))) (<= (1- (length math-t1)) max-degree))))) (defun calcFunc-poly (expr var &optional degree) (if degree (or (natnump degree) (math-reject-arg degree 'fixnatnump)) (setq degree 50)) (let ((p (math-is-polynomial expr var degree 'gen))) (if p (if (equal p '(0)) (list 'vec) (cons 'vec p)) (math-reject-arg expr "Expected a polynomial")))) (defun calcFunc-gpoly (expr var &optional degree) (if degree (or (natnump degree) (math-reject-arg degree 'fixnatnump)) (setq degree 50)) (let* ((math-poly-base-variable var) (d (math-decompose-poly expr var degree nil))) (if d (cons 'vec d) (math-reject-arg expr "Expected a polynomial")))) (defun math-decompose-poly (math-solve-lhs math-solve-var degree sub-rhs) (let ((math-solve-rhs (or sub-rhs 1)) math-t1 math-t2 math-t3) (setq math-t2 (math-polynomial-base math-solve-lhs (function (lambda (math-solve-b) (let ((math-poly-neg-powers '(1)) (math-poly-mult-powers nil) (math-poly-frac-powers 1) (math-poly-exp-base t)) (and (not (equal math-solve-b math-solve-lhs)) (or (not (memq (car-safe math-solve-b) '(+ -))) sub-rhs) (setq math-t3 '(1 0) math-t2 1 math-t1 (math-is-polynomial math-solve-lhs math-solve-b 50)) (if (and (equal math-poly-neg-powers '(1)) (memq math-poly-mult-powers '(nil 1)) (eq math-poly-frac-powers 1) sub-rhs) (setq math-t1 (cons (math-sub (car math-t1) math-solve-rhs) (cdr math-t1))) (math-solve-poly-funny-powers sub-rhs)) (math-solve-crunch-poly degree) (or (math-expr-contains math-solve-b math-solve-var) (math-expr-contains (car math-t3) math-solve-var)))))))) (if math-t2 (list (math-pow math-t2 (car math-t3)) (cons 'vec math-t1) (if sub-rhs (math-pow math-t2 (nth 1 math-t3)) (math-div (math-pow math-t2 (nth 1 math-t3)) math-solve-rhs)))))) (defun math-solve-linear (var sign b a) (math-try-solve-for var (math-div (math-neg b) a) (math-solve-sign sign a) t)) (defun math-solve-quadratic (var c b a) (math-try-solve-for var (if (math-looks-evenp b) (let ((halfb (math-div b 2))) (math-div (math-add (math-neg halfb) (math-solve-get-sign (math-normalize (list 'calcFunc-sqrt (math-add (math-sqr halfb) (math-mul (math-neg c) a)))))) a)) (math-div (math-add (math-neg b) (math-solve-get-sign (math-normalize (list 'calcFunc-sqrt (math-add (math-sqr b) (math-mul 4 (math-mul (math-neg c) a))))))) (math-mul 2 a))) nil t)) (defun math-solve-cubic (var d c b a) (let* ((p (math-div b a)) (q (math-div c a)) (r (math-div d a)) (psqr (math-sqr p)) (aa (math-sub q (math-div psqr 3))) (bb (math-add r (math-div (math-sub (math-mul 2 (math-mul psqr p)) (math-mul 9 (math-mul p q))) 27))) m) (if (Math-zerop aa) (math-try-solve-for (math-pow (math-add var (math-div p 3)) 3) (math-neg bb) nil t) (if (Math-zerop bb) (math-try-solve-for (math-mul (math-add var (math-div p 3)) (math-add (math-sqr (math-add var (math-div p 3))) aa)) 0 nil t) (setq m (math-mul 2 (list 'calcFunc-sqrt (math-div aa -3)))) (math-try-solve-for var (math-sub (math-normalize (math-mul m (list 'calcFunc-cos (math-div (math-sub (list 'calcFunc-arccos (math-div (math-mul 3 bb) (math-mul aa m))) (math-mul 2 (math-mul (math-add 1 (math-solve-get-int 1 3)) (math-half-circle calc-symbolic-mode)))) 3)))) (math-div p 3)) nil t))))) (defun math-solve-quartic (var d c b a aa) (setq a (math-div a aa)) (setq b (math-div b aa)) (setq c (math-div c aa)) (setq d (math-div d aa)) (math-try-solve-for var (let* ((asqr (math-sqr a)) (asqr4 (math-div asqr 4)) (y (let ((math-solve-full nil) calc-next-why) (math-solve-cubic math-solve-var (math-sub (math-sub (math-mul 4 (math-mul b d)) (math-mul asqr d)) (math-sqr c)) (math-sub (math-mul a c) (math-mul 4 d)) (math-neg b) 1))) (rsqr (math-add (math-sub asqr4 b) y)) (r (list 'calcFunc-sqrt rsqr)) (sign1 (math-solve-get-sign 1)) (de (list 'calcFunc-sqrt (math-add (math-sub (math-mul 3 asqr4) (math-mul 2 b)) (if (Math-zerop rsqr) (math-mul 2 (math-mul sign1 (list 'calcFunc-sqrt (math-sub (math-sqr y) (math-mul 4 d))))) (math-sub (math-mul sign1 (math-div (math-sub (math-sub (math-mul 4 (math-mul a b)) (math-mul 8 c)) (math-mul asqr a)) (math-mul 4 r))) rsqr)))))) (math-normalize (math-sub (math-add (math-mul sign1 (math-div r 2)) (math-solve-get-sign (math-div de 2))) (math-div a 4)))) nil t)) (defvar math-symbolic-solve nil) (defvar math-int-coefs nil) ;; The variable math-int-threshold is local to math-poly-all-roots, ;; but is used by math-poly-newton-root. (defvar math-int-threshold) ;; The variables math-int-scale, math-int-factors and math-double-roots ;; are local to math-poly-all-roots, but are used by math-poly-integer-root. (defvar math-int-scale) (defvar math-int-factors) (defvar math-double-roots) (defun math-poly-all-roots (var p &optional math-factoring) (catch 'ouch (let* ((math-symbolic-solve calc-symbolic-mode) (roots nil) (deg (1- (length p))) (orig-p (reverse p)) (math-int-coefs nil) (math-int-scale nil) (math-double-roots nil) (math-int-factors nil) (math-int-threshold nil) (pp p)) ;; If rational coefficients, look for exact rational factors. (while (and pp (Math-ratp (car pp))) (setq pp (cdr pp))) (if pp (if (or math-factoring math-symbolic-solve) (throw 'ouch nil)) (let ((lead (car orig-p)) (calc-prefer-frac t) (scale (apply 'math-lcm-denoms p))) (setq math-int-scale (math-abs (math-mul scale lead)) math-int-threshold (math-div '(float 5 -2) math-int-scale) math-int-coefs (cdr (math-div (cons 'vec orig-p) lead))))) (if (> deg 4) (let ((calc-prefer-frac nil) (calc-symbolic-mode nil) (pp p) (def-p (copy-sequence orig-p))) (while pp (if (Math-numberp (car pp)) (setq pp (cdr pp)) (throw 'ouch nil))) (while (> deg (if math-symbolic-solve 2 4)) (let* ((x (math-poly-any-root def-p '(float 0 0) nil)) b c pp) (if (and (eq (car-safe x) 'cplx) (math-nearly-zerop (nth 2 x) (nth 1 x))) (setq x (calcFunc-re x))) (or math-factoring (setq roots (cons x roots))) (or (math-numberp x) (setq x (math-evaluate-expr x))) (setq pp def-p b (car def-p)) (while (setq pp (cdr pp)) (setq c (car pp)) (setcar pp b) (setq b (math-add (math-mul x b) c))) (setq def-p (cdr def-p) deg (1- deg)))) (setq p (reverse def-p)))) (if (> deg 1) (let ((math-solve-var '(var DUMMY var-DUMMY)) (math-solve-sign nil) (math-solve-ranges nil) (math-solve-full 'all)) (if (= (length p) (length math-int-coefs)) (setq p (reverse math-int-coefs))) (setq roots (append (cdr (apply (cond ((= deg 2) 'math-solve-quadratic) ((= deg 3) 'math-solve-cubic) (t 'math-solve-quartic)) math-solve-var p)) roots))) (if (> deg 0) (setq roots (cons (math-div (math-neg (car p)) (nth 1 p)) roots)))) (if math-factoring (progn (while roots (math-poly-integer-root (car roots)) (setq roots (cdr roots))) (list math-int-factors (nreverse math-int-coefs) math-int-scale)) (let ((vec nil) res) (while roots (let ((root (car roots)) (math-solve-full (and math-solve-full 'all))) (if (math-floatp root) (setq root (math-poly-any-root orig-p root t))) (setq vec (append vec (cdr (or (math-try-solve-for var root nil t) (throw 'ouch nil)))))) (setq roots (cdr roots))) (setq vec (cons 'vec (nreverse vec))) (if math-symbolic-solve (setq vec (math-normalize vec))) (if (eq math-solve-full t) (list 'calcFunc-subscr vec (math-solve-get-int 1 (1- (length orig-p)) 1)) vec)))))) (defun math-lcm-denoms (&rest fracs) (let ((den 1)) (while fracs (if (eq (car-safe (car fracs)) 'frac) (setq den (calcFunc-lcm den (nth 2 (car fracs))))) (setq fracs (cdr fracs))) den)) (defun math-poly-any-root (p x polish) ; p is a reverse poly coeff list (let* ((newt (if (math-zerop x) (math-poly-newton-root p '(cplx (float 123 -6) (float 1 -4)) 4) (math-poly-newton-root p x 4))) (res (if (math-zerop (cdr newt)) (car newt) (if (and (math-lessp (cdr newt) '(float 1 -3)) (not polish)) (setq newt (math-poly-newton-root p (car newt) 30))) (if (math-zerop (cdr newt)) (car newt) (math-poly-laguerre-root p x polish))))) (and math-symbolic-solve (math-floatp res) (throw 'ouch nil)) res)) (defun math-poly-newton-root (p x iters) (let* ((calc-prefer-frac nil) (calc-symbolic-mode nil) (try-integer math-int-coefs) (dx x) b d) (while (and (> (setq iters (1- iters)) 0) (let ((pp p)) (math-working "newton" x) (setq b (car p) d 0) (while (setq pp (cdr pp)) (setq d (math-add (math-mul x d) b) b (math-add (math-mul x b) (car pp)))) (not (math-zerop d))) (progn (setq dx (math-div b d) x (math-sub x dx)) (if try-integer (let ((adx (math-abs-approx dx))) (and (math-lessp adx math-int-threshold) (let ((iroot (math-poly-integer-root x))) (if iroot (setq x iroot dx 0) (setq try-integer nil)))))) (or (not (or (eq dx 0) (math-nearly-zerop dx (math-abs-approx x)))) (progn (setq dx 0) nil))))) (cons x (if (math-zerop x) 1 (math-div (math-abs-approx dx) (math-abs-approx x)))))) (defun math-poly-integer-root (x) (and (math-lessp (calcFunc-xpon (math-abs-approx x)) calc-internal-prec) math-int-coefs (let* ((calc-prefer-frac t) (xre (calcFunc-re x)) (xim (calcFunc-im x)) (xresq (math-sqr xre)) (ximsq (math-sqr xim))) (if (math-lessp ximsq (calcFunc-scf xresq -1)) ;; Look for linear factor (let* ((rnd (math-div (math-round (math-mul xre math-int-scale)) math-int-scale)) (icp math-int-coefs) (rem (car icp)) (newcoef nil)) (while (setq icp (cdr icp)) (setq newcoef (cons rem newcoef) rem (math-add (car icp) (math-mul rem rnd)))) (and (math-zerop rem) (progn (setq math-int-coefs (nreverse newcoef) math-int-factors (cons (list (math-neg rnd)) math-int-factors)) rnd))) ;; Look for irreducible quadratic factor (let* ((rnd1 (math-div (math-round (math-mul xre (math-mul -2 math-int-scale))) math-int-scale)) (sqscale (math-sqr math-int-scale)) (rnd0 (math-div (math-round (math-mul (math-add xresq ximsq) sqscale)) sqscale)) (rem1 (car math-int-coefs)) (icp (cdr math-int-coefs)) (rem0 (car icp)) (newcoef nil) (found (assoc (list rnd0 rnd1 (math-posp xim)) math-double-roots)) this) (if found (setq math-double-roots (delq found math-double-roots) rem0 0 rem1 0) (while (setq icp (cdr icp)) (setq this rem1 newcoef (cons rem1 newcoef) rem1 (math-sub rem0 (math-mul this rnd1)) rem0 (math-sub (car icp) (math-mul this rnd0))))) (and (math-zerop rem0) (math-zerop rem1) (let ((aa (math-div rnd1 -2))) (or found (setq math-int-coefs (reverse newcoef) math-double-roots (cons (list (list rnd0 rnd1 (math-negp xim))) math-double-roots) math-int-factors (cons (cons rnd0 rnd1) math-int-factors))) (math-add aa (let ((calc-symbolic-mode math-symbolic-solve)) (math-mul (math-sqrt (math-sub (math-sqr aa) rnd0)) (if (math-negp xim) -1 1))))))))))) ;;; The following routine is from Numerical Recipes, section 9.5. (defun math-poly-laguerre-root (p x polish) (let* ((calc-prefer-frac nil) (calc-symbolic-mode nil) (iters 0) (m (1- (length p))) (try-newt (not polish)) (tried-newt nil) b d f x1 dx dxold) (while (and (or (< (setq iters (1+ iters)) 50) (math-reject-arg x "*Laguerre's method failed to converge")) (let ((err (math-abs-approx (car p))) (abx (math-abs-approx x)) (pp p)) (setq b (car p) d 0 f 0) (while (setq pp (cdr pp)) (setq f (math-add (math-mul x f) d) d (math-add (math-mul x d) b) b (math-add (math-mul x b) (car pp)) err (math-add (math-abs-approx b) (math-mul abx err)))) (math-lessp (calcFunc-scf err (- -2 calc-internal-prec)) (math-abs-approx b))) (or (not (math-zerop d)) (not (math-zerop f)) (progn (setq x (math-pow (math-neg b) (list 'frac 1 m))) nil)) (let* ((g (math-div d b)) (g2 (math-sqr g)) (h (math-sub g2 (math-mul 2 (math-div f b)))) (sq (math-sqrt (math-mul (1- m) (math-sub (math-mul m h) g2)))) (gp (math-add g sq)) (gm (math-sub g sq))) (if (math-lessp (calcFunc-abssqr gp) (calcFunc-abssqr gm)) (setq gp gm)) (setq dx (math-div m gp) x1 (math-sub x dx)) (if (and try-newt (math-lessp (math-abs-approx dx) (calcFunc-scf (math-abs-approx x) -3))) (let ((newt (math-poly-newton-root p x1 7))) (setq tried-newt t try-newt nil) (if (math-zerop (cdr newt)) (setq x (car newt) x1 x) (if (math-lessp (cdr newt) '(float 1 -6)) (let ((newt2 (math-poly-newton-root p (car newt) 20))) (if (math-zerop (cdr newt2)) (setq x (car newt2) x1 x) (setq x (car newt)))))))) (not (or (eq x x1) (math-nearly-equal x x1)))) (let ((cdx (math-abs-approx dx))) (setq x x1 tried-newt nil) (prog1 (or (<= iters 6) (math-lessp cdx dxold) (progn (if polish (let ((digs (calcFunc-xpon (math-div (math-abs-approx x) cdx)))) (calc-record-why "*Could not attain full precision") (if (natnump digs) (let ((calc-internal-prec (max 3 digs))) (setq x (math-normalize x)))))) nil)) (setq dxold cdx))) (or polish (math-lessp (calcFunc-scf (math-abs-approx x) (- calc-internal-prec)) dxold)))) (or (and (math-floatp x) (math-poly-integer-root x)) x))) (defun math-solve-above-dummy (x) (and (not (Math-primp x)) (if (and (equal (nth 1 x) '(var SOLVEDUM SOLVEDUM)) (= (length x) 2)) x (let ((res nil)) (while (and (setq x (cdr x)) (not (setq res (math-solve-above-dummy (car x)))))) res)))) (defun math-solve-find-root-term (x neg) ; sets "t2", "t3" (if (math-solve-find-root-in-prod x) (setq math-t3 neg math-t1 x) (and (memq (car-safe x) '(+ -)) (or (math-solve-find-root-term (nth 1 x) neg) (math-solve-find-root-term (nth 2 x) (if (eq (car x) '-) (not neg) neg)))))) (defun math-solve-find-root-in-prod (x) (and (consp x) (math-expr-contains x math-solve-var) (or (and (eq (car x) 'calcFunc-sqrt) (setq math-t2 2)) (and (eq (car x) '^) (or (and (memq (math-quarter-integer (nth 2 x)) '(1 2 3)) (setq math-t2 2)) (and (eq (car-safe (nth 2 x)) 'frac) (eq (nth 2 (nth 2 x)) 3) (setq math-t2 3)))) (and (memq (car x) '(* /)) (or (and (not (math-expr-contains (nth 1 x) math-solve-var)) (math-solve-find-root-in-prod (nth 2 x))) (and (not (math-expr-contains (nth 2 x) math-solve-var)) (math-solve-find-root-in-prod (nth 1 x)))))))) ;; The variable math-solve-vars is local to math-solve-system, ;; but is used by math-solve-system-rec. (defvar math-solve-vars) ;; The variable math-solve-simplifying is local to math-solve-system ;; and math-solve-system-rec, but is used by math-solve-system-subst. (defvar math-solve-simplifying) (defun math-solve-system (exprs math-solve-vars math-solve-full) (setq exprs (mapcar 'list (if (Math-vectorp exprs) (cdr exprs) (list exprs))) math-solve-vars (if (Math-vectorp math-solve-vars) (cdr math-solve-vars) (list math-solve-vars))) (or (let ((math-solve-simplifying nil)) (math-solve-system-rec exprs math-solve-vars nil)) (let ((math-solve-simplifying t)) (math-solve-system-rec exprs math-solve-vars nil)))) ;;; The following backtracking solver works by choosing a variable ;;; and equation, and trying to solve the equation for the variable. ;;; If it succeeds it calls itself recursively with that variable and ;;; equation removed from their respective lists, and with the solution ;;; added to solns as well as being substituted into all existing ;;; equations. The algorithm terminates when any solution path ;;; manages to remove all the variables from var-list. ;;; To support calcFunc-roots, entries in eqn-list and solns are ;;; actually lists of equations. ;; The variables math-solve-system-res and math-solve-system-vv are ;; local to math-solve-system-rec, but are used by math-solve-system-subst. (defvar math-solve-system-vv) (defvar math-solve-system-res) (defun math-solve-system-rec (eqn-list var-list solns) (if var-list (let ((v var-list) (math-solve-system-res nil)) ;; Try each variable in turn. (while (and v (let* ((math-solve-system-vv (car v)) (e eqn-list) (elim (eq (car-safe math-solve-system-vv) 'calcFunc-elim))) (if elim (setq math-solve-system-vv (nth 1 math-solve-system-vv))) ;; Try each equation in turn. (while (and e (let ((e2 (car e)) (eprev nil) res2) (setq math-solve-system-res nil) ;; Try to solve for math-solve-system-vv the list of equations e2. (while (and e2 (setq res2 (or (and (eq (car e2) eprev) res2) (math-solve-for (car e2) 0 math-solve-system-vv math-solve-full)))) (setq eprev (car e2) math-solve-system-res (cons (if (eq math-solve-full 'all) (cdr res2) (list res2)) math-solve-system-res) e2 (cdr e2))) (if e2 (setq math-solve-system-res nil) ;; Found a solution. Now try other variables. (setq math-solve-system-res (nreverse math-solve-system-res) math-solve-system-res (math-solve-system-rec (mapcar 'math-solve-system-subst (delq (car e) (copy-sequence eqn-list))) (delq (car v) (copy-sequence var-list)) (let ((math-solve-simplifying nil) (s (mapcar (function (lambda (x) (cons (car x) (math-solve-system-subst (cdr x))))) solns))) (if elim s (cons (cons math-solve-system-vv (apply 'append math-solve-system-res)) s))))) (not math-solve-system-res)))) (setq e (cdr e))) (not math-solve-system-res))) (setq v (cdr v))) math-solve-system-res) ;; Eliminated all variables, so now put solution into the proper format. (setq solns (sort solns (function (lambda (x y) (not (memq (car x) (memq (car y) math-solve-vars))))))) (if (eq math-solve-full 'all) (math-transpose (math-normalize (cons 'vec (if solns (mapcar (function (lambda (x) (cons 'vec (cdr x)))) solns) (mapcar (function (lambda (x) (cons 'vec x))) eqn-list))))) (math-normalize (cons 'vec (if solns (mapcar (function (lambda (x) (cons 'calcFunc-eq x))) solns) (mapcar 'car eqn-list))))))) (defun math-solve-system-subst (x) ; uses "res" and "v" (let ((accum nil) (res2 math-solve-system-res)) (while x (setq accum (nconc accum (mapcar (function (lambda (r) (if math-solve-simplifying (math-simplify (math-expr-subst (car x) math-solve-system-vv r)) (math-expr-subst (car x) math-solve-system-vv r)))) (car res2))) x (cdr x) res2 (cdr res2))) accum)) ;; calc-command-flags is declared in calc.el (defvar calc-command-flags) (defun math-get-from-counter (name) (let ((ctr (assq name calc-command-flags))) (if ctr (setcdr ctr (1+ (cdr ctr))) (setq ctr (cons name 1) calc-command-flags (cons ctr calc-command-flags))) (cdr ctr))) (defvar var-GenCount) (defun math-solve-get-sign (val) (setq val (math-simplify val)) (if (and (eq (car-safe val) '*) (Math-numberp (nth 1 val))) (list '* (nth 1 val) (math-solve-get-sign (nth 2 val))) (and (eq (car-safe val) 'calcFunc-sqrt) (eq (car-safe (nth 1 val)) '^) (setq val (math-normalize (list '^ (nth 1 (nth 1 val)) (math-div (nth 2 (nth 1 val)) 2))))) (if math-solve-full (if (and (calc-var-value 'var-GenCount) (Math-natnump var-GenCount) (not (eq math-solve-full 'all))) (prog1 (math-mul (list 'calcFunc-as var-GenCount) val) (setq var-GenCount (math-add var-GenCount 1)) (calc-refresh-evaltos 'var-GenCount)) (let* ((var (concat "s" (int-to-string (math-get-from-counter 'solve-sign)))) (var2 (list 'var (intern var) (intern (concat "var-" var))))) (if (eq math-solve-full 'all) (setq math-solve-ranges (cons (list var2 1 -1) math-solve-ranges))) (math-mul var2 val))) (calc-record-why "*Choosing positive solution") val))) (defun math-solve-get-int (val &optional range first) (if math-solve-full (if (and (calc-var-value 'var-GenCount) (Math-natnump var-GenCount) (not (eq math-solve-full 'all))) (prog1 (math-mul val (list 'calcFunc-an var-GenCount)) (setq var-GenCount (math-add var-GenCount 1)) (calc-refresh-evaltos 'var-GenCount)) (let* ((var (concat "n" (int-to-string (math-get-from-counter 'solve-int)))) (var2 (list 'var (intern var) (intern (concat "var-" var))))) (if (and range (eq math-solve-full 'all)) (setq math-solve-ranges (cons (cons var2 (cdr (calcFunc-index range (or first 0)))) math-solve-ranges))) (math-mul val var2))) (calc-record-why "*Choosing 0 for arbitrary integer in solution") 0)) (defun math-solve-sign (sign expr) (and sign (let ((s1 (math-possible-signs expr))) (cond ((memq s1 '(4 6)) sign) ((memq s1 '(1 3)) (- sign)))))) (defun math-looks-evenp (expr) (if (Math-integerp expr) (math-evenp expr) (if (memq (car expr) '(* /)) (math-looks-evenp (nth 1 expr))))) (defun math-solve-for (lhs rhs math-solve-var math-solve-full &optional sign) (if (math-expr-contains rhs math-solve-var) (math-solve-for (math-sub lhs rhs) 0 math-solve-var math-solve-full) (and (math-expr-contains lhs math-solve-var) (math-with-extra-prec 1 (let* ((math-poly-base-variable math-solve-var) (res (math-try-solve-for lhs rhs sign))) (if (and (eq math-solve-full 'all) (math-known-realp math-solve-var)) (let ((old-len (length res)) new-len) (setq res (delq nil (mapcar (function (lambda (x) (and (not (memq (car-safe x) '(cplx polar))) x))) res)) new-len (length res)) (if (< new-len old-len) (calc-record-why (if (= new-len 1) "*All solutions were complex" (format "*Omitted %d complex solutions" (- old-len new-len))))))) res))))) (defun math-solve-eqn (expr var full) (if (memq (car-safe expr) '(calcFunc-neq calcFunc-lt calcFunc-gt calcFunc-leq calcFunc-geq)) (let ((res (math-solve-for (cons '- (cdr expr)) 0 var full (if (eq (car expr) 'calcFunc-neq) nil 1)))) (and res (if (eq math-solve-sign 1) (list (car expr) var res) (if (eq math-solve-sign -1) (list (car expr) res var) (or (eq (car expr) 'calcFunc-neq) (calc-record-why "*Can't determine direction of inequality")) (and (memq (car expr) '(calcFunc-neq calcFunc-lt calcFunc-gt)) (list 'calcFunc-neq var res)))))) (let ((res (math-solve-for expr 0 var full))) (and res (list 'calcFunc-eq var res))))) (defun math-reject-solution (expr var func) (if (math-expr-contains expr var) (or (equal (car calc-next-why) '(* "Unable to find a symbolic solution")) (calc-record-why "*Unable to find a solution"))) (list func expr var)) (defun calcFunc-solve (expr var) (or (if (or (Math-vectorp expr) (Math-vectorp var)) (math-solve-system expr var nil) (math-solve-eqn expr var nil)) (math-reject-solution expr var 'calcFunc-solve))) (defun calcFunc-fsolve (expr var) (or (if (or (Math-vectorp expr) (Math-vectorp var)) (math-solve-system expr var t) (math-solve-eqn expr var t)) (math-reject-solution expr var 'calcFunc-fsolve))) (defun calcFunc-roots (expr var) (let ((math-solve-ranges nil)) (or (if (or (Math-vectorp expr) (Math-vectorp var)) (math-solve-system expr var 'all) (math-solve-for expr 0 var 'all)) (math-reject-solution expr var 'calcFunc-roots)))) (defun calcFunc-finv (expr var) (let ((res (math-solve-for expr math-integ-var var nil))) (if res (math-normalize (math-expr-subst res math-integ-var var)) (math-reject-solution expr var 'calcFunc-finv)))) (defun calcFunc-ffinv (expr var) (let ((res (math-solve-for expr math-integ-var var t))) (if res (math-normalize (math-expr-subst res math-integ-var var)) (math-reject-solution expr var 'calcFunc-finv)))) (put 'calcFunc-inv 'math-inverse (function (lambda (x) (math-div 1 x)))) (put 'calcFunc-inv 'math-inverse-sign -1) (put 'calcFunc-sqrt 'math-inverse (function (lambda (x) (math-sqr x)))) (put 'calcFunc-conj 'math-inverse (function (lambda (x) (list 'calcFunc-conj x)))) (put 'calcFunc-abs 'math-inverse (function (lambda (x) (math-solve-get-sign x)))) (put 'calcFunc-deg 'math-inverse (function (lambda (x) (list 'calcFunc-rad x)))) (put 'calcFunc-deg 'math-inverse-sign 1) (put 'calcFunc-rad 'math-inverse (function (lambda (x) (list 'calcFunc-deg x)))) (put 'calcFunc-rad 'math-inverse-sign 1) (put 'calcFunc-ln 'math-inverse (function (lambda (x) (list 'calcFunc-exp x)))) (put 'calcFunc-ln 'math-inverse-sign 1) (put 'calcFunc-log10 'math-inverse (function (lambda (x) (list 'calcFunc-exp10 x)))) (put 'calcFunc-log10 'math-inverse-sign 1) (put 'calcFunc-lnp1 'math-inverse (function (lambda (x) (list 'calcFunc-expm1 x)))) (put 'calcFunc-lnp1 'math-inverse-sign 1) (put 'calcFunc-exp 'math-inverse (function (lambda (x) (math-add (math-normalize (list 'calcFunc-ln x)) (math-mul 2 (math-mul '(var pi var-pi) (math-solve-get-int '(var i var-i)))))))) (put 'calcFunc-exp 'math-inverse-sign 1) (put 'calcFunc-expm1 'math-inverse (function (lambda (x) (math-add (math-normalize (list 'calcFunc-lnp1 x)) (math-mul 2 (math-mul '(var pi var-pi) (math-solve-get-int '(var i var-i)))))))) (put 'calcFunc-expm1 'math-inverse-sign 1) (put 'calcFunc-sin 'math-inverse (function (lambda (x) (let ((n (math-solve-get-int 1))) (math-add (math-mul (math-normalize (list 'calcFunc-arcsin x)) (math-pow -1 n)) (math-mul (math-half-circle t) n)))))) (put 'calcFunc-cos 'math-inverse (function (lambda (x) (math-add (math-solve-get-sign (math-normalize (list 'calcFunc-arccos x))) (math-solve-get-int (math-full-circle t)))))) (put 'calcFunc-tan 'math-inverse (function (lambda (x) (math-add (math-normalize (list 'calcFunc-arctan x)) (math-solve-get-int (math-half-circle t)))))) (put 'calcFunc-arcsin 'math-inverse (function (lambda (x) (math-normalize (list 'calcFunc-sin x))))) (put 'calcFunc-arccos 'math-inverse (function (lambda (x) (math-normalize (list 'calcFunc-cos x))))) (put 'calcFunc-arctan 'math-inverse (function (lambda (x) (math-normalize (list 'calcFunc-tan x))))) (put 'calcFunc-sinh 'math-inverse (function (lambda (x) (let ((n (math-solve-get-int 1))) (math-add (math-mul (math-normalize (list 'calcFunc-arcsinh x)) (math-pow -1 n)) (math-mul (math-half-circle t) (math-mul '(var i var-i) n))))))) (put 'calcFunc-sinh 'math-inverse-sign 1) (put 'calcFunc-cosh 'math-inverse (function (lambda (x) (math-add (math-solve-get-sign (math-normalize (list 'calcFunc-arccosh x))) (math-mul (math-full-circle t) (math-solve-get-int '(var i var-i))))))) (put 'calcFunc-tanh 'math-inverse (function (lambda (x) (math-add (math-normalize (list 'calcFunc-arctanh x)) (math-mul (math-half-circle t) (math-solve-get-int '(var i var-i))))))) (put 'calcFunc-tanh 'math-inverse-sign 1) (put 'calcFunc-arcsinh 'math-inverse (function (lambda (x) (math-normalize (list 'calcFunc-sinh x))))) (put 'calcFunc-arcsinh 'math-inverse-sign 1) (put 'calcFunc-arccosh 'math-inverse (function (lambda (x) (math-normalize (list 'calcFunc-cosh x))))) (put 'calcFunc-arctanh 'math-inverse (function (lambda (x) (math-normalize (list 'calcFunc-tanh x))))) (put 'calcFunc-arctanh 'math-inverse-sign 1) (defun calcFunc-taylor (expr var num) (let ((x0 0) (v var)) (if (memq (car-safe var) '(+ - calcFunc-eq)) (setq x0 (if (eq (car var) '+) (math-neg (nth 2 var)) (nth 2 var)) v (nth 1 var))) (or (and (eq (car-safe v) 'var) (math-expr-contains expr v) (natnump num) (let ((accum (math-expr-subst expr v x0)) (var2 (if (eq (car var) 'calcFunc-eq) (cons '- (cdr var)) var)) (n 0) (nfac 1) (fprime expr)) (while (and (<= (setq n (1+ n)) num) (setq fprime (calcFunc-deriv fprime v nil t))) (setq fprime (math-simplify fprime) nfac (math-mul nfac n) accum (math-add accum (math-div (math-mul (math-pow var2 n) (math-expr-subst fprime v x0)) nfac)))) (and fprime (math-normalize accum)))) (list 'calcFunc-taylor expr var num)))) (provide 'calcalg2) ;;; arch-tag: f2932ec8-dd63-418b-a542-11a644b9d4c4 ;;; calcalg2.el ends here