mirror/emacs
1
0
Fork 0
mirror of https://git.savannah.gnu.org/git/emacs.git synced 2024-12-04 08:47:11 +00:00
Code Issues Releases Activity
d389648023
emacs/lisp/calc/calcalg2.el
Eli Zaretskii 136211a997 Initial import of Calc 2.02f.
2001-11-06 18:59:06 +00:00

3508 lines
109 KiB
EmacsLisp
Raw Blame History

;; Calculator for GNU Emacs, part II [calc-alg-2.el]
;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.
;; Written by Dave Gillespie, daveg@synaptics.com.
;; 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.
;; 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
(and (< 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))
(if (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)))))))
(setq math-integ-var '(var X ---))
(setq math-integ-var-2 '(var Y ---))
(setq math-integ-vars (list 'f math-integ-var math-integ-var-2))
(setq math-integ-var-list (list math-integ-var))
(setq 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))))))
)
(defun math-do-integral (expr)
(let (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 (if math-enable-subst
(math-integ-try-substitutions expr)
(math-integ-try-linear-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)
(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)))))))
(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)))))
)
(setq math-tabulate-initial nil)
(setq math-tabulate-function nil)
(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)))))
)
(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))
)
(setq math-disable-sums nil)
(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))))
)
;; 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))
)
(setq math-sum-int-pow-cache (list '(0 1)))
(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)))))
)
(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))
)
(setq math-disable-prods nil)
(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))))
)
;;; 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)))
)
(setq math-solve-ranges 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)
)
(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)))))
)
(setq math-symbolic-solve nil)
(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))))))))))
)
(setq math-int-coefs nil)
;;; 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" (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" (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)))
)
Reference in New Issue View Git Blame Copy Permalink