Mercurial > emacs
view lisp/calc/calcalg2.el @ 58164:4db6f3a44cd6
fix to use new easymenu
author | Vinicius Jose Latorre <viniciusjl@ig.com.br> |
---|---|
date | Fri, 12 Nov 2004 02:14:03 +0000 |
parents | 771eb065ad11 |
children | 7f5b01c17652 cb7f41387eb3 |
line wrap: on
line source
;;; calcalg2.el --- more algebraic functions for Calc ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc. ;; Author: David Gillespie <daveg@synaptics.com> ;; Maintainers: D. Goel <deego@gnufans.org> ;; Colin Walters <walters@debian.org> ;; This file is part of GNU Emacs. ;; GNU Emacs is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY. No author or distributor ;; accepts responsibility to anyone for the consequences of using it ;; or for whether it serves any particular purpose or works at all, ;; unless he says so in writing. Refer to the GNU Emacs General Public ;; License for full details. ;; Everyone is granted permission to copy, modify and redistribute ;; GNU Emacs, but only under the conditions described in the ;; GNU Emacs General Public License. A copy of this license is ;; supposed to have been given to you along with GNU Emacs so you ;; can know your rights and responsibilities. It should be in a ;; file named COPYING. Among other things, the copyright notice ;; and this notice must be preserved on all copies. ;;; Commentary: ;;; Code: ;; This file is autoloaded from calc-ext.el. (require 'calc-ext) (require 'calc-macs) (defun calc-Need-calc-alg-2 () nil) (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) (interactive "sIntegration variable: ") (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 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)))))) (defun math-derivative (expr) ; uses global values: deriv-var, deriv-total. (cond ((equal expr deriv-var) 1) ((or (Math-scalarp expr) (eq (car expr) 'sdev) (and (eq (car expr) 'var) (or (not 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 deriv-symb 'pre-expand)) (let ((exp (math-expand-formula expr))) (and exp (or (let ((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 deriv-symb (throw 'math-deriv nil) (list (if deriv-total 'calcFunc-tderiv 'calcFunc-deriv) expr 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 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 deriv-var &optional deriv-value deriv-symb) (let* ((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 deriv-var deriv-value) res)))) (defun calcFunc-tderiv (expr deriv-var &optional deriv-value deriv-symb) (math-setup-declarations) (let* ((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 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-div 1 (math-sqr (math-normalize (list 'calcFunc-cos 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-div 1 (math-sqr (math-normalize (list 'calcFunc-cosh 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))) 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))) 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) 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) 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)) (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 cur-record = B. (defun math-integral (expr &optional simplify same-as-above) (let* ((simp cur-record) (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 cur-record (progn (math-tracing-integral "Found " (math-format-value (nth 1 cur-record) 1000)) (and (consp (nth 1 cur-record)) (math-replace-integral-parts cur-record)) (math-tracing-integral " => " (math-format-value (nth 1 cur-record) 1000) "\n"))) (or (and cur-record (not (eq (nth 1 cur-record) 'cancelled)) (or (not (integerp (nth 1 cur-record))) (>= (nth 1 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 math-integ-msg))) (if cur-record (setcar (cdr cur-record) (if same-as-above (vector simp) 'busy)) (setq cur-record (list expr (if same-as-above (vector simp) 'busy)) math-integral-cache (cons 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 math-integ-msg))) (setcar (cdr cur-record) (or val (if (or math-enable-subst (not math-any-substs)) math-integ-level 'cancelled))))) (setq val 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)) (defvar math-integral-cache nil) (defvar math-integral-cache-state nil) (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 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.") (defun math-do-integral (expr) (let ((math-linear-subst-tried nil) t1 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 t1 (math-integral (nth 1 expr))) (setq t2 (math-integral (nth 2 expr))) (math-add t1 t2))) ((eq (car expr) '-) (and (setq t1 (math-integral (nth 1 expr))) (setq t2 (math-integral (nth 2 expr))) (math-sub t1 t2))) ((eq (car expr) 'neg) (and (setq t1 (math-integral (nth 1 expr))) (math-neg t1))) ((eq (car expr) '*) (cond ((not (math-expr-contains (nth 1 expr) math-integ-var)) (and (setq t1 (math-integral (nth 2 expr))) (math-mul (nth 1 expr) t1))) ((not (math-expr-contains (nth 2 expr) math-integ-var)) (and (setq t1 (math-integral (nth 1 expr))) (math-mul 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 t1 (math-integral (math-div 1 (nth 2 expr)))) (math-mul (nth 1 expr) t1))) ((not (math-expr-contains (nth 2 expr) math-integ-var)) (and (setq t1 (math-integral (nth 1 expr))) (math-div 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 t1 (math-integral (math-div (nth 2 (nth 1 expr)) (nth 2 expr)))) (math-mul 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 t1 (math-integral (math-div (nth 1 (nth 1 expr)) (nth 2 expr)))) (math-mul 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 t1 (math-integral (math-div (nth 1 expr) (nth 2 (nth 2 expr))))) (math-div 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 t1 (math-integral (math-div (nth 1 expr) (nth 1 (nth 2 expr))))) (math-div 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 t1 (math-is-polynomial (nth 2 expr) math-integ-var 1)) (math-div expr (math-mul (nth 1 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 t1 (math-is-polynomial (nth 1 expr) math-integ-var 1)) (setq t2 (math-add (nth 2 expr) 1)) (math-div (math-pow (nth 1 expr) t2) (math-mul t2 (nth 1 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 t1 (math-is-polynomial expr math-integ-var 20)) (let ((accum 0) (n 1)) (while t1 (if (setq accum (math-add accum (math-div (math-mul (car t1) (math-pow math-integ-var n)) n)) t1 (cdr t1)) (setq n (1+ n)))) accum)) ;; Try looking it up! (cond ((= (length expr) 2) (and (symbolp (car expr)) (setq t1 (get (car expr) 'math-integral)) (progn (while (and t1 (not (setq t2 (funcall (car t1) (nth 1 expr))))) (setq t1 (cdr t1))) (and t2 (math-normalize t2))))) ((= (length expr) 3) (and (symbolp (car expr)) (setq t1 (get (car expr) 'math-integral-2)) (progn (while (and t1 (not (setq t2 (funcall (car t1) (nth 1 expr) (nth 2 expr))))) (setq t1 (cdr t1))) (and t2 (math-normalize t2)))))) ;; Integral of a rational function. (and (math-ratpoly-p expr math-integ-var) (setq t1 (calcFunc-apart expr math-integ-var)) (not (equal t1 expr)) (math-integral t1)) ;; Try user-defined integration rules. (and 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))) (defun math-do-integral-methods (expr) (let ((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 expr) (math-integ-try-substitutions expr) ;; If function has sines and cosines, try tan(x/2) substitution. (and (let ((p (setq rat-in (math-expr-rational-in expr)))) (while (and p (memq (car (car p)) '(calcFunc-sin calcFunc-cos calcFunc-tan)) (equal (nth 1 (car p)) math-integ-var)) (setq p (cdr p))) (null p)) (or (and (math-integ-parts-easy expr) (math-integ-try-parts expr t)) (math-integrate-by-good-substitution 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-exp)) (equal (nth 1 (car p)) math-integ-var)) (setq p (cdr p))) (null p)) (or (and (math-integ-parts-easy expr) (math-integ-try-parts expr t)) (math-integrate-by-good-substitution 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 t1 nil) (while (and p (or (equal (car p) math-integ-var) (and (eq (car (car p)) 'calcFunc-sqrt) (setq t1 (math-is-polynomial (nth 1 (setq t2 (car p))) math-integ-var 2))))) (setq p (cdr p))) (and (null p) t1)) (if (cdr (cdr t1)) (if (math-guess-if-neg (nth 2 t1)) (let* ((c (math-sqrt (math-neg (nth 2 t1)))) (d (math-div (nth 1 t1) (math-mul -2 c))) (a (math-sqrt (math-add (car t1) (math-sqr d))))) (math-integrate-by-good-substitution expr (list 'calcFunc-arcsin (math-div-thru (math-add (math-mul c math-integ-var) d) a)))) (let* ((c (math-sqrt (nth 2 t1))) (d (math-div (nth 1 t1) (math-mul 2 c))) (aa (math-sub (car t1) (math-sqr d)))) (if (and nil (not (and (eq d 0) (eq c 1)))) (math-integrate-by-good-substitution expr (math-add (math-mul c math-integ-var) d)) (if (math-guess-if-neg aa) (math-integrate-by-good-substitution 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 expr (list 'calcFunc-arcsinh (math-div-thru (math-add (math-mul c math-integ-var) d) (math-sqrt aa)))))))) (math-integrate-by-good-substitution expr t2)) ) ;; Try integration by parts. (math-integ-try-parts 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))) (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 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 cur-record) 'parts) (calcFunc-expand temp) (setq v (list 'var 'PARTS cur-record) var-thing (list 'vec (math-sub v temp) v) temp (let (calc-next-why) (math-solve-for (math-sub v temp) 0 v nil))) (and temp (not (integerp temp)) (math-simplify-extended temp))))) (setcar (cdr 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) so-far) (progn (setq so-far (cons (list (car sub-expr)) so-far)) (not (setq res (math-integrate-by-substitution 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 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 expr)) (or (math-integrate-by-substitution 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 so-far (cons (list sub-expr) so-far)) (while (and (setq sub-expr (cdr sub-expr)) (not (setq res (math-integ-try-substitutions (car sub-expr) allow-rat))))) res)))) (defun math-expr-rational-in (expr) (let ((parts nil)) (math-expr-rational-in-rec expr) (mapcar 'car parts))) (defun math-expr-rational-in-rec (expr) (cond ((Math-primp expr) (and (equal expr math-integ-var) (not (assoc expr parts)) (setq parts (cons (list 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 parts)) (math-expr-contains expr math-integ-var) (setq parts (cons (list 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) cur-record) (display-buffer (get-buffer-create "*Integral Cache*")) (set-buffer (get-buffer "*Integral Cache*")) (erase-buffer) (while p (setq cur-record (car p)) (or arg (math-replace-integral-parts cur-record)) (insert (math-format-flat-expr (car cur-record) 0) " --> " (if (symbolp (nth 1 cur-record)) (concat "(" (symbol-name (nth 1 cur-record)) ")") (math-format-flat-expr (nth 1 cur-record) 0)) "\n") (setq p (cdr p))) (goto-char (point-min))) (set-buffer buf)))) (defun math-try-integral (expr) (let ((math-integ-level math-integral-limit) (math-integ-depth 0) (math-integ-msg "Working...done") (cur-record nil) ; a technicality (math-integrating t) (calc-prefer-frac t) (calc-symbolic-mode t) (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))))) (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 (boundp 'var-IntegLimit) (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-neg (math-from-radians-2 (list 'calcFunc-ln (list 'calcFunc-cos 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-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) (defun calcFunc-table (expr var &optional low high step) (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf))) (or high (setq high low low 1)) (and (or (math-infinitep low) (math-infinitep high)) (not step) (math-scan-for-limits expr)) (and step (math-zerop step) (math-reject-arg step 'nonzerop)) (let ((known (+ (if (Math-objectp low) 1 0) (if (Math-objectp high) 1 0) (if (or (null step) (Math-objectp step)) 1 0))) (count '(var inf var-inf)) vec) (or (= known 2) ; handy optimization (equal high '(var inf var-inf)) (progn (setq count (math-div (math-sub high 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 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))) low (math-add low (or step 1)) count (1- count))) (if math-tabulate-function vec (cons 'vec (nreverse vec)))) (if (Math-integerp count) (calc-record-why 'fixnump high) (if (Math-num-integerp low) (if (Math-num-integerp high) (calc-record-why 'integerp step) (calc-record-why 'integerp high)) (calc-record-why 'integerp low))) (append (list (or math-tabulate-function 'calcFunc-table) expr var) (and (not (and (equal low '(neg (var inf var-inf))) (equal high '(var inf var-inf)))) (list low 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 low (math-max low (math-ceiling low-val)) high (math-min 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 (memq low '(0 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) ;;; Attempt to reduce lhs = rhs to solve-var = rhs', where solve-var appears ;;; in lhs but not in rhs or rhs'; return rhs'. ;;; Uses global values: solve-*. (defun math-try-solve-for (lhs rhs &optional sign no-poly) (let (t1 t2 t3) (cond ((equal lhs solve-var) (setq math-solve-sign sign) (if (eq solve-full 'all) (let ((vec (list 'vec (math-evaluate-expr 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)) rhs)) ((Math-primp lhs) nil) ((and (eq (car lhs) '-) (eq (car-safe (nth 1 lhs)) (car-safe (nth 2 lhs))) (Math-zerop rhs) (= (length (nth 1 lhs)) 2) (= (length (nth 2 lhs)) 2) (setq t1 (get (car (nth 1 lhs)) 'math-inverse)) (setq t2 (funcall t1 '(var SOLVEDUM SOLVEDUM))) (eq (math-expr-contains-count t2 '(var SOLVEDUM SOLVEDUM)) 1) (setq t3 (math-solve-above-dummy t2)) (setq t1 (math-try-solve-for (math-sub (nth 1 (nth 1 lhs)) (math-expr-subst t2 t3 (nth 1 (nth 2 lhs)))) 0))) t1) ((eq (car lhs) 'neg) (math-try-solve-for (nth 1 lhs) (math-neg rhs) (and sign (- sign)))) ((and (not (eq solve-full 't)) (math-try-solve-prod))) ((and (not no-poly) (setq t2 (math-decompose-poly lhs solve-var 15 rhs))) (setq t1 (cdr (nth 1 t2)) t1 (let ((math-solve-ranges math-solve-ranges)) (cond ((= (length t1) 5) (apply 'math-solve-quartic (car t2) t1)) ((= (length t1) 4) (apply 'math-solve-cubic (car t2) t1)) ((= (length t1) 3) (apply 'math-solve-quadratic (car t2) t1)) ((= (length t1) 2) (apply 'math-solve-linear (car t2) sign t1)) (solve-full (math-poly-all-roots (car t2) t1)) (calc-symbolic-mode nil) (t (math-try-solve-for (car t2) (math-poly-any-root (reverse t1) 0 t) nil t))))) (if t1 (if (eq (nth 2 t2) 1) t1 (math-solve-prod t1 (math-try-solve-for (nth 2 t2) 0 nil t))) (calc-record-why "*Unable to find a symbolic solution") nil)) ((and (math-solve-find-root-term lhs nil) (eq (math-expr-contains-count lhs t1) 1)) ; just in case (math-try-solve-for (math-simplify (math-sub (if (or t3 (math-evenp t2)) (math-pow t1 t2) (math-neg (math-pow t1 t2))) (math-expand-power (math-sub (math-normalize (math-expr-subst lhs t1 0)) rhs) t2 solve-var))) 0)) ((eq (car lhs) '+) (cond ((not (math-expr-contains (nth 1 lhs) solve-var)) (math-try-solve-for (nth 2 lhs) (math-sub rhs (nth 1 lhs)) sign)) ((not (math-expr-contains (nth 2 lhs) solve-var)) (math-try-solve-for (nth 1 lhs) (math-sub rhs (nth 2 lhs)) sign)))) ((eq (car lhs) 'calcFunc-eq) (math-try-solve-for (math-sub (nth 1 lhs) (nth 2 lhs)) rhs sign no-poly)) ((eq (car lhs) '-) (cond ((or (and (eq (car-safe (nth 1 lhs)) 'calcFunc-sin) (eq (car-safe (nth 2 lhs)) 'calcFunc-cos)) (and (eq (car-safe (nth 1 lhs)) 'calcFunc-cos) (eq (car-safe (nth 2 lhs)) 'calcFunc-sin))) (math-try-solve-for (math-sub (nth 1 lhs) (list (car (nth 1 lhs)) (math-sub (math-quarter-circle t) (nth 1 (nth 2 lhs))))) rhs)) ((not (math-expr-contains (nth 1 lhs) solve-var)) (math-try-solve-for (nth 2 lhs) (math-sub (nth 1 lhs) rhs) (and sign (- sign)))) ((not (math-expr-contains (nth 2 lhs) solve-var)) (math-try-solve-for (nth 1 lhs) (math-add rhs (nth 2 lhs)) sign)))) ((and (eq solve-full 't) (math-try-solve-prod))) ((and (eq (car lhs) '%) (not (math-expr-contains (nth 2 lhs) solve-var))) (math-try-solve-for (nth 1 lhs) (math-add rhs (math-solve-get-int (nth 2 lhs))))) ((eq (car lhs) 'calcFunc-log) (cond ((not (math-expr-contains (nth 2 lhs) solve-var)) (math-try-solve-for (nth 1 lhs) (math-pow (nth 2 lhs) rhs))) ((not (math-expr-contains (nth 1 lhs) solve-var)) (math-try-solve-for (nth 2 lhs) (math-pow (nth 1 lhs) (math-div 1 rhs)))))) ((and (= (length lhs) 2) (symbolp (car lhs)) (setq t1 (get (car lhs) 'math-inverse)) (setq t2 (funcall t1 rhs))) (setq t1 (get (car lhs) 'math-inverse-sign)) (math-try-solve-for (nth 1 lhs) (math-normalize t2) (and sign t1 (if (integerp t1) (* t1 sign) (funcall t1 lhs sign))))) ((and (symbolp (car lhs)) (setq t1 (get (car lhs) 'math-inverse-n)) (setq t2 (funcall t1 lhs rhs))) t2) ((setq t1 (math-expand-formula lhs)) (math-try-solve-for t1 rhs sign)) (t (calc-record-why "*No inverse known" lhs) nil)))) (defun math-try-solve-prod () (cond ((eq (car lhs) '*) (cond ((not (math-expr-contains (nth 1 lhs) solve-var)) (math-try-solve-for (nth 2 lhs) (math-div rhs (nth 1 lhs)) (math-solve-sign sign (nth 1 lhs)))) ((not (math-expr-contains (nth 2 lhs) solve-var)) (math-try-solve-for (nth 1 lhs) (math-div rhs (nth 2 lhs)) (math-solve-sign sign (nth 2 lhs)))) ((Math-zerop rhs) (math-solve-prod (let ((math-solve-ranges math-solve-ranges)) (math-try-solve-for (nth 2 lhs) 0)) (math-try-solve-for (nth 1 lhs) 0))))) ((eq (car lhs) '/) (cond ((not (math-expr-contains (nth 1 lhs) solve-var)) (math-try-solve-for (nth 2 lhs) (math-div (nth 1 lhs) rhs) (math-solve-sign sign (nth 1 lhs)))) ((not (math-expr-contains (nth 2 lhs) solve-var)) (math-try-solve-for (nth 1 lhs) (math-mul rhs (nth 2 lhs)) (math-solve-sign sign (nth 2 lhs)))) ((setq t1 (math-try-solve-for (math-sub (nth 1 lhs) (math-mul (nth 2 lhs) rhs)) 0)) t1))) ((eq (car lhs) '^) (cond ((not (math-expr-contains (nth 1 lhs) solve-var)) (math-try-solve-for (nth 2 lhs) (math-add (math-normalize (list 'calcFunc-log rhs (nth 1 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 lhs))))))) ((not (math-expr-contains (nth 2 lhs) solve-var)) (cond ((and (integerp (nth 2 lhs)) (>= (nth 2 lhs) 2) (setq t1 (math-integer-log2 (nth 2 lhs)))) (setq t2 rhs) (if (and (eq solve-full t) (math-known-realp (nth 1 lhs))) (progn (while (>= (setq t1 (1- t1)) 0) (setq t2 (list 'calcFunc-sqrt t2))) (setq t2 (math-solve-get-sign t2))) (while (>= (setq t1 (1- t1)) 0) (setq t2 (math-solve-get-sign (math-normalize (list 'calcFunc-sqrt t2)))))) (math-try-solve-for (nth 1 lhs) (math-normalize t2))) ((math-looks-negp (nth 2 lhs)) (math-try-solve-for (list '^ (nth 1 lhs) (math-neg (nth 2 lhs))) (math-div 1 rhs))) ((and (eq solve-full t) (Math-integerp (nth 2 lhs)) (math-known-realp (nth 1 lhs))) (setq t1 (math-normalize (list 'calcFunc-nroot rhs (nth 2 lhs)))) (if (math-evenp (nth 2 lhs)) (setq t1 (math-solve-get-sign t1))) (math-try-solve-for (nth 1 lhs) t1 (and sign (math-oddp (nth 2 lhs)) (math-solve-sign sign (nth 2 lhs))))) (t (math-try-solve-for (nth 1 lhs) (math-mul (math-normalize (list 'calcFunc-exp (if (Math-realp (nth 2 lhs)) (math-div (math-mul '(var pi var-pi) (math-solve-get-int '(var i var-i) (and (integerp (nth 2 lhs)) (math-abs (nth 2 lhs))))) (math-div (nth 2 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 lhs)) (math-abs (nth 2 lhs)))))) (nth 2 lhs))))) (math-normalize (list 'calcFunc-nroot rhs (nth 2 lhs)))) (and sign (math-oddp (nth 2 lhs)) (math-solve-sign sign (nth 2 lhs))))))))) (t nil))) (defun math-solve-prod (lsoln rsoln) (cond ((null lsoln) rsoln) ((null rsoln) lsoln) ((eq solve-full 'all) (cons 'vec (append (cdr lsoln) (cdr rsoln)))) (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". (defun math-solve-poly-funny-powers (sub-rhs) ; uses "t1", "t2" (setq t1 lhs) (let ((pp math-poly-neg-powers) fac) (while pp (setq fac (math-pow (car pp) (or math-poly-mult-powers 1)) t1 (math-mul t1 fac) rhs (math-mul rhs fac) pp (cdr pp)))) (if sub-rhs (setq t1 (math-sub t1 rhs))) (let ((math-poly-neg-powers nil)) (setq t2 (math-mul (or math-poly-mult-powers 1) (let ((calc-prefer-frac t)) (math-div 1 math-poly-frac-powers))) t1 (math-is-polynomial (math-simplify (calcFunc-expand t1)) 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 t1 (Math-zerop (car t1))) (setq t1 (cdr t1) count (1+ count))) (and t1 (let* ((degree (1- (length t1))) (scale degree)) (while (and (> scale 1) (= (car t3) 1)) (and (= (% degree scale) 0) (let ((p 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 t3 (cons scale (cdr t3)) t1 new-t1)))) (setq scale (1- scale))) (setq t3 (list (math-mul (car t3) t2) (math-mul count t2))) (<= (1- (length 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 (lhs solve-var degree sub-rhs) (let ((rhs (or sub-rhs 1)) t1 t2 t3) (setq t2 (math-polynomial-base lhs (function (lambda (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 b lhs)) (or (not (memq (car-safe b) '(+ -))) sub-rhs) (setq t3 '(1 0) t2 1 t1 (math-is-polynomial lhs 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 t1 (cons (math-sub (car t1) rhs) (cdr t1))) (math-solve-poly-funny-powers sub-rhs)) (math-solve-crunch-poly degree) (or (math-expr-contains b solve-var) (math-expr-contains (car t3) solve-var)))))))) (if t2 (list (math-pow t2 (car t3)) (cons 'vec t1) (if sub-rhs (math-pow t2 (nth 1 t3)) (math-div (math-pow t2 (nth 1 t3)) 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 ((solve-full nil) calc-next-why) (math-solve-cubic 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) (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 ((solve-var '(var DUMMY var-DUMMY)) (math-solve-sign nil) (math-solve-ranges nil) (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)) 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)) (solve-full (and 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 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 t3 neg 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 solve-var) (or (and (eq (car x) 'calcFunc-sqrt) (setq t2 2)) (and (eq (car x) '^) (or (and (memq (math-quarter-integer (nth 2 x)) '(1 2 3)) (setq t2 2)) (and (eq (car-safe (nth 2 x)) 'frac) (eq (nth 2 (nth 2 x)) 3) (setq t2 3)))) (and (memq (car x) '(* /)) (or (and (not (math-expr-contains (nth 1 x) solve-var)) (math-solve-find-root-in-prod (nth 2 x))) (and (not (math-expr-contains (nth 2 x) solve-var)) (math-solve-find-root-in-prod (nth 1 x)))))))) (defun math-solve-system (exprs solve-vars solve-full) (setq exprs (mapcar 'list (if (Math-vectorp exprs) (cdr exprs) (list exprs))) solve-vars (if (Math-vectorp solve-vars) (cdr solve-vars) (list solve-vars))) (or (let ((math-solve-simplifying nil)) (math-solve-system-rec exprs solve-vars nil)) (let ((math-solve-simplifying t)) (math-solve-system-rec exprs 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. (defun math-solve-system-rec (eqn-list var-list solns) (if var-list (let ((v var-list) (res nil)) ;; Try each variable in turn. (while (and v (let* ((vv (car v)) (e eqn-list) (elim (eq (car-safe vv) 'calcFunc-elim))) (if elim (setq vv (nth 1 vv))) ;; Try each equation in turn. (while (and e (let ((e2 (car e)) (eprev nil) res2) (setq res nil) ;; Try to solve for vv the list of equations e2. (while (and e2 (setq res2 (or (and (eq (car e2) eprev) res2) (math-solve-for (car e2) 0 vv solve-full)))) (setq eprev (car e2) res (cons (if (eq solve-full 'all) (cdr res2) (list res2)) res) e2 (cdr e2))) (if e2 (setq res nil) ;; Found a solution. Now try other variables. (setq res (nreverse res) 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 vv (apply 'append res)) s))))) (not res)))) (setq e (cdr e))) (not res))) (setq v (cdr v))) 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) solve-vars))))))) (if (eq 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 res)) (while x (setq accum (nconc accum (mapcar (function (lambda (r) (if math-solve-simplifying (math-simplify (math-expr-subst (car x) vv r)) (math-expr-subst (car x) vv r)))) (car res2))) x (cdr x) res2 (cdr res2))) accum)) (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))) (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 solve-full (if (and (calc-var-value 'var-GenCount) (Math-natnump var-GenCount) (not (eq 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 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 solve-full (if (and (calc-var-value 'var-GenCount) (Math-natnump var-GenCount) (not (eq 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 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 solve-var solve-full &optional sign) (if (math-expr-contains rhs solve-var) (math-solve-for (math-sub lhs rhs) 0 solve-var solve-full) (and (math-expr-contains lhs solve-var) (math-with-extra-prec 1 (let* ((math-poly-base-variable solve-var) (res (math-try-solve-for lhs rhs sign))) (if (and (eq solve-full 'all) (math-known-realp 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)))) ;;; arch-tag: f2932ec8-dd63-418b-a542-11a644b9d4c4 ;;; calcalg2.el ends here