1 ;; Calculator for GNU Emacs, part II [calc-alg-2.el]
2 ;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.
3 ;; Written by Dave Gillespie, daveg@synaptics.com.
5 ;; This file is part of GNU Emacs.
7 ;; GNU Emacs is distributed in the hope that it will be useful,
8 ;; but WITHOUT ANY WARRANTY. No author or distributor
9 ;; accepts responsibility to anyone for the consequences of using it
10 ;; or for whether it serves any particular purpose or works at all,
11 ;; unless he says so in writing. Refer to the GNU Emacs General Public
12 ;; License for full details.
14 ;; Everyone is granted permission to copy, modify and redistribute
15 ;; GNU Emacs, but only under the conditions described in the
16 ;; GNU Emacs General Public License. A copy of this license is
17 ;; supposed to have been given to you along with GNU Emacs so you
18 ;; can know your rights and responsibilities. It should be in a
19 ;; file named COPYING. Among other things, the copyright notice
20 ;; and this notice must be preserved on all copies.
24 ;; This file is autoloaded from calc-ext.el.
29 (defun calc-Need-calc-alg-2 () nil)
32 (defun calc-derivative (var num)
33 (interactive "sDifferentiate with respect to: \np")
35 (and (< num 0) (error "Order of derivative must be positive"))
36 (let ((func (if (calc-is-hyperbolic) 'calcFunc-tderiv 'calcFunc-deriv))
38 (if (or (equal var "") (equal var "$"))
42 (setq var (math-read-expr var))
43 (if (eq (car-safe var) 'error)
44 (error "Bad format in expression: %s" (nth 1 var)))
47 (while (>= (setq num (1- num)) 0)
48 (setq expr (list func expr var)))
49 (calc-enter-result n "derv" expr)))
52 (defun calc-integral (var)
53 (interactive "sIntegration variable: ")
55 (if (or (equal var "") (equal var "$"))
56 (calc-enter-result 2 "intg" (list 'calcFunc-integ
59 (let ((var (math-read-expr var)))
60 (if (eq (car-safe var) 'error)
61 (error "Bad format in expression: %s" (nth 1 var)))
62 (calc-enter-result 1 "intg" (list 'calcFunc-integ
67 (defun calc-num-integral (&optional varname lowname highname)
68 (interactive "sIntegration variable: ")
69 (calc-tabular-command 'calcFunc-ninteg "Integration" "nint"
70 nil varname lowname highname)
73 (defun calc-summation (arg &optional varname lowname highname)
74 (interactive "P\nsSummation variable: ")
75 (calc-tabular-command 'calcFunc-sum "Summation" "sum"
76 arg varname lowname highname)
79 (defun calc-alt-summation (arg &optional varname lowname highname)
80 (interactive "P\nsSummation variable: ")
81 (calc-tabular-command 'calcFunc-asum "Summation" "asum"
82 arg varname lowname highname)
85 (defun calc-product (arg &optional varname lowname highname)
86 (interactive "P\nsIndex variable: ")
87 (calc-tabular-command 'calcFunc-prod "Index" "prod"
88 arg varname lowname highname)
91 (defun calc-tabulate (arg &optional varname lowname highname)
92 (interactive "P\nsIndex variable: ")
93 (calc-tabular-command 'calcFunc-table "Index" "tabl"
94 arg varname lowname highname)
97 (defun calc-tabular-command (func prompt prefix arg varname lowname highname)
99 (let (var (low nil) (high nil) (step nil) stepname stepnum (num 1) expr)
103 (if (or (equal varname "") (equal varname "$") (null varname))
104 (setq high (calc-top-n (+ stepnum 1))
105 low (calc-top-n (+ stepnum 2))
106 var (calc-top-n (+ stepnum 3))
108 (setq var (if (stringp varname) (math-read-expr varname) varname))
109 (if (eq (car-safe var) 'error)
110 (error "Bad format in expression: %s" (nth 1 var)))
112 (setq lowname (read-string (concat prompt " variable: " varname
114 (if (or (equal lowname "") (equal lowname "$"))
115 (setq high (calc-top-n (+ stepnum 1))
116 low (calc-top-n (+ stepnum 2))
118 (setq low (if (stringp lowname) (math-read-expr lowname) lowname))
119 (if (eq (car-safe low) 'error)
120 (error "Bad format in expression: %s" (nth 1 low)))
122 (setq highname (read-string (concat prompt " variable: " varname
125 (if (or (equal highname "") (equal highname "$"))
126 (setq high (calc-top-n (+ stepnum 1))
128 (setq high (if (stringp highname) (math-read-expr highname)
130 (if (eq (car-safe high) 'error)
131 (error "Bad format in expression: %s" (nth 1 high)))
134 (setq stepname (read-string (concat prompt " variable: "
139 (if (or (equal stepname "") (equal stepname "$"))
140 (setq step (calc-top-n 1)
142 (setq step (math-read-expr stepname))
143 (if (eq (car-safe step) 'error)
144 (error "Bad format in expression: %s"
148 (setq step (calc-top-n 1))
150 (setq step (prefix-numeric-value arg)))))
151 (setq expr (calc-top-n num))
152 (calc-enter-result num prefix (append (list func expr var low high)
153 (and step (list step))))))
156 (defun calc-solve-for (var)
157 (interactive "sVariable to solve for: ")
159 (let ((func (if (calc-is-inverse)
160 (if (calc-is-hyperbolic) 'calcFunc-ffinv 'calcFunc-finv)
161 (if (calc-is-hyperbolic) 'calcFunc-fsolve 'calcFunc-solve))))
162 (if (or (equal var "") (equal var "$"))
163 (calc-enter-result 2 "solv" (list func
166 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
167 (not (string-match "\\[" var)))
168 (math-read-expr (concat "[" var "]"))
169 (math-read-expr var))))
170 (if (eq (car-safe var) 'error)
171 (error "Bad format in expression: %s" (nth 1 var)))
172 (calc-enter-result 1 "solv" (list func
177 (defun calc-poly-roots (var)
178 (interactive "sVariable to solve for: ")
180 (if (or (equal var "") (equal var "$"))
181 (calc-enter-result 2 "prts" (list 'calcFunc-roots
184 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
185 (not (string-match "\\[" var)))
186 (math-read-expr (concat "[" var "]"))
187 (math-read-expr var))))
188 (if (eq (car-safe var) 'error)
189 (error "Bad format in expression: %s" (nth 1 var)))
190 (calc-enter-result 1 "prts" (list 'calcFunc-roots
195 (defun calc-taylor (var nterms)
196 (interactive "sTaylor expansion variable: \nNNumber of terms: ")
198 (let ((var (math-read-expr var)))
199 (if (eq (car-safe var) 'error)
200 (error "Bad format in expression: %s" (nth 1 var)))
201 (calc-enter-result 1 "tylr" (list 'calcFunc-taylor
204 (prefix-numeric-value nterms)))))
208 (defun math-derivative (expr) ; uses global values: deriv-var, deriv-total.
209 (cond ((equal expr deriv-var)
211 ((or (Math-scalarp expr)
212 (eq (car expr) 'sdev)
213 (and (eq (car expr) 'var)
214 (or (not deriv-total)
215 (math-const-var expr)
217 (math-setup-declarations)
218 (memq 'const (nth 1 (or (assq (nth 2 expr)
220 math-decls-all)))))))
223 (math-add (math-derivative (nth 1 expr))
224 (math-derivative (nth 2 expr))))
226 (math-sub (math-derivative (nth 1 expr))
227 (math-derivative (nth 2 expr))))
228 ((memq (car expr) '(calcFunc-eq calcFunc-neq calcFunc-lt
229 calcFunc-gt calcFunc-leq calcFunc-geq))
231 (math-derivative (nth 1 expr))
232 (math-derivative (nth 2 expr))))
233 ((eq (car expr) 'neg)
234 (math-neg (math-derivative (nth 1 expr))))
236 (math-add (math-mul (nth 2 expr)
237 (math-derivative (nth 1 expr)))
238 (math-mul (nth 1 expr)
239 (math-derivative (nth 2 expr)))))
241 (math-sub (math-div (math-derivative (nth 1 expr))
243 (math-div (math-mul (nth 1 expr)
244 (math-derivative (nth 2 expr)))
245 (math-sqr (nth 2 expr)))))
247 (let ((du (math-derivative (nth 1 expr)))
248 (dv (math-derivative (nth 2 expr))))
250 (setq du (math-mul (nth 2 expr)
251 (math-mul (math-normalize
254 (math-add (nth 2 expr) -1)))
257 (setq dv (math-mul (math-normalize
258 (list 'calcFunc-ln (nth 1 expr)))
259 (math-mul expr dv))))
262 (math-derivative (nth 1 expr))) ; a reasonable definition
263 ((eq (car expr) 'vec)
264 (math-map-vec 'math-derivative expr))
265 ((and (memq (car expr) '(calcFunc-conj calcFunc-re calcFunc-im))
267 (list (car expr) (math-derivative (nth 1 expr))))
268 ((and (memq (car expr) '(calcFunc-subscr calcFunc-mrow calcFunc-mcol))
270 (let ((d (math-derivative (nth 1 expr))))
272 0 ; assume x and x_1 are independent vars
273 (list (car expr) d (nth 2 expr)))))
274 (t (or (and (symbolp (car expr))
275 (if (= (length expr) 2)
276 (let ((handler (get (car expr) 'math-derivative)))
278 (let ((deriv (math-derivative (nth 1 expr))))
279 (if (Math-zerop deriv)
281 (math-mul (funcall handler (nth 1 expr))
283 (let ((handler (get (car expr) 'math-derivative-n)))
285 (funcall handler expr)))))
286 (and (not (eq deriv-symb 'pre-expand))
287 (let ((exp (math-expand-formula expr)))
289 (or (let ((deriv-symb 'pre-expand))
290 (catch 'math-deriv (math-derivative expr)))
291 (math-derivative exp)))))
292 (if (or (Math-objvecp expr)
294 (not (symbolp (car expr))))
296 (throw 'math-deriv nil)
297 (list (if deriv-total 'calcFunc-tderiv 'calcFunc-deriv)
304 (while (setq arg (cdr arg))
305 (or (Math-zerop (setq derv (math-derivative (car arg))))
306 (let ((func (intern (concat (symbol-name (car expr))
311 (prop (cond ((= (length expr) 2)
320 'math-derivative-5))))
326 (let ((handler (get func prop)))
327 (or (and prop handler
328 (apply handler (cdr expr)))
332 (throw 'math-deriv nil)
333 (cons func (cdr expr))))))))))
338 (defun calcFunc-deriv (expr deriv-var &optional deriv-value deriv-symb)
339 (let* ((deriv-total nil)
340 (res (catch 'math-deriv (math-derivative expr))))
341 (or (eq (car-safe res) 'calcFunc-deriv)
343 (setq res (math-normalize res)))
346 (math-expr-subst res deriv-var deriv-value)
350 (defun calcFunc-tderiv (expr deriv-var &optional deriv-value deriv-symb)
351 (math-setup-declarations)
352 (let* ((deriv-total t)
353 (res (catch 'math-deriv (math-derivative expr))))
354 (or (eq (car-safe res) 'calcFunc-tderiv)
356 (setq res (math-normalize res)))
359 (math-expr-subst res deriv-var deriv-value)
363 (put 'calcFunc-inv\' 'math-derivative-1
364 (function (lambda (u) (math-neg (math-div 1 (math-sqr u))))))
366 (put 'calcFunc-sqrt\' 'math-derivative-1
367 (function (lambda (u) (math-div 1 (math-mul 2 (list 'calcFunc-sqrt u))))))
369 (put 'calcFunc-deg\' 'math-derivative-1
370 (function (lambda (u) (math-div-float '(float 18 1) (math-pi)))))
372 (put 'calcFunc-rad\' 'math-derivative-1
373 (function (lambda (u) (math-pi-over-180))))
375 (put 'calcFunc-ln\' 'math-derivative-1
376 (function (lambda (u) (math-div 1 u))))
378 (put 'calcFunc-log10\' 'math-derivative-1
379 (function (lambda (u)
380 (math-div (math-div 1 (math-normalize '(calcFunc-ln 10)))
383 (put 'calcFunc-lnp1\' 'math-derivative-1
384 (function (lambda (u) (math-div 1 (math-add u 1)))))
386 (put 'calcFunc-log\' 'math-derivative-2
387 (function (lambda (x b)
388 (and (not (Math-zerop b))
389 (let ((lnv (math-normalize
390 (list 'calcFunc-ln b))))
391 (math-div 1 (math-mul lnv x)))))))
393 (put 'calcFunc-log\'2 'math-derivative-2
394 (function (lambda (x b)
395 (let ((lnv (list 'calcFunc-ln b)))
396 (math-neg (math-div (list 'calcFunc-log x b)
397 (math-mul lnv b)))))))
399 (put 'calcFunc-exp\' 'math-derivative-1
400 (function (lambda (u) (math-normalize (list 'calcFunc-exp u)))))
402 (put 'calcFunc-expm1\' 'math-derivative-1
403 (function (lambda (u) (math-normalize (list 'calcFunc-expm1 u)))))
405 (put 'calcFunc-sin\' 'math-derivative-1
406 (function (lambda (u) (math-to-radians-2 (math-normalize
407 (list 'calcFunc-cos u))))))
409 (put 'calcFunc-cos\' 'math-derivative-1
410 (function (lambda (u) (math-neg (math-to-radians-2
412 (list 'calcFunc-sin u)))))))
414 (put 'calcFunc-tan\' 'math-derivative-1
415 (function (lambda (u) (math-to-radians-2
416 (math-div 1 (math-sqr
418 (list 'calcFunc-cos u))))))))
420 (put 'calcFunc-arcsin\' 'math-derivative-1
421 (function (lambda (u)
423 (math-div 1 (math-normalize
425 (math-sub 1 (math-sqr u)))))))))
427 (put 'calcFunc-arccos\' 'math-derivative-1
428 (function (lambda (u)
430 (math-div -1 (math-normalize
432 (math-sub 1 (math-sqr u)))))))))
434 (put 'calcFunc-arctan\' 'math-derivative-1
435 (function (lambda (u) (math-from-radians-2
436 (math-div 1 (math-add 1 (math-sqr u)))))))
438 (put 'calcFunc-sinh\' 'math-derivative-1
439 (function (lambda (u) (math-normalize (list 'calcFunc-cosh u)))))
441 (put 'calcFunc-cosh\' 'math-derivative-1
442 (function (lambda (u) (math-normalize (list 'calcFunc-sinh u)))))
444 (put 'calcFunc-tanh\' 'math-derivative-1
445 (function (lambda (u) (math-div 1 (math-sqr
447 (list 'calcFunc-cosh u)))))))
449 (put 'calcFunc-arcsinh\' 'math-derivative-1
450 (function (lambda (u)
451 (math-div 1 (math-normalize
453 (math-add (math-sqr u) 1)))))))
455 (put 'calcFunc-arccosh\' 'math-derivative-1
456 (function (lambda (u)
457 (math-div 1 (math-normalize
459 (math-add (math-sqr u) -1)))))))
461 (put 'calcFunc-arctanh\' 'math-derivative-1
462 (function (lambda (u) (math-div 1 (math-sub 1 (math-sqr u))))))
464 (put 'calcFunc-bern\'2 'math-derivative-2
465 (function (lambda (n x)
466 (math-mul n (list 'calcFunc-bern (math-add n -1) x)))))
468 (put 'calcFunc-euler\'2 'math-derivative-2
469 (function (lambda (n x)
470 (math-mul n (list 'calcFunc-euler (math-add n -1) x)))))
472 (put 'calcFunc-gammag\'2 'math-derivative-2
473 (function (lambda (a x) (math-deriv-gamma a x 1))))
475 (put 'calcFunc-gammaG\'2 'math-derivative-2
476 (function (lambda (a x) (math-deriv-gamma a x -1))))
478 (put 'calcFunc-gammaP\'2 'math-derivative-2
479 (function (lambda (a x) (math-deriv-gamma a x
482 (list 'calcFunc-gamma
485 (put 'calcFunc-gammaQ\'2 'math-derivative-2
486 (function (lambda (a x) (math-deriv-gamma a x
489 (list 'calcFunc-gamma
492 (defun math-deriv-gamma (a x scale)
494 (math-mul (math-pow x (math-add a -1))
495 (list 'calcFunc-exp (math-neg x))))
498 (put 'calcFunc-betaB\' 'math-derivative-3
499 (function (lambda (x a b) (math-deriv-beta x a b 1))))
501 (put 'calcFunc-betaI\' 'math-derivative-3
502 (function (lambda (x a b) (math-deriv-beta x a b
504 1 (list 'calcFunc-beta
507 (defun math-deriv-beta (x a b scale)
508 (math-mul (math-mul (math-pow x (math-add a -1))
509 (math-pow (math-sub 1 x) (math-add b -1)))
513 (put 'calcFunc-erf\' 'math-derivative-1
514 (function (lambda (x) (math-div 2
515 (math-mul (list 'calcFunc-exp
517 (if calc-symbolic-mode
522 (put 'calcFunc-erfc\' 'math-derivative-1
523 (function (lambda (x) (math-div -2
524 (math-mul (list 'calcFunc-exp
526 (if calc-symbolic-mode
531 (put 'calcFunc-besJ\'2 'math-derivative-2
532 (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besJ
540 (put 'calcFunc-besY\'2 'math-derivative-2
541 (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besY
549 (put 'calcFunc-sum 'math-derivative-n
552 (if (math-expr-contains (cons 'vec (cdr (cdr expr))) deriv-var)
553 (throw 'math-deriv nil)
555 (cons (math-derivative (nth 1 expr))
556 (cdr (cdr expr))))))))
558 (put 'calcFunc-prod 'math-derivative-n
561 (if (math-expr-contains (cons 'vec (cdr (cdr expr))) deriv-var)
562 (throw 'math-deriv nil)
565 (cons (math-div (math-derivative (nth 1 expr))
567 (cdr (cdr expr)))))))))
569 (put 'calcFunc-integ 'math-derivative-n
572 (if (= (length expr) 3)
573 (if (equal (nth 2 expr) deriv-var)
576 (list 'calcFunc-integ
577 (math-derivative (nth 1 expr))
579 (if (= (length expr) 5)
580 (let ((lower (math-expr-subst (nth 1 expr) (nth 2 expr)
582 (upper (math-expr-subst (nth 1 expr) (nth 2 expr)
584 (math-add (math-sub (math-mul upper
585 (math-derivative (nth 4 expr)))
587 (math-derivative (nth 3 expr))))
588 (if (equal (nth 2 expr) deriv-var)
591 (list 'calcFunc-integ
592 (math-derivative (nth 1 expr)) (nth 2 expr)
593 (nth 3 expr) (nth 4 expr)))))))))))
595 (put 'calcFunc-if 'math-derivative-n
598 (and (= (length expr) 4)
599 (list 'calcFunc-if (nth 1 expr)
600 (math-derivative (nth 2 expr))
601 (math-derivative (nth 3 expr)))))))
603 (put 'calcFunc-subscr 'math-derivative-n
606 (and (= (length expr) 3)
607 (list 'calcFunc-subscr (nth 1 expr)
608 (math-derivative (nth 2 expr)))))))
614 (setq math-integ-var '(var X ---))
615 (setq math-integ-var-2 '(var Y ---))
616 (setq math-integ-vars (list 'f math-integ-var math-integ-var-2))
617 (setq math-integ-var-list (list math-integ-var))
618 (setq math-integ-var-list-list (list math-integ-var-list))
620 (defmacro math-tracing-integral (&rest parts)
623 (list 'save-excursion
624 '(set-buffer trace-buffer)
625 '(goto-char (point-max))
628 '(insert (make-string (- math-integral-limit
629 math-integ-level) 32)
630 (format "%2d " math-integ-depth)
631 (make-string math-integ-level 32)))
632 ;;(list 'condition-case 'err
634 ;; '(error (insert (prin1-to-string err))))
638 ;;; The following wrapper caches results and avoids infinite recursion.
639 ;;; Each cache entry is: ( A B ) Integral of A is B;
640 ;;; ( A N ) Integral of A failed at level N;
641 ;;; ( A busy ) Currently working on integral of A;
642 ;;; ( A parts ) Currently working, integ-by-parts;
643 ;;; ( A parts2 ) Currently working, integ-by-parts;
644 ;;; ( A cancelled ) Ignore this cache entry;
645 ;;; ( A [B] ) Same result as for cur-record = B.
646 (defun math-integral (expr &optional simplify same-as-above)
647 (let* ((simp cur-record)
648 (cur-record (assoc expr math-integral-cache))
649 (math-integ-depth (1+ math-integ-depth))
651 (math-tracing-integral "Integrating "
652 (math-format-value expr 1000)
656 (math-tracing-integral "Found "
657 (math-format-value (nth 1 cur-record) 1000))
658 (and (consp (nth 1 cur-record))
659 (math-replace-integral-parts cur-record))
660 (math-tracing-integral " => "
661 (math-format-value (nth 1 cur-record) 1000)
664 (not (eq (nth 1 cur-record) 'cancelled))
665 (or (not (integerp (nth 1 cur-record)))
666 (>= (nth 1 cur-record) math-integ-level)))
667 (and (math-integral-contains-parts expr)
673 (let (math-integ-msg)
674 (if (eq calc-display-working-message 'lots)
676 (calc-set-command-flag 'clear-message)
677 (setq math-integ-msg (format
678 "Working... Integrating %s"
679 (math-format-flat-expr expr 0)))
680 (message math-integ-msg)))
682 (setcar (cdr cur-record)
683 (if same-as-above (vector simp) 'busy))
685 (list expr (if same-as-above (vector simp) 'busy))
686 math-integral-cache (cons cur-record
687 math-integral-cache)))
688 (if (eq simplify 'yes)
690 (math-tracing-integral "Simplifying...")
691 (setq simp (math-simplify expr))
692 (setq val (if (equal simp expr)
694 (math-tracing-integral " no change\n")
695 (math-do-integral expr))
696 (math-tracing-integral " simplified\n")
697 (math-integral simp 'no t))))
698 (or (setq val (math-do-integral expr))
700 (let ((simp (math-simplify expr)))
701 (or (equal simp expr)
703 (math-tracing-integral "Trying again after "
704 "simplification...\n")
705 (setq val (math-integral simp 'no t))))))))
706 (if (eq calc-display-working-message 'lots)
707 (message math-integ-msg)))
708 (setcar (cdr cur-record) (or val
709 (if (or math-enable-subst
710 (not math-any-substs))
713 (setq val cur-record)
714 (while (vectorp (nth 1 val))
715 (setq val (aref (nth 1 val) 0)))
716 (setq val (if (memq (nth 1 val) '(parts parts2))
718 (setcar (cdr val) 'parts2)
719 (list 'var 'PARTS val))
720 (and (consp (nth 1 val))
722 (math-tracing-integral "Integral of "
723 (math-format-value expr 1000)
725 (math-format-value val 1000)
729 (defvar math-integral-cache nil)
730 (defvar math-integral-cache-state nil)
732 (defun math-integral-contains-parts (expr)
733 (if (Math-primp expr)
734 (and (eq (car-safe expr) 'var)
735 (eq (nth 1 expr) 'PARTS)
736 (listp (nth 2 expr)))
737 (while (and (setq expr (cdr expr))
738 (not (math-integral-contains-parts (car expr)))))
742 (defun math-replace-integral-parts (expr)
743 (or (Math-primp expr)
744 (while (setq expr (cdr expr))
745 (and (consp (car expr))
746 (if (eq (car (car expr)) 'var)
747 (and (eq (nth 1 (car expr)) 'PARTS)
748 (consp (nth 2 (car expr)))
749 (if (listp (nth 1 (nth 2 (car expr))))
751 (setcar expr (nth 1 (nth 2 (car expr))))
752 (math-replace-integral-parts (cons 'foo expr)))
753 (setcar (cdr cur-record) 'cancelled)))
754 (math-replace-integral-parts (car expr))))))
757 (defun math-do-integral (expr)
759 (or (cond ((not (math-expr-contains expr math-integ-var))
760 (math-mul expr math-integ-var))
761 ((equal expr math-integ-var)
762 (math-div (math-sqr expr) 2))
764 (and (setq t1 (math-integral (nth 1 expr)))
765 (setq t2 (math-integral (nth 2 expr)))
768 (and (setq t1 (math-integral (nth 1 expr)))
769 (setq t2 (math-integral (nth 2 expr)))
771 ((eq (car expr) 'neg)
772 (and (setq t1 (math-integral (nth 1 expr)))
775 (cond ((not (math-expr-contains (nth 1 expr) math-integ-var))
776 (and (setq t1 (math-integral (nth 2 expr)))
777 (math-mul (nth 1 expr) t1)))
778 ((not (math-expr-contains (nth 2 expr) math-integ-var))
779 (and (setq t1 (math-integral (nth 1 expr)))
780 (math-mul t1 (nth 2 expr))))
781 ((memq (car-safe (nth 1 expr)) '(+ -))
782 (math-integral (list (car (nth 1 expr))
783 (math-mul (nth 1 (nth 1 expr))
785 (math-mul (nth 2 (nth 1 expr))
788 ((memq (car-safe (nth 2 expr)) '(+ -))
789 (math-integral (list (car (nth 2 expr))
790 (math-mul (nth 1 (nth 2 expr))
792 (math-mul (nth 2 (nth 2 expr))
796 (cond ((and (not (math-expr-contains (nth 1 expr)
798 (not (math-equal-int (nth 1 expr) 1)))
799 (and (setq t1 (math-integral (math-div 1 (nth 2 expr))))
800 (math-mul (nth 1 expr) t1)))
801 ((not (math-expr-contains (nth 2 expr) math-integ-var))
802 (and (setq t1 (math-integral (nth 1 expr)))
803 (math-div t1 (nth 2 expr))))
804 ((and (eq (car-safe (nth 1 expr)) '*)
805 (not (math-expr-contains (nth 1 (nth 1 expr))
807 (and (setq t1 (math-integral
808 (math-div (nth 2 (nth 1 expr))
810 (math-mul t1 (nth 1 (nth 1 expr)))))
811 ((and (eq (car-safe (nth 1 expr)) '*)
812 (not (math-expr-contains (nth 2 (nth 1 expr))
814 (and (setq t1 (math-integral
815 (math-div (nth 1 (nth 1 expr))
817 (math-mul t1 (nth 2 (nth 1 expr)))))
818 ((and (eq (car-safe (nth 2 expr)) '*)
819 (not (math-expr-contains (nth 1 (nth 2 expr))
821 (and (setq t1 (math-integral
822 (math-div (nth 1 expr)
823 (nth 2 (nth 2 expr)))))
824 (math-div t1 (nth 1 (nth 2 expr)))))
825 ((and (eq (car-safe (nth 2 expr)) '*)
826 (not (math-expr-contains (nth 2 (nth 2 expr))
828 (and (setq t1 (math-integral
829 (math-div (nth 1 expr)
830 (nth 1 (nth 2 expr)))))
831 (math-div t1 (nth 2 (nth 2 expr)))))
832 ((eq (car-safe (nth 2 expr)) 'calcFunc-exp)
834 (math-mul (nth 1 expr)
836 (math-neg (nth 1 (nth 2 expr)))))))))
838 (cond ((not (math-expr-contains (nth 1 expr) math-integ-var))
839 (or (and (setq t1 (math-is-polynomial (nth 2 expr)
848 (math-mul (nth 2 expr)
853 ((not (math-expr-contains (nth 2 expr) math-integ-var))
854 (if (and (integerp (nth 2 expr)) (< (nth 2 expr) 0))
856 (list '/ 1 (math-pow (nth 1 expr) (- (nth 2 expr))))
858 (or (and (setq t1 (math-is-polynomial (nth 1 expr)
861 (setq t2 (math-add (nth 2 expr) 1))
862 (math-div (math-pow (nth 1 expr) t2)
863 (math-mul t2 (nth 1 t1))))
864 (and (Math-negp (nth 2 expr))
867 (math-pow (nth 1 expr)
873 ;; Integral of a polynomial.
874 (and (setq t1 (math-is-polynomial expr math-integ-var 20))
878 (if (setq accum (math-add accum
879 (math-div (math-mul (car t1)
888 ;; Try looking it up!
889 (cond ((= (length expr) 2)
890 (and (symbolp (car expr))
891 (setq t1 (get (car expr) 'math-integral))
894 (not (setq t2 (funcall (car t1)
897 (and t2 (math-normalize t2)))))
899 (and (symbolp (car expr))
900 (setq t1 (get (car expr) 'math-integral-2))
903 (not (setq t2 (funcall (car t1)
907 (and t2 (math-normalize t2))))))
909 ;; Integral of a rational function.
910 (and (math-ratpoly-p expr math-integ-var)
911 (setq t1 (calcFunc-apart expr math-integ-var))
912 (not (equal t1 expr))
915 ;; Try user-defined integration rules.
917 (let ((math-old-integ (symbol-function 'calcFunc-integ))
918 (input (list 'calcFunc-integtry expr math-integ-var))
922 (fset 'calcFunc-integ 'math-sub-integration)
923 (setq res (math-rewrite input
924 '(var IntegRules var-IntegRules)
926 (fset 'calcFunc-integ math-old-integ)
927 (and (not (equal res input))
928 (if (setq part (math-expr-calls
929 res '(calcFunc-integsubst)))
930 (and (memq (length part) '(3 4 5))
938 (math-integrate-by-substitution
941 (list 'calcFunc-integfailed
944 (if (not (math-expr-calls res
946 calcFunc-integfailed)))
948 (fset 'calcFunc-integ math-old-integ))))
950 ;; See if the function is a symbolic derivative.
951 (and (string-match "'" (symbol-name (car expr)))
952 (let ((name (symbol-name (car expr)))
953 (p expr) (n 0) (which nil) (bad nil))
954 (while (setq n (1+ n) p (cdr p))
955 (if (equal (car p) math-integ-var)
956 (if which (setq bad t) (setq which n))
957 (if (math-expr-contains (car p) math-integ-var)
960 (let ((prime (if (= which 1) "'" (format "'%d" which))))
961 (and (string-match (concat prime "\\('['0-9]*\\|$\\)")
965 (substring name 0 (match-beginning 0))
966 (substring name (+ (match-beginning 0)
970 ;; Try transformation methods (parts, substitutions).
971 (and (> math-integ-level 0)
972 (math-do-integral-methods expr))
974 ;; Try expanding the function's definition.
975 (let ((res (math-expand-formula expr)))
977 (math-integral res)))))
980 (defun math-sub-integration (expr &rest rest)
981 (or (if (or (not rest)
982 (and (< math-integ-level math-integral-limit)
983 (eq (car rest) math-integ-var)))
985 (let ((res (apply math-old-integ expr rest)))
986 (and (or (= math-integ-level math-integral-limit)
987 (not (math-expr-calls res 'calcFunc-integ)))
989 (list 'calcFunc-integfailed expr))
992 (defun math-do-integral-methods (expr)
993 (let ((so-far math-integ-var-list-list)
996 ;; Integration by substitution, for various likely sub-expressions.
997 ;; (In first pass, we look only for sub-exprs that are linear in X.)
998 (or (if math-enable-subst
999 (math-integ-try-substitutions expr)
1000 (math-integ-try-linear-substitutions expr))
1002 ;; If function has sines and cosines, try tan(x/2) substitution.
1003 (and (let ((p (setq rat-in (math-expr-rational-in expr))))
1005 (memq (car (car p)) '(calcFunc-sin
1008 (equal (nth 1 (car p)) math-integ-var))
1011 (or (and (math-integ-parts-easy expr)
1012 (math-integ-try-parts expr t))
1013 (math-integrate-by-good-substitution
1014 expr (list 'calcFunc-tan (math-div math-integ-var 2)))))
1016 ;; If function has sinh and cosh, try tanh(x/2) substitution.
1017 (and (let ((p rat-in))
1019 (memq (car (car p)) '(calcFunc-sinh
1023 (equal (nth 1 (car p)) math-integ-var))
1026 (or (and (math-integ-parts-easy expr)
1027 (math-integ-try-parts expr t))
1028 (math-integrate-by-good-substitution
1029 expr (list 'calcFunc-tanh (math-div math-integ-var 2)))))
1031 ;; If function has square roots, try sin, tan, or sec substitution.
1032 (and (let ((p rat-in))
1035 (or (equal (car p) math-integ-var)
1036 (and (eq (car (car p)) 'calcFunc-sqrt)
1037 (setq t1 (math-is-polynomial
1038 (nth 1 (setq t2 (car p)))
1039 math-integ-var 2)))))
1043 (if (math-guess-if-neg (nth 2 t1))
1044 (let* ((c (math-sqrt (math-neg (nth 2 t1))))
1045 (d (math-div (nth 1 t1) (math-mul -2 c)))
1046 (a (math-sqrt (math-add (car t1) (math-sqr d)))))
1047 (math-integrate-by-good-substitution
1048 expr (list 'calcFunc-arcsin
1050 (math-add (math-mul c math-integ-var) d)
1052 (let* ((c (math-sqrt (nth 2 t1)))
1053 (d (math-div (nth 1 t1) (math-mul 2 c)))
1054 (aa (math-sub (car t1) (math-sqr d))))
1055 (if (and nil (not (and (eq d 0) (eq c 1))))
1056 (math-integrate-by-good-substitution
1057 expr (math-add (math-mul c math-integ-var) d))
1058 (if (math-guess-if-neg aa)
1059 (math-integrate-by-good-substitution
1060 expr (list 'calcFunc-arccosh
1062 (math-add (math-mul c math-integ-var)
1064 (math-sqrt (math-neg aa)))))
1065 (math-integrate-by-good-substitution
1066 expr (list 'calcFunc-arcsinh
1068 (math-add (math-mul c math-integ-var)
1070 (math-sqrt aa))))))))
1071 (math-integrate-by-good-substitution expr t2)) )
1073 ;; Try integration by parts.
1074 (math-integ-try-parts expr)
1080 (defun math-integ-parts-easy (expr)
1081 (cond ((Math-primp expr) t)
1082 ((memq (car expr) '(+ - *))
1083 (and (math-integ-parts-easy (nth 1 expr))
1084 (math-integ-parts-easy (nth 2 expr))))
1086 (and (math-integ-parts-easy (nth 1 expr))
1087 (math-atomic-factorp (nth 2 expr))))
1089 (and (natnump (nth 2 expr))
1090 (math-integ-parts-easy (nth 1 expr))))
1091 ((eq (car expr) 'neg)
1092 (math-integ-parts-easy (nth 1 expr)))
1096 (defun math-integ-try-parts (expr &optional math-good-parts)
1097 ;; Integration by parts:
1098 ;; integ(f(x) g(x),x) = f(x) h(x) - integ(h(x) f'(x),x)
1099 ;; where h(x) = integ(g(x),x).
1100 (or (let ((exp (calcFunc-expand expr)))
1101 (and (not (equal exp expr))
1102 (math-integral exp)))
1103 (and (eq (car expr) '*)
1104 (let ((first-bad (or (math-polynomial-p (nth 1 expr)
1106 (equal (nth 2 expr) math-prev-parts-v))))
1107 (or (and first-bad ; so try this one first
1108 (math-integrate-by-parts (nth 1 expr) (nth 2 expr)))
1109 (math-integrate-by-parts (nth 2 expr) (nth 1 expr))
1110 (and (not first-bad)
1111 (math-integrate-by-parts (nth 1 expr) (nth 2 expr))))))
1112 (and (eq (car expr) '/)
1113 (math-expr-contains (nth 1 expr) math-integ-var)
1114 (let ((recip (math-div 1 (nth 2 expr))))
1115 (or (math-integrate-by-parts (nth 1 expr) recip)
1116 (math-integrate-by-parts recip (nth 1 expr)))))
1117 (and (eq (car expr) '^)
1118 (math-integrate-by-parts (math-pow (nth 1 expr)
1119 (math-sub (nth 2 expr) 1))
1123 (defun math-integrate-by-parts (u vprime)
1124 (let ((math-integ-level (if (or math-good-parts
1125 (math-polynomial-p u math-integ-var))
1127 (1- math-integ-level)))
1128 (math-doing-parts t)
1130 (and (>= math-integ-level 0)
1133 (setcar (cdr cur-record) 'parts)
1134 (math-tracing-integral "Integrating by parts, u = "
1135 (math-format-value u 1000)
1137 (math-format-value vprime 1000)
1139 (and (setq v (math-integral vprime))
1140 (setq temp (calcFunc-deriv u math-integ-var nil t))
1141 (setq temp (let ((math-prev-parts-v v))
1142 (math-integral (math-mul v temp) 'yes)))
1143 (setq temp (math-sub (math-mul u v) temp))
1144 (if (eq (nth 1 cur-record) 'parts)
1145 (calcFunc-expand temp)
1146 (setq v (list 'var 'PARTS cur-record)
1147 var-thing (list 'vec (math-sub v temp) v)
1148 temp (let (calc-next-why)
1149 (math-solve-for (math-sub v temp) 0 v nil)))
1150 (and temp (not (integerp temp))
1151 (math-simplify-extended temp)))))
1152 (setcar (cdr cur-record) 'busy))))
1155 ;;; This tries two different formulations, hoping the algebraic simplifier
1156 ;;; will be strong enough to handle at least one.
1157 (defun math-integrate-by-substitution (expr u &optional user uinv uinvprime)
1158 (and (> math-integ-level 0)
1159 (let ((math-integ-level (max (- math-integ-level 2) 0)))
1160 (math-integrate-by-good-substitution expr u user uinv uinvprime)))
1163 (defun math-integrate-by-good-substitution (expr u &optional user
1165 (let ((math-living-dangerously t)
1167 (and (setq uinv (if uinv
1168 (math-expr-subst uinv math-integ-var
1170 (let (calc-next-why)
1173 math-integ-var nil))))
1175 (math-tracing-integral "Integrating by substitution, u = "
1176 (math-format-value u 1000)
1178 (or (and (setq deriv (calcFunc-deriv u
1181 (setq temp (math-integral (math-expr-subst
1184 (math-div expr deriv)
1192 (and (setq deriv (or uinvprime
1193 (calcFunc-deriv uinv
1197 (setq temp (math-integral (math-mul
1210 (math-simplify-extended
1211 (math-expr-subst temp math-integ-var u))))
1214 ;;; Look for substitutions of the form u = a x + b.
1215 (defun math-integ-try-linear-substitutions (sub-expr)
1216 (and (not (Math-primp sub-expr))
1217 (or (and (not (memq (car sub-expr) '(+ - * / neg)))
1218 (not (and (eq (car sub-expr) '^)
1219 (integerp (nth 2 sub-expr))))
1220 (math-expr-contains sub-expr math-integ-var)
1222 (while (and (setq sub-expr (cdr sub-expr))
1223 (or (not (math-linear-in (car sub-expr)
1225 (assoc (car sub-expr) so-far)
1227 (setq so-far (cons (list (car sub-expr))
1230 (math-integrate-by-substitution
1231 expr (car sub-expr))))))))
1234 (while (and (setq sub-expr (cdr sub-expr))
1235 (not (setq res (math-integ-try-linear-substitutions
1240 ;;; Recursively try different substitutions based on various sub-expressions.
1241 (defun math-integ-try-substitutions (sub-expr &optional allow-rat)
1242 (and (not (Math-primp sub-expr))
1243 (not (assoc sub-expr so-far))
1244 (math-expr-contains sub-expr math-integ-var)
1245 (or (and (if (and (not (memq (car sub-expr) '(+ - * / neg)))
1246 (not (and (eq (car sub-expr) '^)
1247 (integerp (nth 2 sub-expr)))))
1249 (prog1 allow-rat (setq allow-rat nil)))
1250 (not (eq sub-expr expr))
1251 (or (math-integrate-by-substitution expr sub-expr)
1252 (and (eq (car sub-expr) '^)
1253 (integerp (nth 2 sub-expr))
1254 (< (nth 2 sub-expr) 0)
1255 (math-integ-try-substitutions
1256 (math-pow (nth 1 sub-expr) (- (nth 2 sub-expr)))
1259 (setq so-far (cons (list sub-expr) so-far))
1260 (while (and (setq sub-expr (cdr sub-expr))
1261 (not (setq res (math-integ-try-substitutions
1262 (car sub-expr) allow-rat)))))
1266 (defun math-expr-rational-in (expr)
1268 (math-expr-rational-in-rec expr)
1269 (mapcar 'car parts))
1272 (defun math-expr-rational-in-rec (expr)
1273 (cond ((Math-primp expr)
1274 (and (equal expr math-integ-var)
1275 (not (assoc expr parts))
1276 (setq parts (cons (list expr) parts))))
1277 ((or (memq (car expr) '(+ - * / neg))
1278 (and (eq (car expr) '^) (integerp (nth 2 expr))))
1279 (math-expr-rational-in-rec (nth 1 expr))
1280 (and (nth 2 expr) (math-expr-rational-in-rec (nth 2 expr))))
1281 ((and (eq (car expr) '^)
1282 (eq (math-quarter-integer (nth 2 expr)) 2))
1283 (math-expr-rational-in-rec (list 'calcFunc-sqrt (nth 1 expr))))
1285 (and (not (assoc expr parts))
1286 (math-expr-contains expr math-integ-var)
1287 (setq parts (cons (list expr) parts)))))
1290 (defun math-expr-calls (expr funcs &optional arg-contains)
1292 (if (or (memq (car expr) funcs)
1293 (and (eq (car expr) '^) (eq (car funcs) 'calcFunc-sqrt)
1294 (eq (math-quarter-integer (nth 2 expr)) 2)))
1295 (and (or (not arg-contains)
1296 (math-expr-contains expr arg-contains))
1298 (and (not (Math-primp expr))
1300 (while (and (setq expr (cdr expr))
1301 (not (setq res (math-expr-calls
1302 (car expr) funcs arg-contains)))))
1306 (defun math-fix-const-terms (expr except-vars)
1307 (cond ((not (math-expr-depends expr except-vars)) 0)
1308 ((Math-primp expr) expr)
1310 (math-add (math-fix-const-terms (nth 1 expr) except-vars)
1311 (math-fix-const-terms (nth 2 expr) except-vars)))
1313 (math-sub (math-fix-const-terms (nth 1 expr) except-vars)
1314 (math-fix-const-terms (nth 2 expr) except-vars)))
1318 ;; Command for debugging the Calculator's symbolic integrator.
1319 (defun calc-dump-integral-cache (&optional arg)
1321 (let ((buf (current-buffer)))
1323 (let ((p math-integral-cache)
1325 (display-buffer (get-buffer-create "*Integral Cache*"))
1326 (set-buffer (get-buffer "*Integral Cache*"))
1329 (setq cur-record (car p))
1330 (or arg (math-replace-integral-parts cur-record))
1331 (insert (math-format-flat-expr (car cur-record) 0)
1333 (if (symbolp (nth 1 cur-record))
1334 (concat "(" (symbol-name (nth 1 cur-record)) ")")
1335 (math-format-flat-expr (nth 1 cur-record) 0))
1338 (goto-char (point-min)))
1342 (defun math-try-integral (expr)
1343 (let ((math-integ-level math-integral-limit)
1344 (math-integ-depth 0)
1345 (math-integ-msg "Working...done")
1346 (cur-record nil) ; a technicality
1347 (math-integrating t)
1348 (calc-prefer-frac t)
1349 (calc-symbolic-mode t)
1350 (has-rules (calc-has-rules 'var-IntegRules)))
1351 (or (math-integral expr 'yes)
1352 (and math-any-substs
1353 (setq math-enable-subst t)
1354 (math-integral expr 'yes))
1355 (and (> math-max-integral-limit math-integral-limit)
1356 (setq math-integral-limit math-max-integral-limit
1357 math-integ-level math-integral-limit)
1358 (math-integral expr 'yes))))
1361 (defun calcFunc-integ (expr var &optional low high)
1363 ;; Do these even if the parts turn out not to be integrable.
1364 ((eq (car-safe expr) '+)
1365 (math-add (calcFunc-integ (nth 1 expr) var low high)
1366 (calcFunc-integ (nth 2 expr) var low high)))
1367 ((eq (car-safe expr) '-)
1368 (math-sub (calcFunc-integ (nth 1 expr) var low high)
1369 (calcFunc-integ (nth 2 expr) var low high)))
1370 ((eq (car-safe expr) 'neg)
1371 (math-neg (calcFunc-integ (nth 1 expr) var low high)))
1372 ((and (eq (car-safe expr) '*)
1373 (not (math-expr-contains (nth 1 expr) var)))
1374 (math-mul (nth 1 expr) (calcFunc-integ (nth 2 expr) var low high)))
1375 ((and (eq (car-safe expr) '*)
1376 (not (math-expr-contains (nth 2 expr) var)))
1377 (math-mul (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr)))
1378 ((and (eq (car-safe expr) '/)
1379 (not (math-expr-contains (nth 1 expr) var))
1380 (not (math-equal-int (nth 1 expr) 1)))
1381 (math-mul (nth 1 expr)
1382 (calcFunc-integ (math-div 1 (nth 2 expr)) var low high)))
1383 ((and (eq (car-safe expr) '/)
1384 (not (math-expr-contains (nth 2 expr) var)))
1385 (math-div (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr)))
1386 ((and (eq (car-safe expr) '/)
1387 (eq (car-safe (nth 1 expr)) '*)
1388 (not (math-expr-contains (nth 1 (nth 1 expr)) var)))
1389 (math-mul (nth 1 (nth 1 expr))
1390 (calcFunc-integ (math-div (nth 2 (nth 1 expr)) (nth 2 expr))
1392 ((and (eq (car-safe expr) '/)
1393 (eq (car-safe (nth 1 expr)) '*)
1394 (not (math-expr-contains (nth 2 (nth 1 expr)) var)))
1395 (math-mul (nth 2 (nth 1 expr))
1396 (calcFunc-integ (math-div (nth 1 (nth 1 expr)) (nth 2 expr))
1398 ((and (eq (car-safe expr) '/)
1399 (eq (car-safe (nth 2 expr)) '*)
1400 (not (math-expr-contains (nth 1 (nth 2 expr)) var)))
1401 (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 2 (nth 2 expr)))
1403 (nth 1 (nth 2 expr))))
1404 ((and (eq (car-safe expr) '/)
1405 (eq (car-safe (nth 2 expr)) '*)
1406 (not (math-expr-contains (nth 2 (nth 2 expr)) var)))
1407 (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 1 (nth 2 expr)))
1409 (nth 2 (nth 2 expr))))
1410 ((eq (car-safe expr) 'vec)
1411 (cons 'vec (mapcar (function (lambda (x) (calcFunc-integ x var low high)))
1414 (let ((state (list calc-angle-mode
1415 ;;calc-symbolic-mode
1418 (calc-var-value 'var-IntegRules)
1419 (calc-var-value 'var-IntegSimpRules))))
1420 (or (equal state math-integral-cache-state)
1421 (setq math-integral-cache-state state
1422 math-integral-cache nil)))
1423 (let* ((math-max-integral-limit (or (and (boundp 'var-IntegLimit)
1424 (natnump var-IntegLimit)
1427 (math-integral-limit 1)
1428 (sexpr (math-expr-subst expr var math-integ-var))
1429 (trace-buffer (get-buffer "*Trace*"))
1430 (calc-language (if (eq calc-language 'big) nil calc-language))
1432 (math-enable-subst nil)
1433 (math-prev-parts-v nil)
1434 (math-doing-parts nil)
1435 (math-good-parts nil)
1438 (let ((calcbuf (current-buffer))
1439 (calcwin (selected-window)))
1442 (if (get-buffer-window trace-buffer)
1443 (select-window (get-buffer-window trace-buffer)))
1444 (set-buffer trace-buffer)
1445 (goto-char (point-max))
1446 (or (assq 'scroll-stop (buffer-local-variables))
1448 (make-local-variable 'scroll-step)
1449 (setq scroll-step 3)))
1451 (set-buffer calcbuf)
1452 (math-try-integral sexpr))
1453 (select-window calcwin)
1454 (set-buffer calcbuf)))
1455 (math-try-integral sexpr))))
1458 (if (calc-has-rules 'var-IntegAfterRules)
1459 (setq res (math-rewrite res '(var IntegAfterRules
1460 var-IntegAfterRules))))
1463 (math-sub (math-expr-subst res math-integ-var high)
1464 (math-expr-subst res math-integ-var low))
1465 (setq res (math-fix-const-terms res math-integ-vars))
1467 (math-expr-subst res math-integ-var low)
1468 (math-expr-subst res math-integ-var var)))))
1469 (append (list 'calcFunc-integ expr var)
1470 (and low (list low))
1471 (and high (list high)))))))
1475 (math-defintegral calcFunc-inv
1476 (math-integral (math-div 1 u)))
1478 (math-defintegral calcFunc-conj
1479 (let ((int (math-integral u)))
1481 (list 'calcFunc-conj int))))
1483 (math-defintegral calcFunc-deg
1484 (let ((int (math-integral u)))
1486 (list 'calcFunc-deg int))))
1488 (math-defintegral calcFunc-rad
1489 (let ((int (math-integral u)))
1491 (list 'calcFunc-rad int))))
1493 (math-defintegral calcFunc-re
1494 (let ((int (math-integral u)))
1496 (list 'calcFunc-re int))))
1498 (math-defintegral calcFunc-im
1499 (let ((int (math-integral u)))
1501 (list 'calcFunc-im int))))
1503 (math-defintegral calcFunc-sqrt
1504 (and (equal u math-integ-var)
1505 (math-mul '(frac 2 3)
1506 (list 'calcFunc-sqrt (math-pow u 3)))))
1508 (math-defintegral calcFunc-exp
1509 (or (and (equal u math-integ-var)
1510 (list 'calcFunc-exp u))
1511 (let ((p (math-is-polynomial u math-integ-var 2)))
1513 (let ((sqa (math-sqrt (math-neg (nth 2 p)))))
1516 (math-mul (math-div (list 'calcFunc-sqrt '(var pi var-pi))
1520 (math-div (math-sub (math-mul (car p)
1523 (math-sqr (nth 1 p))
1527 (math-sub (math-mul sqa math-integ-var)
1528 (math-div (nth 1 p) (math-mul 2 sqa)))))
1531 (math-defintegral calcFunc-ln
1532 (or (and (equal u math-integ-var)
1533 (math-sub (math-mul u (list 'calcFunc-ln u)) u))
1534 (and (eq (car u) '*)
1535 (math-integral (math-add (list 'calcFunc-ln (nth 1 u))
1536 (list 'calcFunc-ln (nth 2 u)))))
1537 (and (eq (car u) '/)
1538 (math-integral (math-sub (list 'calcFunc-ln (nth 1 u))
1539 (list 'calcFunc-ln (nth 2 u)))))
1540 (and (eq (car u) '^)
1541 (math-integral (math-mul (nth 2 u)
1542 (list 'calcFunc-ln (nth 1 u)))))))
1544 (math-defintegral calcFunc-log10
1545 (and (equal u math-integ-var)
1546 (math-sub (math-mul u (list 'calcFunc-ln u))
1547 (math-div u (list 'calcFunc-ln 10)))))
1549 (math-defintegral-2 calcFunc-log
1550 (math-integral (math-div (list 'calcFunc-ln u)
1551 (list 'calcFunc-ln v))))
1553 (math-defintegral calcFunc-sin
1554 (or (and (equal u math-integ-var)
1555 (math-neg (math-from-radians-2 (list 'calcFunc-cos u))))
1556 (and (nth 2 (math-is-polynomial u math-integ-var 2))
1557 (math-integral (math-to-exponentials (list 'calcFunc-sin u))))))
1559 (math-defintegral calcFunc-cos
1560 (or (and (equal u math-integ-var)
1561 (math-from-radians-2 (list 'calcFunc-sin u)))
1562 (and (nth 2 (math-is-polynomial u math-integ-var 2))
1563 (math-integral (math-to-exponentials (list 'calcFunc-cos u))))))
1565 (math-defintegral calcFunc-tan
1566 (and (equal u math-integ-var)
1567 (math-neg (math-from-radians-2
1568 (list 'calcFunc-ln (list 'calcFunc-cos u))))))
1570 (math-defintegral calcFunc-arcsin
1571 (and (equal u math-integ-var)
1572 (math-add (math-mul u (list 'calcFunc-arcsin u))
1573 (math-from-radians-2
1574 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))))
1576 (math-defintegral calcFunc-arccos
1577 (and (equal u math-integ-var)
1578 (math-sub (math-mul u (list 'calcFunc-arccos u))
1579 (math-from-radians-2
1580 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))))
1582 (math-defintegral calcFunc-arctan
1583 (and (equal u math-integ-var)
1584 (math-sub (math-mul u (list 'calcFunc-arctan u))
1585 (math-from-radians-2
1586 (math-div (list 'calcFunc-ln (math-add 1 (math-sqr u)))
1589 (math-defintegral calcFunc-sinh
1590 (and (equal u math-integ-var)
1591 (list 'calcFunc-cosh u)))
1593 (math-defintegral calcFunc-cosh
1594 (and (equal u math-integ-var)
1595 (list 'calcFunc-sinh u)))
1597 (math-defintegral calcFunc-tanh
1598 (and (equal u math-integ-var)
1599 (list 'calcFunc-ln (list 'calcFunc-cosh u))))
1601 (math-defintegral calcFunc-arcsinh
1602 (and (equal u math-integ-var)
1603 (math-sub (math-mul u (list 'calcFunc-arcsinh u))
1604 (list 'calcFunc-sqrt (math-add (math-sqr u) 1)))))
1606 (math-defintegral calcFunc-arccosh
1607 (and (equal u math-integ-var)
1608 (math-sub (math-mul u (list 'calcFunc-arccosh u))
1609 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))
1611 (math-defintegral calcFunc-arctanh
1612 (and (equal u math-integ-var)
1613 (math-sub (math-mul u (list 'calcFunc-arctan u))
1614 (math-div (list 'calcFunc-ln
1615 (math-add 1 (math-sqr u)))
1618 ;;; (Ax + B) / (ax^2 + bx + c)^n forms.
1619 (math-defintegral-2 /
1620 (math-integral-rational-funcs u v))
1622 (defun math-integral-rational-funcs (u v)
1623 (let ((pu (math-is-polynomial u math-integ-var 1))
1627 (if (and (eq (car-safe v) '^) (natnump (nth 2 v)))
1628 (setq vpow (nth 2 v)
1630 (and (setq pv (math-is-polynomial v math-integ-var 2))
1631 (let ((int (math-mul-thru
1633 (math-integral-q02 (car pv) (nth 1 pv)
1634 (nth 2 pv) v vpow))))
1636 (setq int (math-add int
1641 (nth 2 pv) v vpow)))))
1644 (defun math-integral-q12 (a b c v vpow)
1648 (math-sub (math-div math-integ-var b)
1649 (math-mul (math-div a (math-sqr b))
1650 (list 'calcFunc-ln v))))
1652 (math-div (math-add (list 'calcFunc-ln v)
1656 (let ((nm1 (math-sub vpow 1))
1657 (nm2 (math-sub vpow 2)))
1659 (math-div a (math-mul nm1 (math-pow v nm1)))
1660 (math-div 1 (math-mul nm2 (math-pow v nm2))))
1663 (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b))))
1664 (let ((part (math-div b (math-mul 2 c))))
1665 (math-mul-thru (math-pow c vpow)
1666 (math-integral-q12 part 1 nil
1667 (math-add math-integ-var part)
1670 (and (math-ratp q) (math-negp q)
1671 (let ((calc-symbolic-mode t))
1672 (math-ratp (math-sqrt (math-neg q))))
1673 (throw 'int-rat nil)) ; should have used calcFunc-apart first
1674 (math-sub (math-div (list 'calcFunc-ln v) (math-mul 2 c))
1675 (math-mul-thru (math-div b (math-mul 2 c))
1676 (math-integral-q02 a b c v 1))))
1678 (let ((n (1- vpow)))
1679 (math-sub (math-neg (math-div
1680 (math-add (math-mul b math-integ-var)
1682 (math-mul n (math-mul q (math-pow v n)))))
1683 (math-mul-thru (math-div (math-mul b (1- (* 2 n)))
1685 (math-integral-q02 a b c v n)))))))
1688 (defun math-integral-q02 (a b c v vpow)
1692 (math-div (list 'calcFunc-ln v) b))
1694 (math-div (math-pow v (- 1 vpow))
1695 (math-mul (- 1 vpow) b)))))
1697 (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b))))
1698 (let ((part (math-div b (math-mul 2 c))))
1699 (math-mul-thru (math-pow c vpow)
1700 (math-integral-q02 part 1 nil
1701 (math-add math-integ-var part)
1704 (setq part (math-add (math-mul 2 (math-mul c math-integ-var)) b))
1706 (let ((n (1- vpow)))
1707 (math-add (math-div part (math-mul n (math-mul q (math-pow v n))))
1708 (math-mul-thru (math-div (math-mul (- (* 4 n) 2) c)
1710 (math-integral-q02 a b c v n)))))
1711 ((math-guess-if-neg q)
1712 (setq rq (list 'calcFunc-sqrt (math-neg q)))
1713 ;;(math-div-thru (list 'calcFunc-ln
1714 ;; (math-div (math-sub part rq)
1715 ;; (math-add part rq)))
1717 (math-div (math-mul -2 (list 'calcFunc-arctanh
1718 (math-div part rq)))
1721 (setq rq (list 'calcFunc-sqrt q))
1722 (math-div (math-mul 2 (math-to-radians-2
1723 (list 'calcFunc-arctan
1724 (math-div part rq))))
1729 (math-defintegral calcFunc-erf
1730 (and (equal u math-integ-var)
1731 (math-add (math-mul u (list 'calcFunc-erf u))
1732 (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u))
1733 (list 'calcFunc-sqrt
1734 '(var pi var-pi)))))))
1736 (math-defintegral calcFunc-erfc
1737 (and (equal u math-integ-var)
1738 (math-sub (math-mul u (list 'calcFunc-erfc u))
1739 (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u))
1740 (list 'calcFunc-sqrt
1741 '(var pi var-pi)))))))
1746 (defun calcFunc-table (expr var &optional low high step)
1747 (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
1748 (or high (setq high low low 1))
1749 (and (or (math-infinitep low) (math-infinitep high))
1751 (math-scan-for-limits expr))
1752 (and step (math-zerop step) (math-reject-arg step 'nonzerop))
1753 (let ((known (+ (if (Math-objectp low) 1 0)
1754 (if (Math-objectp high) 1 0)
1755 (if (or (null step) (Math-objectp step)) 1 0)))
1756 (count '(var inf var-inf))
1758 (or (= known 2) ; handy optimization
1759 (equal high '(var inf var-inf))
1761 (setq count (math-div (math-sub high low) (or step 1)))
1762 (or (Math-objectp count)
1763 (setq count (math-simplify count)))
1764 (if (Math-messy-integerp count)
1765 (setq count (math-trunc count)))))
1766 (if (Math-negp count)
1768 (if (integerp count)
1769 (let ((var-DUMMY nil)
1770 (vec math-tabulate-initial)
1771 (math-working-step-2 (1+ count))
1772 (math-working-step 0))
1773 (setq expr (math-evaluate-expr
1774 (math-expr-subst expr var '(var DUMMY var-DUMMY))))
1776 (setq math-working-step (1+ math-working-step)
1778 vec (cond ((eq math-tabulate-function 'calcFunc-sum)
1779 (math-add vec (math-evaluate-expr expr)))
1780 ((eq math-tabulate-function 'calcFunc-prod)
1781 (math-mul vec (math-evaluate-expr expr)))
1783 (cons (math-evaluate-expr expr) vec)))
1784 low (math-add low (or step 1))
1786 (if math-tabulate-function
1788 (cons 'vec (nreverse vec))))
1789 (if (Math-integerp count)
1790 (calc-record-why 'fixnump high)
1791 (if (Math-num-integerp low)
1792 (if (Math-num-integerp high)
1793 (calc-record-why 'integerp step)
1794 (calc-record-why 'integerp high))
1795 (calc-record-why 'integerp low)))
1796 (append (list (or math-tabulate-function 'calcFunc-table)
1798 (and (not (and (equal low '(neg (var inf var-inf)))
1799 (equal high '(var inf var-inf))))
1801 (and step (list step)))))
1804 (setq math-tabulate-initial nil)
1805 (setq math-tabulate-function nil)
1807 (defun math-scan-for-limits (x)
1808 (cond ((Math-primp x))
1809 ((and (eq (car x) 'calcFunc-subscr)
1810 (Math-vectorp (nth 1 x))
1811 (math-expr-contains (nth 2 x) var))
1812 (let* ((calc-next-why nil)
1813 (low-val (math-solve-for (nth 2 x) 1 var nil))
1814 (high-val (math-solve-for (nth 2 x) (1- (length (nth 1 x)))
1817 (and low-val (math-realp low-val)
1818 high-val (math-realp high-val))
1819 (and (Math-lessp high-val low-val)
1820 (setq temp low-val low-val high-val high-val temp))
1821 (setq low (math-max low (math-ceiling low-val))
1822 high (math-min high (math-floor high-val)))))
1824 (while (setq x (cdr x))
1825 (math-scan-for-limits (car x)))))
1829 (defun calcFunc-sum (expr var &optional low high step)
1830 (if math-disable-sums (math-reject-arg))
1831 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2)))
1832 (math-sum-rec expr var low high step)))
1833 (math-disable-sums t))
1834 (math-normalize res))
1836 (setq math-disable-sums nil)
1838 (defun math-sum-rec (expr var &optional low high step)
1839 (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
1840 (and low (not high) (setq high low low 1))
1844 ((not (math-expr-contains expr var))
1845 (math-mul expr (math-add (math-div (math-sub high low) (or step 1))
1847 ((and step (not (math-equal-int step 1)))
1848 (if (math-negp step)
1849 (math-sum-rec expr var high low (math-neg step))
1850 (let ((lo (math-simplify (math-div low step))))
1851 (if (math-known-num-integerp lo)
1852 (math-sum-rec (math-normalize
1853 (math-expr-subst expr var
1854 (math-mul step var)))
1855 var lo (math-simplify (math-div high step)))
1856 (math-sum-rec (math-normalize
1857 (math-expr-subst expr var
1858 (math-add (math-mul step var)
1861 (math-simplify (math-div (math-sub high low)
1863 ((memq (setq t1 (math-compare low high)) '(0 1))
1865 (math-expr-subst expr var low)
1867 ((setq t1 (math-is-polynomial expr var 20))
1871 (setq poly (math-poly-mix poly 1
1872 (math-sum-integer-power n) (car t1))
1875 (setq n (math-build-polynomial-expr poly high))
1876 (if (memq low '(0 1))
1878 (math-sub n (math-build-polynomial-expr poly
1879 (math-sub low 1))))))
1880 ((and (memq (car expr) '(+ -))
1881 (setq t1 (math-sum-rec (nth 1 expr) var low high)
1882 t2 (math-sum-rec (nth 2 expr) var low high))
1883 (not (and (math-expr-calls t1 '(calcFunc-sum))
1884 (math-expr-calls t2 '(calcFunc-sum)))))
1885 (list (car expr) t1 t2))
1886 ((and (eq (car expr) '*)
1887 (setq t1 (math-sum-const-factors expr var)))
1888 (math-mul (car t1) (math-sum-rec (cdr t1) var low high)))
1889 ((and (eq (car expr) '*) (memq (car-safe (nth 1 expr)) '(+ -)))
1890 (math-sum-rec (math-add-or-sub (math-mul (nth 1 (nth 1 expr))
1892 (math-mul (nth 2 (nth 1 expr))
1894 nil (eq (car (nth 1 expr)) '-))
1896 ((and (eq (car expr) '*) (memq (car-safe (nth 2 expr)) '(+ -)))
1897 (math-sum-rec (math-add-or-sub (math-mul (nth 1 expr)
1898 (nth 1 (nth 2 expr)))
1899 (math-mul (nth 1 expr)
1900 (nth 2 (nth 2 expr)))
1901 nil (eq (car (nth 2 expr)) '-))
1903 ((and (eq (car expr) '/)
1904 (not (math-primp (nth 1 expr)))
1905 (setq t1 (math-sum-const-factors (nth 1 expr) var)))
1907 (math-sum-rec (math-div (cdr t1) (nth 2 expr))
1909 ((and (eq (car expr) '/)
1910 (setq t1 (math-sum-const-factors (nth 2 expr) var)))
1911 (math-div (math-sum-rec (math-div (nth 1 expr) (cdr t1))
1914 ((eq (car expr) 'neg)
1915 (math-neg (math-sum-rec (nth 1 expr) var low high)))
1916 ((and (eq (car expr) '^)
1917 (not (math-expr-contains (nth 1 expr) var))
1918 (setq t1 (math-is-polynomial (nth 2 expr) var 1)))
1919 (let ((x (math-pow (nth 1 expr) (nth 1 t1))))
1920 (math-div (math-mul (math-sub (math-pow x (math-add 1 high))
1922 (math-pow (nth 1 expr) (car t1)))
1924 ((and (setq t1 (math-to-exponentials expr))
1925 (setq t1 (math-sum-rec t1 var low high))
1926 (not (math-expr-calls t1 '(calcFunc-sum))))
1928 ((memq (car expr) '(calcFunc-ln calcFunc-log10))
1929 (list (car expr) (calcFunc-prod (nth 1 expr) var low high)))
1930 ((and (eq (car expr) 'calcFunc-log)
1932 (not (math-expr-contains (nth 2 expr) var)))
1934 (calcFunc-prod (nth 1 expr) var low high)
1936 (if (equal val '(var nan var-nan)) (setq val nil))
1938 (let* ((math-tabulate-initial 0)
1939 (math-tabulate-function 'calcFunc-sum))
1940 (calcFunc-table expr var low high))))
1943 (defun calcFunc-asum (expr var low &optional high step no-mul-flag)
1944 (or high (setq high low low 1))
1945 (if (and step (not (math-equal-int step 1)))
1946 (if (math-negp step)
1947 (math-mul (math-pow -1 low)
1948 (calcFunc-asum expr var high low (math-neg step) t))
1949 (let ((lo (math-simplify (math-div low step))))
1950 (if (math-num-integerp lo)
1951 (calcFunc-asum (math-normalize
1952 (math-expr-subst expr var
1953 (math-mul step var)))
1954 var lo (math-simplify (math-div high step)))
1955 (calcFunc-asum (math-normalize
1956 (math-expr-subst expr var
1957 (math-add (math-mul step var)
1960 (math-simplify (math-div (math-sub high low)
1962 (math-mul (if no-mul-flag 1 (math-pow -1 low))
1963 (calcFunc-sum (math-mul (math-pow -1 var) expr) var low high)))
1966 (defun math-sum-const-factors (expr var)
1970 (while (eq (car-safe p) '*)
1971 (if (math-expr-contains (nth 1 p) var)
1972 (setq not-const (cons (nth 1 p) not-const))
1973 (setq const (cons (nth 1 p) const)))
1975 (if (math-expr-contains p var)
1976 (setq not-const (cons p not-const))
1977 (setq const (cons p const)))
1979 (cons (let ((temp (car const)))
1980 (while (setq const (cdr const))
1981 (setq temp (list '* (car const) temp)))
1983 (let ((temp (or (car not-const) 1)))
1984 (while (setq not-const (cdr not-const))
1985 (setq temp (list '* (car not-const) temp)))
1989 ;; Following is from CRC Math Tables, 27th ed, pp. 52-53.
1990 (defun math-sum-integer-power (pow)
1991 (let ((calc-prefer-frac t)
1992 (n (length math-sum-int-pow-cache)))
1994 (let* ((new (list 0 0))
1996 (pp (cdr (nth (1- n) math-sum-int-pow-cache)))
2001 (setq q (math-div (car pp) p)
2002 new (cons (math-mul q n) new)
2003 sum (math-add sum q)
2006 (setcar lin (math-sub 1 (math-mul n sum)))
2007 (setq math-sum-int-pow-cache
2008 (nconc math-sum-int-pow-cache (list (nreverse new)))
2010 (nth pow math-sum-int-pow-cache))
2012 (setq math-sum-int-pow-cache (list '(0 1)))
2014 (defun math-to-exponentials (expr)
2017 (let ((x (nth 1 expr))
2018 (pi (if calc-symbolic-mode '(var pi var-pi) (math-pi)))
2019 (i (if calc-symbolic-mode '(var i var-i) '(cplx 0 1))))
2020 (cond ((eq (car expr) 'calcFunc-exp)
2021 (list '^ '(var e var-e) x))
2022 ((eq (car expr) 'calcFunc-sin)
2023 (or (eq calc-angle-mode 'rad)
2024 (setq x (list '/ (list '* x pi) 180)))
2026 (list '^ '(var e var-e) (list '* x i))
2027 (list '^ '(var e var-e)
2028 (list 'neg (list '* x i))))
2030 ((eq (car expr) 'calcFunc-cos)
2031 (or (eq calc-angle-mode 'rad)
2032 (setq x (list '/ (list '* x pi) 180)))
2034 (list '^ '(var e var-e)
2036 (list '^ '(var e var-e)
2037 (list 'neg (list '* x i))))
2039 ((eq (car expr) 'calcFunc-sinh)
2041 (list '^ '(var e var-e) x)
2042 (list '^ '(var e var-e) (list 'neg x)))
2044 ((eq (car expr) 'calcFunc-cosh)
2046 (list '^ '(var e var-e) x)
2047 (list '^ '(var e var-e) (list 'neg x)))
2052 (defun math-to-exps (expr)
2053 (cond (calc-symbolic-mode expr)
2055 (if (equal expr '(var e var-e)) (math-e) expr))
2056 ((and (eq (car expr) '^)
2057 (equal (nth 1 expr) '(var e var-e)))
2058 (list 'calcFunc-exp (nth 2 expr)))
2060 (cons (car expr) (mapcar 'math-to-exps (cdr expr)))))
2064 (defun calcFunc-prod (expr var &optional low high step)
2065 (if math-disable-prods (math-reject-arg))
2066 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2)))
2067 (math-prod-rec expr var low high step)))
2068 (math-disable-prods t))
2069 (math-normalize res))
2071 (setq math-disable-prods nil)
2073 (defun math-prod-rec (expr var &optional low high step)
2074 (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
2075 (and low (not high) (setq high '(var inf var-inf)))
2079 ((not (math-expr-contains expr var))
2080 (math-pow expr (math-add (math-div (math-sub high low) (or step 1))
2082 ((and step (not (math-equal-int step 1)))
2083 (if (math-negp step)
2084 (math-prod-rec expr var high low (math-neg step))
2085 (let ((lo (math-simplify (math-div low step))))
2086 (if (math-known-num-integerp lo)
2087 (math-prod-rec (math-normalize
2088 (math-expr-subst expr var
2089 (math-mul step var)))
2090 var lo (math-simplify (math-div high step)))
2091 (math-prod-rec (math-normalize
2092 (math-expr-subst expr var
2093 (math-add (math-mul step
2097 (math-simplify (math-div (math-sub high low)
2099 ((and (memq (car expr) '(* /))
2100 (setq t1 (math-prod-rec (nth 1 expr) var low high)
2101 t2 (math-prod-rec (nth 2 expr) var low high))
2102 (not (and (math-expr-calls t1 '(calcFunc-prod))
2103 (math-expr-calls t2 '(calcFunc-prod)))))
2104 (list (car expr) t1 t2))
2105 ((and (eq (car expr) '^)
2106 (not (math-expr-contains (nth 2 expr) var)))
2107 (math-pow (math-prod-rec (nth 1 expr) var low high)
2109 ((and (eq (car expr) '^)
2110 (not (math-expr-contains (nth 1 expr) var)))
2111 (math-pow (nth 1 expr)
2112 (calcFunc-sum (nth 2 expr) var low high)))
2113 ((eq (car expr) 'sqrt)
2114 (math-normalize (list 'calcFunc-sqrt
2115 (list 'calcFunc-prod (nth 1 expr)
2117 ((eq (car expr) 'neg)
2118 (math-mul (math-pow -1 (math-add (math-sub high low) 1))
2119 (math-prod-rec (nth 1 expr) var low high)))
2120 ((eq (car expr) 'calcFunc-exp)
2121 (list 'calcFunc-exp (calcFunc-sum (nth 1 expr) var low high)))
2122 ((and (setq t1 (math-is-polynomial expr var 1))
2125 ((or (and (math-equal-int (nth 1 t1) 1)
2126 (setq low (math-simplify
2127 (math-add low (car t1)))
2129 (math-add high (car t1)))))
2130 (and (math-equal-int (nth 1 t1) -1)
2133 (math-sub (car t1) high))
2135 (math-sub (car t1) t2)))))
2136 (if (or (math-zerop low) (math-zerop high))
2138 (if (and (or (math-negp low) (math-negp high))
2139 (or (math-num-integerp low)
2140 (math-num-integerp high)))
2141 (if (math-posp high)
2143 (math-mul (math-pow -1
2145 (math-add low high) 1))
2147 (list 'calcFunc-fact
2149 (list 'calcFunc-fact
2150 (math-sub -1 high)))))
2152 (list 'calcFunc-fact high)
2153 (list 'calcFunc-fact (math-sub low 1))))))
2154 ((and (or (and (math-equal-int (nth 1 t1) 2)
2155 (setq t2 (math-simplify
2156 (math-add (math-mul low 2)
2159 (math-add (math-mul high 2)
2161 (and (math-equal-int (nth 1 t1) -2)
2162 (setq t2 (math-simplify
2169 (or (math-integerp t2)
2170 (and (math-messy-integerp t2)
2171 (setq t2 (math-trunc t2)))
2173 (and (math-messy-integerp t3)
2174 (setq t3 (math-trunc t3)))))
2175 (if (or (math-zerop t2) (math-zerop t3))
2177 (if (or (math-evenp t2) (math-evenp t3))
2178 (if (or (math-negp t2) (math-negp t3))
2179 (if (math-posp high)
2182 (list 'calcFunc-dfact
2184 (list 'calcFunc-dfact
2187 (list 'calcFunc-dfact t3)
2188 (list 'calcFunc-dfact
2193 (list '/ (list '- (list '- t2 t3)
2197 (list 'calcFunc-dfact
2199 (list 'calcFunc-dfact
2203 (list 'calcFunc-dfact t3)
2204 (list 'calcFunc-dfact
2208 (if (equal val '(var nan var-nan)) (setq val nil))
2210 (let* ((math-tabulate-initial 1)
2211 (math-tabulate-function 'calcFunc-prod))
2212 (calcFunc-table expr var low high))))
2218 ;;; Attempt to reduce lhs = rhs to solve-var = rhs', where solve-var appears
2219 ;;; in lhs but not in rhs or rhs'; return rhs'.
2220 ;;; Uses global values: solve-*.
2221 (defun math-try-solve-for (lhs rhs &optional sign no-poly)
2223 (cond ((equal lhs solve-var)
2224 (setq math-solve-sign sign)
2225 (if (eq solve-full 'all)
2226 (let ((vec (list 'vec (math-evaluate-expr rhs)))
2228 (while math-solve-ranges
2229 (setq p (car math-solve-ranges)
2232 (while (setq p (cdr p))
2233 (setq newvec (nconc newvec
2234 (cdr (math-expr-subst
2235 vec var (car p))))))
2237 math-solve-ranges (cdr math-solve-ranges)))
2238 (math-normalize vec))
2242 ((and (eq (car lhs) '-)
2243 (eq (car-safe (nth 1 lhs)) (car-safe (nth 2 lhs)))
2245 (= (length (nth 1 lhs)) 2)
2246 (= (length (nth 2 lhs)) 2)
2247 (setq t1 (get (car (nth 1 lhs)) 'math-inverse))
2248 (setq t2 (funcall t1 '(var SOLVEDUM SOLVEDUM)))
2249 (eq (math-expr-contains-count t2 '(var SOLVEDUM SOLVEDUM)) 1)
2250 (setq t3 (math-solve-above-dummy t2))
2251 (setq t1 (math-try-solve-for (math-sub (nth 1 (nth 1 lhs))
2254 (nth 1 (nth 2 lhs))))
2257 ((eq (car lhs) 'neg)
2258 (math-try-solve-for (nth 1 lhs) (math-neg rhs)
2259 (and sign (- sign))))
2260 ((and (not (eq solve-full 't)) (math-try-solve-prod)))
2262 (setq t2 (math-decompose-poly lhs solve-var 15 rhs)))
2263 (setq t1 (cdr (nth 1 t2))
2264 t1 (let ((math-solve-ranges math-solve-ranges))
2265 (cond ((= (length t1) 5)
2266 (apply 'math-solve-quartic (car t2) t1))
2268 (apply 'math-solve-cubic (car t2) t1))
2270 (apply 'math-solve-quadratic (car t2) t1))
2272 (apply 'math-solve-linear (car t2) sign t1))
2274 (math-poly-all-roots (car t2) t1))
2275 (calc-symbolic-mode nil)
2279 (math-poly-any-root (reverse t1) 0 t)
2282 (if (eq (nth 2 t2) 1)
2284 (math-solve-prod t1 (math-try-solve-for (nth 2 t2) 0 nil t)))
2285 (calc-record-why "*Unable to find a symbolic solution")
2287 ((and (math-solve-find-root-term lhs nil)
2288 (eq (math-expr-contains-count lhs t1) 1)) ; just in case
2289 (math-try-solve-for (math-simplify
2290 (math-sub (if (or t3 (math-evenp t2))
2292 (math-neg (math-pow t1 t2)))
2294 (math-sub (math-normalize
2301 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
2302 (math-try-solve-for (nth 2 lhs)
2303 (math-sub rhs (nth 1 lhs))
2305 ((not (math-expr-contains (nth 2 lhs) solve-var))
2306 (math-try-solve-for (nth 1 lhs)
2307 (math-sub rhs (nth 2 lhs))
2309 ((eq (car lhs) 'calcFunc-eq)
2310 (math-try-solve-for (math-sub (nth 1 lhs) (nth 2 lhs))
2313 (cond ((or (and (eq (car-safe (nth 1 lhs)) 'calcFunc-sin)
2314 (eq (car-safe (nth 2 lhs)) 'calcFunc-cos))
2315 (and (eq (car-safe (nth 1 lhs)) 'calcFunc-cos)
2316 (eq (car-safe (nth 2 lhs)) 'calcFunc-sin)))
2317 (math-try-solve-for (math-sub (nth 1 lhs)
2318 (list (car (nth 1 lhs))
2320 (math-quarter-circle t)
2321 (nth 1 (nth 2 lhs)))))
2323 ((not (math-expr-contains (nth 1 lhs) solve-var))
2324 (math-try-solve-for (nth 2 lhs)
2325 (math-sub (nth 1 lhs) rhs)
2326 (and sign (- sign))))
2327 ((not (math-expr-contains (nth 2 lhs) solve-var))
2328 (math-try-solve-for (nth 1 lhs)
2329 (math-add rhs (nth 2 lhs))
2331 ((and (eq solve-full 't) (math-try-solve-prod)))
2332 ((and (eq (car lhs) '%)
2333 (not (math-expr-contains (nth 2 lhs) solve-var)))
2334 (math-try-solve-for (nth 1 lhs) (math-add rhs
2337 ((eq (car lhs) 'calcFunc-log)
2338 (cond ((not (math-expr-contains (nth 2 lhs) solve-var))
2339 (math-try-solve-for (nth 1 lhs) (math-pow (nth 2 lhs) rhs)))
2340 ((not (math-expr-contains (nth 1 lhs) solve-var))
2341 (math-try-solve-for (nth 2 lhs) (math-pow
2343 (math-div 1 rhs))))))
2344 ((and (= (length lhs) 2)
2346 (setq t1 (get (car lhs) 'math-inverse))
2347 (setq t2 (funcall t1 rhs)))
2348 (setq t1 (get (car lhs) 'math-inverse-sign))
2349 (math-try-solve-for (nth 1 lhs) (math-normalize t2)
2353 (funcall t1 lhs sign)))))
2354 ((and (symbolp (car lhs))
2355 (setq t1 (get (car lhs) 'math-inverse-n))
2356 (setq t2 (funcall t1 lhs rhs)))
2358 ((setq t1 (math-expand-formula lhs))
2359 (math-try-solve-for t1 rhs sign))
2361 (calc-record-why "*No inverse known" lhs)
2365 (setq math-solve-ranges nil)
2367 (defun math-try-solve-prod ()
2368 (cond ((eq (car lhs) '*)
2369 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
2370 (math-try-solve-for (nth 2 lhs)
2371 (math-div rhs (nth 1 lhs))
2372 (math-solve-sign sign (nth 1 lhs))))
2373 ((not (math-expr-contains (nth 2 lhs) solve-var))
2374 (math-try-solve-for (nth 1 lhs)
2375 (math-div rhs (nth 2 lhs))
2376 (math-solve-sign sign (nth 2 lhs))))
2378 (math-solve-prod (let ((math-solve-ranges math-solve-ranges))
2379 (math-try-solve-for (nth 2 lhs) 0))
2380 (math-try-solve-for (nth 1 lhs) 0)))))
2382 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
2383 (math-try-solve-for (nth 2 lhs)
2384 (math-div (nth 1 lhs) rhs)
2385 (math-solve-sign sign (nth 1 lhs))))
2386 ((not (math-expr-contains (nth 2 lhs) solve-var))
2387 (math-try-solve-for (nth 1 lhs)
2388 (math-mul rhs (nth 2 lhs))
2389 (math-solve-sign sign (nth 2 lhs))))
2390 ((setq t1 (math-try-solve-for (math-sub (nth 1 lhs)
2391 (math-mul (nth 2 lhs)
2396 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
2399 (math-add (math-normalize
2400 (list 'calcFunc-log rhs (nth 1 lhs)))
2403 (math-mul '(var pi var-pi)
2407 (list 'calcFunc-ln (nth 1 lhs)))))))
2408 ((not (math-expr-contains (nth 2 lhs) solve-var))
2409 (cond ((and (integerp (nth 2 lhs))
2411 (setq t1 (math-integer-log2 (nth 2 lhs))))
2413 (if (and (eq solve-full t)
2414 (math-known-realp (nth 1 lhs)))
2416 (while (>= (setq t1 (1- t1)) 0)
2417 (setq t2 (list 'calcFunc-sqrt t2)))
2418 (setq t2 (math-solve-get-sign t2)))
2419 (while (>= (setq t1 (1- t1)) 0)
2420 (setq t2 (math-solve-get-sign
2422 (list 'calcFunc-sqrt t2))))))
2425 (math-normalize t2)))
2426 ((math-looks-negp (nth 2 lhs))
2428 (list '^ (nth 1 lhs) (math-neg (nth 2 lhs)))
2430 ((and (eq solve-full t)
2431 (Math-integerp (nth 2 lhs))
2432 (math-known-realp (nth 1 lhs)))
2433 (setq t1 (math-normalize
2434 (list 'calcFunc-nroot rhs (nth 2 lhs))))
2435 (if (math-evenp (nth 2 lhs))
2436 (setq t1 (math-solve-get-sign t1)))
2440 (math-oddp (nth 2 lhs))
2441 (math-solve-sign sign (nth 2 lhs)))))
2442 (t (math-try-solve-for
2447 (if (Math-realp (nth 2 lhs))
2452 (and (integerp (nth 2 lhs))
2455 (math-div (nth 2 lhs) 2))
2462 (and (integerp (nth 2 lhs))
2467 (list 'calcFunc-nroot
2471 (math-oddp (nth 2 lhs))
2472 (math-solve-sign sign (nth 2 lhs)))))))))
2476 (defun math-solve-prod (lsoln rsoln)
2481 ((eq solve-full 'all)
2482 (cons 'vec (append (cdr lsoln) (cdr rsoln))))
2485 (list 'calcFunc-gt (math-solve-get-sign 1) 0)
2491 ;;; This deals with negative, fractional, and symbolic powers of "x".
2492 (defun math-solve-poly-funny-powers (sub-rhs) ; uses "t1", "t2"
2494 (let ((pp math-poly-neg-powers)
2497 (setq fac (math-pow (car pp) (or math-poly-mult-powers 1))
2498 t1 (math-mul t1 fac)
2499 rhs (math-mul rhs fac)
2501 (if sub-rhs (setq t1 (math-sub t1 rhs)))
2502 (let ((math-poly-neg-powers nil))
2503 (setq t2 (math-mul (or math-poly-mult-powers 1)
2504 (let ((calc-prefer-frac t))
2505 (math-div 1 math-poly-frac-powers)))
2506 t1 (math-is-polynomial (math-simplify (calcFunc-expand t1)) b 50)))
2509 ;;; This converts "a x^8 + b x^5 + c x^2" to "(a (x^3)^2 + b (x^3) + c) * x^2".
2510 (defun math-solve-crunch-poly (max-degree) ; uses "t1", "t3"
2512 (while (and t1 (Math-zerop (car t1)))
2516 (let* ((degree (1- (length t1)))
2518 (while (and (> scale 1) (= (car t3) 1))
2519 (and (= (% degree scale) 0)
2525 (if (= (% n scale) 0)
2526 (setq new-t1 (nconc new-t1 (list (car p))))
2527 (or (Math-zerop (car p))
2532 (setq t3 (cons scale (cdr t3))
2534 (setq scale (1- scale)))
2535 (setq t3 (list (math-mul (car t3) t2) (math-mul count t2)))
2536 (<= (1- (length t1)) max-degree))))
2539 (defun calcFunc-poly (expr var &optional degree)
2541 (or (natnump degree) (math-reject-arg degree 'fixnatnump))
2543 (let ((p (math-is-polynomial expr var degree 'gen)))
2548 (math-reject-arg expr "Expected a polynomial")))
2551 (defun calcFunc-gpoly (expr var &optional degree)
2553 (or (natnump degree) (math-reject-arg degree 'fixnatnump))
2555 (let* ((math-poly-base-variable var)
2556 (d (math-decompose-poly expr var degree nil)))
2559 (math-reject-arg expr "Expected a polynomial")))
2562 (defun math-decompose-poly (lhs solve-var degree sub-rhs)
2563 (let ((rhs (or sub-rhs 1))
2565 (setq t2 (math-polynomial-base
2569 (let ((math-poly-neg-powers '(1))
2570 (math-poly-mult-powers nil)
2571 (math-poly-frac-powers 1)
2572 (math-poly-exp-base t))
2573 (and (not (equal b lhs))
2574 (or (not (memq (car-safe b) '(+ -))) sub-rhs)
2575 (setq t3 '(1 0) t2 1
2576 t1 (math-is-polynomial lhs b 50))
2577 (if (and (equal math-poly-neg-powers '(1))
2578 (memq math-poly-mult-powers '(nil 1))
2579 (eq math-poly-frac-powers 1)
2581 (setq t1 (cons (math-sub (car t1) rhs)
2583 (math-solve-poly-funny-powers sub-rhs))
2584 (math-solve-crunch-poly degree)
2585 (or (math-expr-contains b solve-var)
2586 (math-expr-contains (car t3) solve-var))))))))
2588 (list (math-pow t2 (car t3))
2591 (math-pow t2 (nth 1 t3))
2592 (math-div (math-pow t2 (nth 1 t3)) rhs)))))
2595 (defun math-solve-linear (var sign b a)
2596 (math-try-solve-for var
2597 (math-div (math-neg b) a)
2598 (math-solve-sign sign a)
2602 (defun math-solve-quadratic (var c b a)
2605 (if (math-looks-evenp b)
2606 (let ((halfb (math-div b 2)))
2610 (math-solve-get-sign
2612 (list 'calcFunc-sqrt
2613 (math-add (math-sqr halfb)
2614 (math-mul (math-neg c) a))))))
2619 (math-solve-get-sign
2621 (list 'calcFunc-sqrt
2622 (math-add (math-sqr b)
2623 (math-mul 4 (math-mul (math-neg c) a)))))))
2628 (defun math-solve-cubic (var d c b a)
2629 (let* ((p (math-div b a))
2633 (aa (math-sub q (math-div psqr 3)))
2635 (math-div (math-sub (math-mul 2 (math-mul psqr p))
2636 (math-mul 9 (math-mul p q)))
2640 (math-try-solve-for (math-pow (math-add var (math-div p 3)) 3)
2641 (math-neg bb) nil t)
2644 (math-mul (math-add var (math-div p 3))
2645 (math-add (math-sqr (math-add var (math-div p 3)))
2648 (setq m (math-mul 2 (list 'calcFunc-sqrt (math-div aa -3))))
2657 (math-sub (list 'calcFunc-arccos
2658 (math-div (math-mul 3 bb)
2662 (math-add 1 (math-solve-get-int
2665 calc-symbolic-mode))))
2671 (defun math-solve-quartic (var d c b a aa)
2672 (setq a (math-div a aa))
2673 (setq b (math-div b aa))
2674 (setq c (math-div c aa))
2675 (setq d (math-div d aa))
2678 (let* ((asqr (math-sqr a))
2679 (asqr4 (math-div asqr 4))
2680 (y (let ((solve-full nil)
2682 (math-solve-cubic solve-var
2684 (math-mul 4 (math-mul b d))
2687 (math-sub (math-mul a c)
2691 (rsqr (math-add (math-sub asqr4 b) y))
2692 (r (list 'calcFunc-sqrt rsqr))
2693 (sign1 (math-solve-get-sign 1))
2694 (de (list 'calcFunc-sqrt
2696 (math-sub (math-mul 3 asqr4)
2698 (if (Math-zerop rsqr)
2702 (list 'calcFunc-sqrt
2703 (math-sub (math-sqr y)
2709 (math-mul 4 (math-mul a b))
2715 (math-sub (math-add (math-mul sign1 (math-div r 2))
2716 (math-solve-get-sign (math-div de 2)))
2721 (defun math-poly-all-roots (var p &optional math-factoring)
2723 (let* ((math-symbolic-solve calc-symbolic-mode)
2725 (deg (1- (length p)))
2726 (orig-p (reverse p))
2727 (math-int-coefs nil)
2728 (math-int-scale nil)
2729 (math-double-roots nil)
2730 (math-int-factors nil)
2731 (math-int-threshold nil)
2733 ;; If rational coefficients, look for exact rational factors.
2734 (while (and pp (Math-ratp (car pp)))
2737 (if (or math-factoring math-symbolic-solve)
2739 (let ((lead (car orig-p))
2740 (calc-prefer-frac t)
2741 (scale (apply 'math-lcm-denoms p)))
2742 (setq math-int-scale (math-abs (math-mul scale lead))
2743 math-int-threshold (math-div '(float 5 -2) math-int-scale)
2744 math-int-coefs (cdr (math-div (cons 'vec orig-p) lead)))))
2746 (let ((calc-prefer-frac nil)
2747 (calc-symbolic-mode nil)
2749 (def-p (copy-sequence orig-p)))
2751 (if (Math-numberp (car pp))
2754 (while (> deg (if math-symbolic-solve 2 4))
2755 (let* ((x (math-poly-any-root def-p '(float 0 0) nil))
2757 (if (and (eq (car-safe x) 'cplx)
2758 (math-nearly-zerop (nth 2 x) (nth 1 x)))
2759 (setq x (calcFunc-re x)))
2761 (setq roots (cons x roots)))
2762 (or (math-numberp x)
2763 (setq x (math-evaluate-expr x)))
2766 (while (setq pp (cdr pp))
2769 (setq b (math-add (math-mul x b) c)))
2770 (setq def-p (cdr def-p)
2772 (setq p (reverse def-p))))
2774 (let ((solve-var '(var DUMMY var-DUMMY))
2775 (math-solve-sign nil)
2776 (math-solve-ranges nil)
2778 (if (= (length p) (length math-int-coefs))
2779 (setq p (reverse math-int-coefs)))
2780 (setq roots (append (cdr (apply (cond ((= deg 2)
2781 'math-solve-quadratic)
2785 'math-solve-quartic))
2789 (setq roots (cons (math-div (math-neg (car p)) (nth 1 p))
2794 (math-poly-integer-root (car roots))
2795 (setq roots (cdr roots)))
2796 (list math-int-factors (nreverse math-int-coefs) math-int-scale))
2797 (let ((vec nil) res)
2799 (let ((root (car roots))
2800 (solve-full (and solve-full 'all)))
2801 (if (math-floatp root)
2802 (setq root (math-poly-any-root orig-p root t)))
2803 (setq vec (append vec
2804 (cdr (or (math-try-solve-for var root nil t)
2805 (throw 'ouch nil))))))
2806 (setq roots (cdr roots)))
2807 (setq vec (cons 'vec (nreverse vec)))
2808 (if math-symbolic-solve
2809 (setq vec (math-normalize vec)))
2810 (if (eq solve-full t)
2811 (list 'calcFunc-subscr
2813 (math-solve-get-int 1 (1- (length orig-p)) 1))
2816 (setq math-symbolic-solve nil)
2818 (defun math-lcm-denoms (&rest fracs)
2821 (if (eq (car-safe (car fracs)) 'frac)
2822 (setq den (calcFunc-lcm den (nth 2 (car fracs)))))
2823 (setq fracs (cdr fracs)))
2827 (defun math-poly-any-root (p x polish) ; p is a reverse poly coeff list
2828 (let* ((newt (if (math-zerop x)
2829 (math-poly-newton-root
2830 p '(cplx (float 123 -6) (float 1 -4)) 4)
2831 (math-poly-newton-root p x 4)))
2832 (res (if (math-zerop (cdr newt))
2834 (if (and (math-lessp (cdr newt) '(float 1 -3)) (not polish))
2835 (setq newt (math-poly-newton-root p (car newt) 30)))
2836 (if (math-zerop (cdr newt))
2838 (math-poly-laguerre-root p x polish)))))
2839 (and math-symbolic-solve (math-floatp res)
2844 (defun math-poly-newton-root (p x iters)
2845 (let* ((calc-prefer-frac nil)
2846 (calc-symbolic-mode nil)
2847 (try-integer math-int-coefs)
2849 (while (and (> (setq iters (1- iters)) 0)
2851 (math-working "newton" x)
2854 (while (setq pp (cdr pp))
2855 (setq d (math-add (math-mul x d) b)
2856 b (math-add (math-mul x b) (car pp))))
2857 (not (math-zerop d)))
2859 (setq dx (math-div b d)
2862 (let ((adx (math-abs-approx dx)))
2863 (and (math-lessp adx math-int-threshold)
2864 (let ((iroot (math-poly-integer-root x)))
2867 (setq try-integer nil))))))
2868 (or (not (or (eq dx 0)
2869 (math-nearly-zerop dx (math-abs-approx x))))
2870 (progn (setq dx 0) nil)))))
2871 (cons x (if (math-zerop x)
2872 1 (math-div (math-abs-approx dx) (math-abs-approx x)))))
2875 (defun math-poly-integer-root (x)
2876 (and (math-lessp (calcFunc-xpon (math-abs-approx x)) calc-internal-prec)
2878 (let* ((calc-prefer-frac t)
2879 (xre (calcFunc-re x))
2880 (xim (calcFunc-im x))
2881 (xresq (math-sqr xre))
2882 (ximsq (math-sqr xim)))
2883 (if (math-lessp ximsq (calcFunc-scf xresq -1))
2884 ;; Look for linear factor
2885 (let* ((rnd (math-div (math-round (math-mul xre math-int-scale))
2887 (icp math-int-coefs)
2890 (while (setq icp (cdr icp))
2891 (setq newcoef (cons rem newcoef)
2892 rem (math-add (car icp)
2893 (math-mul rem rnd))))
2894 (and (math-zerop rem)
2896 (setq math-int-coefs (nreverse newcoef)
2897 math-int-factors (cons (list (math-neg rnd))
2900 ;; Look for irreducible quadratic factor
2901 (let* ((rnd1 (math-div (math-round
2902 (math-mul xre (math-mul -2 math-int-scale)))
2904 (sqscale (math-sqr math-int-scale))
2905 (rnd0 (math-div (math-round (math-mul (math-add xresq ximsq)
2908 (rem1 (car math-int-coefs))
2909 (icp (cdr math-int-coefs))
2912 (found (assoc (list rnd0 rnd1 (math-posp xim))
2916 (setq math-double-roots (delq found math-double-roots)
2918 (while (setq icp (cdr icp))
2920 newcoef (cons rem1 newcoef)
2921 rem1 (math-sub rem0 (math-mul this rnd1))
2922 rem0 (math-sub (car icp) (math-mul this rnd0)))))
2923 (and (math-zerop rem0)
2925 (let ((aa (math-div rnd1 -2)))
2926 (or found (setq math-int-coefs (reverse newcoef)
2927 math-double-roots (cons (list
2932 math-int-factors (cons (cons rnd0 rnd1)
2935 (let ((calc-symbolic-mode math-symbolic-solve))
2936 (math-mul (math-sqrt (math-sub (math-sqr aa)
2938 (if (math-negp xim) -1 1))))))))))
2940 (setq math-int-coefs nil)
2942 ;;; The following routine is from Numerical Recipes, section 9.5.
2943 (defun math-poly-laguerre-root (p x polish)
2944 (let* ((calc-prefer-frac nil)
2945 (calc-symbolic-mode nil)
2948 (try-newt (not polish))
2952 (and (or (< (setq iters (1+ iters)) 50)
2953 (math-reject-arg x "*Laguerre's method failed to converge"))
2954 (let ((err (math-abs-approx (car p)))
2955 (abx (math-abs-approx x))
2959 (while (setq pp (cdr pp))
2960 (setq f (math-add (math-mul x f) d)
2961 d (math-add (math-mul x d) b)
2962 b (math-add (math-mul x b) (car pp))
2963 err (math-add (math-abs-approx b) (math-mul abx err))))
2964 (math-lessp (calcFunc-scf err (- -2 calc-internal-prec))
2965 (math-abs-approx b)))
2966 (or (not (math-zerop d))
2967 (not (math-zerop f))
2969 (setq x (math-pow (math-neg b) (list 'frac 1 m)))
2971 (let* ((g (math-div d b))
2973 (h (math-sub g2 (math-mul 2 (math-div f b))))
2975 (math-mul (1- m) (math-sub (math-mul m h) g2))))
2976 (gp (math-add g sq))
2977 (gm (math-sub g sq)))
2978 (if (math-lessp (calcFunc-abssqr gp) (calcFunc-abssqr gm))
2980 (setq dx (math-div m gp)
2983 (math-lessp (math-abs-approx dx)
2984 (calcFunc-scf (math-abs-approx x) -3)))
2985 (let ((newt (math-poly-newton-root p x1 7)))
2988 (if (math-zerop (cdr newt))
2989 (setq x (car newt) x1 x)
2990 (if (math-lessp (cdr newt) '(float 1 -6))
2991 (let ((newt2 (math-poly-newton-root
2993 (if (math-zerop (cdr newt2))
2994 (setq x (car newt2) x1 x)
2995 (setq x (car newt))))))))
2997 (math-nearly-equal x x1))))
2998 (let ((cdx (math-abs-approx dx)))
3003 (math-lessp cdx dxold)
3006 (let ((digs (calcFunc-xpon
3007 (math-div (math-abs-approx x) cdx))))
3009 "*Could not attain full precision")
3011 (let ((calc-internal-prec (max 3 digs)))
3012 (setq x (math-normalize x))))))
3016 (math-lessp (calcFunc-scf (math-abs-approx x)
3017 (- calc-internal-prec))
3019 (or (and (math-floatp x)
3020 (math-poly-integer-root x))
3024 (defun math-solve-above-dummy (x)
3025 (and (not (Math-primp x))
3026 (if (and (equal (nth 1 x) '(var SOLVEDUM SOLVEDUM))
3030 (while (and (setq x (cdr x))
3031 (not (setq res (math-solve-above-dummy (car x))))))
3035 (defun math-solve-find-root-term (x neg) ; sets "t2", "t3"
3036 (if (math-solve-find-root-in-prod x)
3039 (and (memq (car-safe x) '(+ -))
3040 (or (math-solve-find-root-term (nth 1 x) neg)
3041 (math-solve-find-root-term (nth 2 x)
3042 (if (eq (car x) '-) (not neg) neg)))))
3045 (defun math-solve-find-root-in-prod (x)
3047 (math-expr-contains x solve-var)
3048 (or (and (eq (car x) 'calcFunc-sqrt)
3050 (and (eq (car x) '^)
3051 (or (and (memq (math-quarter-integer (nth 2 x)) '(1 2 3))
3053 (and (eq (car-safe (nth 2 x)) 'frac)
3054 (eq (nth 2 (nth 2 x)) 3)
3056 (and (memq (car x) '(* /))
3057 (or (and (not (math-expr-contains (nth 1 x) solve-var))
3058 (math-solve-find-root-in-prod (nth 2 x)))
3059 (and (not (math-expr-contains (nth 2 x) solve-var))
3060 (math-solve-find-root-in-prod (nth 1 x)))))))
3064 (defun math-solve-system (exprs solve-vars solve-full)
3065 (setq exprs (mapcar 'list (if (Math-vectorp exprs)
3068 solve-vars (if (Math-vectorp solve-vars)
3071 (or (let ((math-solve-simplifying nil))
3072 (math-solve-system-rec exprs solve-vars nil))
3073 (let ((math-solve-simplifying t))
3074 (math-solve-system-rec exprs solve-vars nil)))
3077 ;;; The following backtracking solver works by choosing a variable
3078 ;;; and equation, and trying to solve the equation for the variable.
3079 ;;; If it succeeds it calls itself recursively with that variable and
3080 ;;; equation removed from their respective lists, and with the solution
3081 ;;; added to solns as well as being substituted into all existing
3082 ;;; equations. The algorithm terminates when any solution path
3083 ;;; manages to remove all the variables from var-list.
3085 ;;; To support calcFunc-roots, entries in eqn-list and solns are
3086 ;;; actually lists of equations.
3088 (defun math-solve-system-rec (eqn-list var-list solns)
3093 ;; Try each variable in turn.
3099 (elim (eq (car-safe vv) 'calcFunc-elim)))
3101 (setq vv (nth 1 vv)))
3103 ;; Try each equation in turn.
3112 ;; Try to solve for vv the list of equations e2.
3114 (setq res2 (or (and (eq (car e2) eprev)
3116 (math-solve-for (car e2) 0 vv
3118 (setq eprev (car e2)
3119 res (cons (if (eq solve-full 'all)
3127 ;; Found a solution. Now try other variables.
3128 (setq res (nreverse res)
3129 res (math-solve-system-rec
3131 'math-solve-system-subst
3133 (copy-sequence eqn-list)))
3134 (delq (car v) (copy-sequence var-list))
3135 (let ((math-solve-simplifying nil)
3141 (math-solve-system-subst
3146 (cons (cons vv (apply 'append res))
3154 ;; Eliminated all variables, so now put solution into the proper format.
3155 (setq solns (sort solns
3158 (not (memq (car x) (memq (car y) solve-vars)))))))
3159 (if (eq solve-full 'all)
3164 (mapcar (function (lambda (x) (cons 'vec (cdr x)))) solns)
3165 (mapcar (function (lambda (x) (cons 'vec x))) eqn-list)))))
3169 (mapcar (function (lambda (x) (cons 'calcFunc-eq x))) solns)
3170 (mapcar 'car eqn-list))))))
3173 (defun math-solve-system-subst (x) ; uses "res" and "v"
3177 (setq accum (nconc accum
3180 (if math-solve-simplifying
3182 (math-expr-subst (car x) vv r))
3183 (math-expr-subst (car x) vv r))))
3191 (defun math-get-from-counter (name)
3192 (let ((ctr (assq name calc-command-flags)))
3194 (setcdr ctr (1+ (cdr ctr)))
3195 (setq ctr (cons name 1)
3196 calc-command-flags (cons ctr calc-command-flags)))
3200 (defun math-solve-get-sign (val)
3201 (setq val (math-simplify val))
3202 (if (and (eq (car-safe val) '*)
3203 (Math-numberp (nth 1 val)))
3204 (list '* (nth 1 val) (math-solve-get-sign (nth 2 val)))
3205 (and (eq (car-safe val) 'calcFunc-sqrt)
3206 (eq (car-safe (nth 1 val)) '^)
3207 (setq val (math-normalize (list '^
3209 (math-div (nth 2 (nth 1 val)) 2)))))
3211 (if (and (calc-var-value 'var-GenCount)
3212 (Math-natnump var-GenCount)
3213 (not (eq solve-full 'all)))
3215 (math-mul (list 'calcFunc-as var-GenCount) val)
3216 (setq var-GenCount (math-add var-GenCount 1))
3217 (calc-refresh-evaltos 'var-GenCount))
3218 (let* ((var (concat "s" (int-to-string (math-get-from-counter 'solve-sign))))
3219 (var2 (list 'var (intern var) (intern (concat "var-" var)))))
3220 (if (eq solve-full 'all)
3221 (setq math-solve-ranges (cons (list var2 1 -1)
3222 math-solve-ranges)))
3223 (math-mul var2 val)))
3224 (calc-record-why "*Choosing positive solution")
3228 (defun math-solve-get-int (val &optional range first)
3230 (if (and (calc-var-value 'var-GenCount)
3231 (Math-natnump var-GenCount)
3232 (not (eq solve-full 'all)))
3234 (math-mul val (list 'calcFunc-an var-GenCount))
3235 (setq var-GenCount (math-add var-GenCount 1))
3236 (calc-refresh-evaltos 'var-GenCount))
3237 (let* ((var (concat "n" (int-to-string (math-get-from-counter 'solve-int))))
3238 (var2 (list 'var (intern var) (intern (concat "var-" var)))))
3239 (if (and range (eq solve-full 'all))
3240 (setq math-solve-ranges (cons (cons var2
3241 (cdr (calcFunc-index
3242 range (or first 0))))
3243 math-solve-ranges)))
3244 (math-mul val var2)))
3245 (calc-record-why "*Choosing 0 for arbitrary integer in solution")
3249 (defun math-solve-sign (sign expr)
3251 (let ((s1 (math-possible-signs expr)))
3252 (cond ((memq s1 '(4 6))
3258 (defun math-looks-evenp (expr)
3259 (if (Math-integerp expr)
3261 (if (memq (car expr) '(* /))
3262 (math-looks-evenp (nth 1 expr))))
3265 (defun math-solve-for (lhs rhs solve-var solve-full &optional sign)
3266 (if (math-expr-contains rhs solve-var)
3267 (math-solve-for (math-sub lhs rhs) 0 solve-var solve-full)
3268 (and (math-expr-contains lhs solve-var)
3269 (math-with-extra-prec 1
3270 (let* ((math-poly-base-variable solve-var)
3271 (res (math-try-solve-for lhs rhs sign)))
3272 (if (and (eq solve-full 'all)
3273 (math-known-realp solve-var))
3274 (let ((old-len (length res))
3279 (and (not (memq (car-safe x)
3283 new-len (length res))
3284 (if (< new-len old-len)
3285 (calc-record-why (if (= new-len 1)
3286 "*All solutions were complex"
3288 "*Omitted %d complex solutions"
3289 (- old-len new-len)))))))
3293 (defun math-solve-eqn (expr var full)
3294 (if (memq (car-safe expr) '(calcFunc-neq calcFunc-lt calcFunc-gt
3295 calcFunc-leq calcFunc-geq))
3296 (let ((res (math-solve-for (cons '- (cdr expr))
3298 (if (eq (car expr) 'calcFunc-neq) nil 1))))
3300 (if (eq math-solve-sign 1)
3301 (list (car expr) var res)
3302 (if (eq math-solve-sign -1)
3303 (list (car expr) res var)
3304 (or (eq (car expr) 'calcFunc-neq)
3306 "*Can't determine direction of inequality"))
3307 (and (memq (car expr) '(calcFunc-neq calcFunc-lt calcFunc-gt))
3308 (list 'calcFunc-neq var res))))))
3309 (let ((res (math-solve-for expr 0 var full)))
3311 (list 'calcFunc-eq var res))))
3314 (defun math-reject-solution (expr var func)
3315 (if (math-expr-contains expr var)
3316 (or (equal (car calc-next-why) '(* "Unable to find a symbolic solution"))
3317 (calc-record-why "*Unable to find a solution")))
3318 (list func expr var)
3321 (defun calcFunc-solve (expr var)
3322 (or (if (or (Math-vectorp expr) (Math-vectorp var))
3323 (math-solve-system expr var nil)
3324 (math-solve-eqn expr var nil))
3325 (math-reject-solution expr var 'calcFunc-solve))
3328 (defun calcFunc-fsolve (expr var)
3329 (or (if (or (Math-vectorp expr) (Math-vectorp var))
3330 (math-solve-system expr var t)
3331 (math-solve-eqn expr var t))
3332 (math-reject-solution expr var 'calcFunc-fsolve))
3335 (defun calcFunc-roots (expr var)
3336 (let ((math-solve-ranges nil))
3337 (or (if (or (Math-vectorp expr) (Math-vectorp var))
3338 (math-solve-system expr var 'all)
3339 (math-solve-for expr 0 var 'all))
3340 (math-reject-solution expr var 'calcFunc-roots)))
3343 (defun calcFunc-finv (expr var)
3344 (let ((res (math-solve-for expr math-integ-var var nil)))
3346 (math-normalize (math-expr-subst res math-integ-var var))
3347 (math-reject-solution expr var 'calcFunc-finv)))
3350 (defun calcFunc-ffinv (expr var)
3351 (let ((res (math-solve-for expr math-integ-var var t)))
3353 (math-normalize (math-expr-subst res math-integ-var var))
3354 (math-reject-solution expr var 'calcFunc-finv)))
3358 (put 'calcFunc-inv 'math-inverse
3359 (function (lambda (x) (math-div 1 x))))
3360 (put 'calcFunc-inv 'math-inverse-sign -1)
3362 (put 'calcFunc-sqrt 'math-inverse
3363 (function (lambda (x) (math-sqr x))))
3365 (put 'calcFunc-conj 'math-inverse
3366 (function (lambda (x) (list 'calcFunc-conj x))))
3368 (put 'calcFunc-abs 'math-inverse
3369 (function (lambda (x) (math-solve-get-sign x))))
3371 (put 'calcFunc-deg 'math-inverse
3372 (function (lambda (x) (list 'calcFunc-rad x))))
3373 (put 'calcFunc-deg 'math-inverse-sign 1)
3375 (put 'calcFunc-rad 'math-inverse
3376 (function (lambda (x) (list 'calcFunc-deg x))))
3377 (put 'calcFunc-rad 'math-inverse-sign 1)
3379 (put 'calcFunc-ln 'math-inverse
3380 (function (lambda (x) (list 'calcFunc-exp x))))
3381 (put 'calcFunc-ln 'math-inverse-sign 1)
3383 (put 'calcFunc-log10 'math-inverse
3384 (function (lambda (x) (list 'calcFunc-exp10 x))))
3385 (put 'calcFunc-log10 'math-inverse-sign 1)
3387 (put 'calcFunc-lnp1 'math-inverse
3388 (function (lambda (x) (list 'calcFunc-expm1 x))))
3389 (put 'calcFunc-lnp1 'math-inverse-sign 1)
3391 (put 'calcFunc-exp 'math-inverse
3392 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-ln x))
3394 (math-mul '(var pi var-pi)
3396 '(var i var-i))))))))
3397 (put 'calcFunc-exp 'math-inverse-sign 1)
3399 (put 'calcFunc-expm1 'math-inverse
3400 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-lnp1 x))
3402 (math-mul '(var pi var-pi)
3404 '(var i var-i))))))))
3405 (put 'calcFunc-expm1 'math-inverse-sign 1)
3407 (put 'calcFunc-sin 'math-inverse
3408 (function (lambda (x) (let ((n (math-solve-get-int 1)))
3409 (math-add (math-mul (math-normalize
3410 (list 'calcFunc-arcsin x))
3412 (math-mul (math-half-circle t)
3415 (put 'calcFunc-cos 'math-inverse
3416 (function (lambda (x) (math-add (math-solve-get-sign
3418 (list 'calcFunc-arccos x)))
3420 (math-full-circle t))))))
3422 (put 'calcFunc-tan 'math-inverse
3423 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-arctan x))
3425 (math-half-circle t))))))
3427 (put 'calcFunc-arcsin 'math-inverse
3428 (function (lambda (x) (math-normalize (list 'calcFunc-sin x)))))
3430 (put 'calcFunc-arccos 'math-inverse
3431 (function (lambda (x) (math-normalize (list 'calcFunc-cos x)))))
3433 (put 'calcFunc-arctan 'math-inverse
3434 (function (lambda (x) (math-normalize (list 'calcFunc-tan x)))))
3436 (put 'calcFunc-sinh 'math-inverse
3437 (function (lambda (x) (let ((n (math-solve-get-int 1)))
3438 (math-add (math-mul (math-normalize
3439 (list 'calcFunc-arcsinh x))
3441 (math-mul (math-half-circle t)
3445 (put 'calcFunc-sinh 'math-inverse-sign 1)
3447 (put 'calcFunc-cosh 'math-inverse
3448 (function (lambda (x) (math-add (math-solve-get-sign
3450 (list 'calcFunc-arccosh x)))
3451 (math-mul (math-full-circle t)
3453 '(var i var-i)))))))
3455 (put 'calcFunc-tanh 'math-inverse
3456 (function (lambda (x) (math-add (math-normalize
3457 (list 'calcFunc-arctanh x))
3458 (math-mul (math-half-circle t)
3460 '(var i var-i)))))))
3461 (put 'calcFunc-tanh 'math-inverse-sign 1)
3463 (put 'calcFunc-arcsinh 'math-inverse
3464 (function (lambda (x) (math-normalize (list 'calcFunc-sinh x)))))
3465 (put 'calcFunc-arcsinh 'math-inverse-sign 1)
3467 (put 'calcFunc-arccosh 'math-inverse
3468 (function (lambda (x) (math-normalize (list 'calcFunc-cosh x)))))
3470 (put 'calcFunc-arctanh 'math-inverse
3471 (function (lambda (x) (math-normalize (list 'calcFunc-tanh x)))))
3472 (put 'calcFunc-arctanh 'math-inverse-sign 1)
3476 (defun calcFunc-taylor (expr var num)
3477 (let ((x0 0) (v var))
3478 (if (memq (car-safe var) '(+ - calcFunc-eq))
3479 (setq x0 (if (eq (car var) '+) (math-neg (nth 2 var)) (nth 2 var))
3481 (or (and (eq (car-safe v) 'var)
3482 (math-expr-contains expr v)
3484 (let ((accum (math-expr-subst expr v x0))
3485 (var2 (if (eq (car var) 'calcFunc-eq)
3491 (while (and (<= (setq n (1+ n)) num)
3492 (setq fprime (calcFunc-deriv fprime v nil t)))
3493 (setq fprime (math-simplify fprime)
3494 nfac (math-mul nfac n)
3495 accum (math-add accum
3496 (math-div (math-mul (math-pow var2 n)
3501 (math-normalize accum))))
3502 (list 'calcFunc-taylor expr var num)))