1 ;;; calcalg2.el --- more algebraic functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
6 ;; Maintainers: D. Goel <deego@gnufans.org>
7 ;; Colin Walters <walters@debian.org>
9 ;; This file is part of GNU Emacs.
11 ;; GNU Emacs is distributed in the hope that it will be useful,
12 ;; but WITHOUT ANY WARRANTY. No author or distributor
13 ;; accepts responsibility to anyone for the consequences of using it
14 ;; or for whether it serves any particular purpose or works at all,
15 ;; unless he says so in writing. Refer to the GNU Emacs General Public
16 ;; License for full details.
18 ;; Everyone is granted permission to copy, modify and redistribute
19 ;; GNU Emacs, but only under the conditions described in the
20 ;; GNU Emacs General Public License. A copy of this license is
21 ;; supposed to have been given to you along with GNU Emacs so you
22 ;; can know your rights and responsibilities. It should be in a
23 ;; file named COPYING. Among other things, the copyright notice
24 ;; and this notice must be preserved on all copies.
30 ;; This file is autoloaded from calc-ext.el.
35 (defun calc-Need-calc-alg-2 () nil)
38 (defun calc-derivative (var num)
39 (interactive "sDifferentiate with respect to: \np")
42 (error "Order of derivative must be positive"))
43 (let ((func (if (calc-is-hyperbolic) 'calcFunc-tderiv 'calcFunc-deriv))
45 (if (or (equal var "") (equal var "$"))
49 (setq var (math-read-expr var))
50 (when (eq (car-safe var) 'error)
51 (error "Bad format in expression: %s" (nth 1 var)))
54 (while (>= (setq num (1- num)) 0)
55 (setq expr (list func expr var)))
56 (calc-enter-result n "derv" expr))))
58 (defun calc-integral (var)
59 (interactive "sIntegration variable: ")
61 (if (or (equal var "") (equal var "$"))
62 (calc-enter-result 2 "intg" (list 'calcFunc-integ
65 (let ((var (math-read-expr var)))
66 (if (eq (car-safe var) 'error)
67 (error "Bad format in expression: %s" (nth 1 var)))
68 (calc-enter-result 1 "intg" (list 'calcFunc-integ
72 (defun calc-num-integral (&optional varname lowname highname)
73 (interactive "sIntegration variable: ")
74 (calc-tabular-command 'calcFunc-ninteg "Integration" "nint"
75 nil varname lowname highname))
77 (defun calc-summation (arg &optional varname lowname highname)
78 (interactive "P\nsSummation variable: ")
79 (calc-tabular-command 'calcFunc-sum "Summation" "sum"
80 arg varname lowname highname))
82 (defun calc-alt-summation (arg &optional varname lowname highname)
83 (interactive "P\nsSummation variable: ")
84 (calc-tabular-command 'calcFunc-asum "Summation" "asum"
85 arg varname lowname highname))
87 (defun calc-product (arg &optional varname lowname highname)
88 (interactive "P\nsIndex variable: ")
89 (calc-tabular-command 'calcFunc-prod "Index" "prod"
90 arg varname lowname highname))
92 (defun calc-tabulate (arg &optional varname lowname highname)
93 (interactive "P\nsIndex variable: ")
94 (calc-tabular-command 'calcFunc-table "Index" "tabl"
95 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)))))))
155 (defun calc-solve-for (var)
156 (interactive "sVariable to solve for: ")
158 (let ((func (if (calc-is-inverse)
159 (if (calc-is-hyperbolic) 'calcFunc-ffinv 'calcFunc-finv)
160 (if (calc-is-hyperbolic) 'calcFunc-fsolve 'calcFunc-solve))))
161 (if (or (equal var "") (equal var "$"))
162 (calc-enter-result 2 "solv" (list func
165 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
166 (not (string-match "\\[" var)))
167 (math-read-expr (concat "[" var "]"))
168 (math-read-expr var))))
169 (if (eq (car-safe var) 'error)
170 (error "Bad format in expression: %s" (nth 1 var)))
171 (calc-enter-result 1 "solv" (list func
175 (defun calc-poly-roots (var)
176 (interactive "sVariable to solve for: ")
178 (if (or (equal var "") (equal var "$"))
179 (calc-enter-result 2 "prts" (list 'calcFunc-roots
182 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
183 (not (string-match "\\[" var)))
184 (math-read-expr (concat "[" var "]"))
185 (math-read-expr var))))
186 (if (eq (car-safe var) 'error)
187 (error "Bad format in expression: %s" (nth 1 var)))
188 (calc-enter-result 1 "prts" (list 'calcFunc-roots
192 (defun calc-taylor (var nterms)
193 (interactive "sTaylor expansion variable: \nNNumber of terms: ")
195 (let ((var (math-read-expr var)))
196 (if (eq (car-safe var) 'error)
197 (error "Bad format in expression: %s" (nth 1 var)))
198 (calc-enter-result 1 "tylr" (list 'calcFunc-taylor
201 (prefix-numeric-value nterms))))))
204 (defun math-derivative (expr) ; uses global values: deriv-var, deriv-total.
205 (cond ((equal expr deriv-var)
207 ((or (Math-scalarp expr)
208 (eq (car expr) 'sdev)
209 (and (eq (car expr) 'var)
210 (or (not deriv-total)
211 (math-const-var expr)
213 (math-setup-declarations)
214 (memq 'const (nth 1 (or (assq (nth 2 expr)
216 math-decls-all)))))))
219 (math-add (math-derivative (nth 1 expr))
220 (math-derivative (nth 2 expr))))
222 (math-sub (math-derivative (nth 1 expr))
223 (math-derivative (nth 2 expr))))
224 ((memq (car expr) '(calcFunc-eq calcFunc-neq calcFunc-lt
225 calcFunc-gt calcFunc-leq calcFunc-geq))
227 (math-derivative (nth 1 expr))
228 (math-derivative (nth 2 expr))))
229 ((eq (car expr) 'neg)
230 (math-neg (math-derivative (nth 1 expr))))
232 (math-add (math-mul (nth 2 expr)
233 (math-derivative (nth 1 expr)))
234 (math-mul (nth 1 expr)
235 (math-derivative (nth 2 expr)))))
237 (math-sub (math-div (math-derivative (nth 1 expr))
239 (math-div (math-mul (nth 1 expr)
240 (math-derivative (nth 2 expr)))
241 (math-sqr (nth 2 expr)))))
243 (let ((du (math-derivative (nth 1 expr)))
244 (dv (math-derivative (nth 2 expr))))
246 (setq du (math-mul (nth 2 expr)
247 (math-mul (math-normalize
250 (math-add (nth 2 expr) -1)))
253 (setq dv (math-mul (math-normalize
254 (list 'calcFunc-ln (nth 1 expr)))
255 (math-mul expr dv))))
258 (math-derivative (nth 1 expr))) ; a reasonable definition
259 ((eq (car expr) 'vec)
260 (math-map-vec 'math-derivative expr))
261 ((and (memq (car expr) '(calcFunc-conj calcFunc-re calcFunc-im))
263 (list (car expr) (math-derivative (nth 1 expr))))
264 ((and (memq (car expr) '(calcFunc-subscr calcFunc-mrow calcFunc-mcol))
266 (let ((d (math-derivative (nth 1 expr))))
268 0 ; assume x and x_1 are independent vars
269 (list (car expr) d (nth 2 expr)))))
270 (t (or (and (symbolp (car expr))
271 (if (= (length expr) 2)
272 (let ((handler (get (car expr) 'math-derivative)))
274 (let ((deriv (math-derivative (nth 1 expr))))
275 (if (Math-zerop deriv)
277 (math-mul (funcall handler (nth 1 expr))
279 (let ((handler (get (car expr) 'math-derivative-n)))
281 (funcall handler expr)))))
282 (and (not (eq deriv-symb 'pre-expand))
283 (let ((exp (math-expand-formula expr)))
285 (or (let ((deriv-symb 'pre-expand))
286 (catch 'math-deriv (math-derivative expr)))
287 (math-derivative exp)))))
288 (if (or (Math-objvecp expr)
290 (not (symbolp (car expr))))
292 (throw 'math-deriv nil)
293 (list (if deriv-total 'calcFunc-tderiv 'calcFunc-deriv)
300 (while (setq arg (cdr arg))
301 (or (Math-zerop (setq derv (math-derivative (car arg))))
302 (let ((func (intern (concat (symbol-name (car expr))
307 (prop (cond ((= (length expr) 2)
316 'math-derivative-5))))
322 (let ((handler (get func prop)))
323 (or (and prop handler
324 (apply handler (cdr expr)))
328 (throw 'math-deriv nil)
329 (cons func (cdr expr))))))))))
333 (defun calcFunc-deriv (expr deriv-var &optional deriv-value deriv-symb)
334 (let* ((deriv-total nil)
335 (res (catch 'math-deriv (math-derivative expr))))
336 (or (eq (car-safe res) 'calcFunc-deriv)
338 (setq res (math-normalize res)))
341 (math-expr-subst res deriv-var deriv-value)
344 (defun calcFunc-tderiv (expr deriv-var &optional deriv-value deriv-symb)
345 (math-setup-declarations)
346 (let* ((deriv-total t)
347 (res (catch 'math-deriv (math-derivative expr))))
348 (or (eq (car-safe res) 'calcFunc-tderiv)
350 (setq res (math-normalize res)))
353 (math-expr-subst res deriv-var deriv-value)
356 (put 'calcFunc-inv\' 'math-derivative-1
357 (function (lambda (u) (math-neg (math-div 1 (math-sqr u))))))
359 (put 'calcFunc-sqrt\' 'math-derivative-1
360 (function (lambda (u) (math-div 1 (math-mul 2 (list 'calcFunc-sqrt u))))))
362 (put 'calcFunc-deg\' 'math-derivative-1
363 (function (lambda (u) (math-div-float '(float 18 1) (math-pi)))))
365 (put 'calcFunc-rad\' 'math-derivative-1
366 (function (lambda (u) (math-pi-over-180))))
368 (put 'calcFunc-ln\' 'math-derivative-1
369 (function (lambda (u) (math-div 1 u))))
371 (put 'calcFunc-log10\' 'math-derivative-1
372 (function (lambda (u)
373 (math-div (math-div 1 (math-normalize '(calcFunc-ln 10)))
376 (put 'calcFunc-lnp1\' 'math-derivative-1
377 (function (lambda (u) (math-div 1 (math-add u 1)))))
379 (put 'calcFunc-log\' 'math-derivative-2
380 (function (lambda (x b)
381 (and (not (Math-zerop b))
382 (let ((lnv (math-normalize
383 (list 'calcFunc-ln b))))
384 (math-div 1 (math-mul lnv x)))))))
386 (put 'calcFunc-log\'2 'math-derivative-2
387 (function (lambda (x b)
388 (let ((lnv (list 'calcFunc-ln b)))
389 (math-neg (math-div (list 'calcFunc-log x b)
390 (math-mul lnv b)))))))
392 (put 'calcFunc-exp\' 'math-derivative-1
393 (function (lambda (u) (math-normalize (list 'calcFunc-exp u)))))
395 (put 'calcFunc-expm1\' 'math-derivative-1
396 (function (lambda (u) (math-normalize (list 'calcFunc-expm1 u)))))
398 (put 'calcFunc-sin\' 'math-derivative-1
399 (function (lambda (u) (math-to-radians-2 (math-normalize
400 (list 'calcFunc-cos u))))))
402 (put 'calcFunc-cos\' 'math-derivative-1
403 (function (lambda (u) (math-neg (math-to-radians-2
405 (list 'calcFunc-sin u)))))))
407 (put 'calcFunc-tan\' 'math-derivative-1
408 (function (lambda (u) (math-to-radians-2
409 (math-div 1 (math-sqr
411 (list 'calcFunc-cos u))))))))
413 (put 'calcFunc-arcsin\' 'math-derivative-1
414 (function (lambda (u)
416 (math-div 1 (math-normalize
418 (math-sub 1 (math-sqr u)))))))))
420 (put 'calcFunc-arccos\' 'math-derivative-1
421 (function (lambda (u)
423 (math-div -1 (math-normalize
425 (math-sub 1 (math-sqr u)))))))))
427 (put 'calcFunc-arctan\' 'math-derivative-1
428 (function (lambda (u) (math-from-radians-2
429 (math-div 1 (math-add 1 (math-sqr u)))))))
431 (put 'calcFunc-sinh\' 'math-derivative-1
432 (function (lambda (u) (math-normalize (list 'calcFunc-cosh u)))))
434 (put 'calcFunc-cosh\' 'math-derivative-1
435 (function (lambda (u) (math-normalize (list 'calcFunc-sinh u)))))
437 (put 'calcFunc-tanh\' 'math-derivative-1
438 (function (lambda (u) (math-div 1 (math-sqr
440 (list 'calcFunc-cosh u)))))))
442 (put 'calcFunc-arcsinh\' 'math-derivative-1
443 (function (lambda (u)
444 (math-div 1 (math-normalize
446 (math-add (math-sqr u) 1)))))))
448 (put 'calcFunc-arccosh\' 'math-derivative-1
449 (function (lambda (u)
450 (math-div 1 (math-normalize
452 (math-add (math-sqr u) -1)))))))
454 (put 'calcFunc-arctanh\' 'math-derivative-1
455 (function (lambda (u) (math-div 1 (math-sub 1 (math-sqr u))))))
457 (put 'calcFunc-bern\'2 'math-derivative-2
458 (function (lambda (n x)
459 (math-mul n (list 'calcFunc-bern (math-add n -1) x)))))
461 (put 'calcFunc-euler\'2 'math-derivative-2
462 (function (lambda (n x)
463 (math-mul n (list 'calcFunc-euler (math-add n -1) x)))))
465 (put 'calcFunc-gammag\'2 'math-derivative-2
466 (function (lambda (a x) (math-deriv-gamma a x 1))))
468 (put 'calcFunc-gammaG\'2 'math-derivative-2
469 (function (lambda (a x) (math-deriv-gamma a x -1))))
471 (put 'calcFunc-gammaP\'2 'math-derivative-2
472 (function (lambda (a x) (math-deriv-gamma a x
475 (list 'calcFunc-gamma
478 (put 'calcFunc-gammaQ\'2 'math-derivative-2
479 (function (lambda (a x) (math-deriv-gamma a x
482 (list 'calcFunc-gamma
485 (defun math-deriv-gamma (a x scale)
487 (math-mul (math-pow x (math-add a -1))
488 (list 'calcFunc-exp (math-neg x)))))
490 (put 'calcFunc-betaB\' 'math-derivative-3
491 (function (lambda (x a b) (math-deriv-beta x a b 1))))
493 (put 'calcFunc-betaI\' 'math-derivative-3
494 (function (lambda (x a b) (math-deriv-beta x a b
496 1 (list 'calcFunc-beta
499 (defun math-deriv-beta (x a b scale)
500 (math-mul (math-mul (math-pow x (math-add a -1))
501 (math-pow (math-sub 1 x) (math-add b -1)))
504 (put 'calcFunc-erf\' 'math-derivative-1
505 (function (lambda (x) (math-div 2
506 (math-mul (list 'calcFunc-exp
508 (if calc-symbolic-mode
513 (put 'calcFunc-erfc\' 'math-derivative-1
514 (function (lambda (x) (math-div -2
515 (math-mul (list 'calcFunc-exp
517 (if calc-symbolic-mode
522 (put 'calcFunc-besJ\'2 'math-derivative-2
523 (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besJ
531 (put 'calcFunc-besY\'2 'math-derivative-2
532 (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besY
540 (put 'calcFunc-sum 'math-derivative-n
543 (if (math-expr-contains (cons 'vec (cdr (cdr expr))) deriv-var)
544 (throw 'math-deriv nil)
546 (cons (math-derivative (nth 1 expr))
547 (cdr (cdr expr))))))))
549 (put 'calcFunc-prod 'math-derivative-n
552 (if (math-expr-contains (cons 'vec (cdr (cdr expr))) deriv-var)
553 (throw 'math-deriv nil)
556 (cons (math-div (math-derivative (nth 1 expr))
558 (cdr (cdr expr)))))))))
560 (put 'calcFunc-integ 'math-derivative-n
563 (if (= (length expr) 3)
564 (if (equal (nth 2 expr) deriv-var)
567 (list 'calcFunc-integ
568 (math-derivative (nth 1 expr))
570 (if (= (length expr) 5)
571 (let ((lower (math-expr-subst (nth 1 expr) (nth 2 expr)
573 (upper (math-expr-subst (nth 1 expr) (nth 2 expr)
575 (math-add (math-sub (math-mul upper
576 (math-derivative (nth 4 expr)))
578 (math-derivative (nth 3 expr))))
579 (if (equal (nth 2 expr) deriv-var)
582 (list 'calcFunc-integ
583 (math-derivative (nth 1 expr)) (nth 2 expr)
584 (nth 3 expr) (nth 4 expr)))))))))))
586 (put 'calcFunc-if 'math-derivative-n
589 (and (= (length expr) 4)
590 (list 'calcFunc-if (nth 1 expr)
591 (math-derivative (nth 2 expr))
592 (math-derivative (nth 3 expr)))))))
594 (put 'calcFunc-subscr 'math-derivative-n
597 (and (= (length expr) 3)
598 (list 'calcFunc-subscr (nth 1 expr)
599 (math-derivative (nth 2 expr)))))))
602 (defvar math-integ-var '(var X ---))
603 (defvar math-integ-var-2 '(var Y ---))
604 (defvar math-integ-vars (list 'f math-integ-var math-integ-var-2))
605 (defvar math-integ-var-list (list math-integ-var))
606 (defvar math-integ-var-list-list (list math-integ-var-list))
608 (defmacro math-tracing-integral (&rest parts)
611 (list 'save-excursion
612 '(set-buffer trace-buffer)
613 '(goto-char (point-max))
616 '(insert (make-string (- math-integral-limit
617 math-integ-level) 32)
618 (format "%2d " math-integ-depth)
619 (make-string math-integ-level 32)))
620 ;;(list 'condition-case 'err
622 ;; '(error (insert (prin1-to-string err))))
625 ;;; The following wrapper caches results and avoids infinite recursion.
626 ;;; Each cache entry is: ( A B ) Integral of A is B;
627 ;;; ( A N ) Integral of A failed at level N;
628 ;;; ( A busy ) Currently working on integral of A;
629 ;;; ( A parts ) Currently working, integ-by-parts;
630 ;;; ( A parts2 ) Currently working, integ-by-parts;
631 ;;; ( A cancelled ) Ignore this cache entry;
632 ;;; ( A [B] ) Same result as for cur-record = B.
633 (defun math-integral (expr &optional simplify same-as-above)
634 (let* ((simp cur-record)
635 (cur-record (assoc expr math-integral-cache))
636 (math-integ-depth (1+ math-integ-depth))
638 (math-tracing-integral "Integrating "
639 (math-format-value expr 1000)
643 (math-tracing-integral "Found "
644 (math-format-value (nth 1 cur-record) 1000))
645 (and (consp (nth 1 cur-record))
646 (math-replace-integral-parts cur-record))
647 (math-tracing-integral " => "
648 (math-format-value (nth 1 cur-record) 1000)
651 (not (eq (nth 1 cur-record) 'cancelled))
652 (or (not (integerp (nth 1 cur-record)))
653 (>= (nth 1 cur-record) math-integ-level)))
654 (and (math-integral-contains-parts expr)
660 (let (math-integ-msg)
661 (if (eq calc-display-working-message 'lots)
663 (calc-set-command-flag 'clear-message)
664 (setq math-integ-msg (format
665 "Working... Integrating %s"
666 (math-format-flat-expr expr 0)))
667 (message math-integ-msg)))
669 (setcar (cdr cur-record)
670 (if same-as-above (vector simp) 'busy))
672 (list expr (if same-as-above (vector simp) 'busy))
673 math-integral-cache (cons cur-record
674 math-integral-cache)))
675 (if (eq simplify 'yes)
677 (math-tracing-integral "Simplifying...")
678 (setq simp (math-simplify expr))
679 (setq val (if (equal simp expr)
681 (math-tracing-integral " no change\n")
682 (math-do-integral expr))
683 (math-tracing-integral " simplified\n")
684 (math-integral simp 'no t))))
685 (or (setq val (math-do-integral expr))
687 (let ((simp (math-simplify expr)))
688 (or (equal simp expr)
690 (math-tracing-integral "Trying again after "
691 "simplification...\n")
692 (setq val (math-integral simp 'no t))))))))
693 (if (eq calc-display-working-message 'lots)
694 (message math-integ-msg)))
695 (setcar (cdr cur-record) (or val
696 (if (or math-enable-subst
697 (not math-any-substs))
700 (setq val cur-record)
701 (while (vectorp (nth 1 val))
702 (setq val (aref (nth 1 val) 0)))
703 (setq val (if (memq (nth 1 val) '(parts parts2))
705 (setcar (cdr val) 'parts2)
706 (list 'var 'PARTS val))
707 (and (consp (nth 1 val))
709 (math-tracing-integral "Integral of "
710 (math-format-value expr 1000)
712 (math-format-value val 1000)
715 (defvar math-integral-cache nil)
716 (defvar math-integral-cache-state nil)
718 (defun math-integral-contains-parts (expr)
719 (if (Math-primp expr)
720 (and (eq (car-safe expr) 'var)
721 (eq (nth 1 expr) 'PARTS)
722 (listp (nth 2 expr)))
723 (while (and (setq expr (cdr expr))
724 (not (math-integral-contains-parts (car expr)))))
727 (defun math-replace-integral-parts (expr)
728 (or (Math-primp expr)
729 (while (setq expr (cdr expr))
730 (and (consp (car expr))
731 (if (eq (car (car expr)) 'var)
732 (and (eq (nth 1 (car expr)) 'PARTS)
733 (consp (nth 2 (car expr)))
734 (if (listp (nth 1 (nth 2 (car expr))))
736 (setcar expr (nth 1 (nth 2 (car expr))))
737 (math-replace-integral-parts (cons 'foo expr)))
738 (setcar (cdr cur-record) 'cancelled)))
739 (math-replace-integral-parts (car expr)))))))
741 (defun math-do-integral (expr)
743 (or (cond ((not (math-expr-contains expr math-integ-var))
744 (math-mul expr math-integ-var))
745 ((equal expr math-integ-var)
746 (math-div (math-sqr expr) 2))
748 (and (setq t1 (math-integral (nth 1 expr)))
749 (setq t2 (math-integral (nth 2 expr)))
752 (and (setq t1 (math-integral (nth 1 expr)))
753 (setq t2 (math-integral (nth 2 expr)))
755 ((eq (car expr) 'neg)
756 (and (setq t1 (math-integral (nth 1 expr)))
759 (cond ((not (math-expr-contains (nth 1 expr) math-integ-var))
760 (and (setq t1 (math-integral (nth 2 expr)))
761 (math-mul (nth 1 expr) t1)))
762 ((not (math-expr-contains (nth 2 expr) math-integ-var))
763 (and (setq t1 (math-integral (nth 1 expr)))
764 (math-mul t1 (nth 2 expr))))
765 ((memq (car-safe (nth 1 expr)) '(+ -))
766 (math-integral (list (car (nth 1 expr))
767 (math-mul (nth 1 (nth 1 expr))
769 (math-mul (nth 2 (nth 1 expr))
772 ((memq (car-safe (nth 2 expr)) '(+ -))
773 (math-integral (list (car (nth 2 expr))
774 (math-mul (nth 1 (nth 2 expr))
776 (math-mul (nth 2 (nth 2 expr))
780 (cond ((and (not (math-expr-contains (nth 1 expr)
782 (not (math-equal-int (nth 1 expr) 1)))
783 (and (setq t1 (math-integral (math-div 1 (nth 2 expr))))
784 (math-mul (nth 1 expr) t1)))
785 ((not (math-expr-contains (nth 2 expr) math-integ-var))
786 (and (setq t1 (math-integral (nth 1 expr)))
787 (math-div t1 (nth 2 expr))))
788 ((and (eq (car-safe (nth 1 expr)) '*)
789 (not (math-expr-contains (nth 1 (nth 1 expr))
791 (and (setq t1 (math-integral
792 (math-div (nth 2 (nth 1 expr))
794 (math-mul t1 (nth 1 (nth 1 expr)))))
795 ((and (eq (car-safe (nth 1 expr)) '*)
796 (not (math-expr-contains (nth 2 (nth 1 expr))
798 (and (setq t1 (math-integral
799 (math-div (nth 1 (nth 1 expr))
801 (math-mul t1 (nth 2 (nth 1 expr)))))
802 ((and (eq (car-safe (nth 2 expr)) '*)
803 (not (math-expr-contains (nth 1 (nth 2 expr))
805 (and (setq t1 (math-integral
806 (math-div (nth 1 expr)
807 (nth 2 (nth 2 expr)))))
808 (math-div t1 (nth 1 (nth 2 expr)))))
809 ((and (eq (car-safe (nth 2 expr)) '*)
810 (not (math-expr-contains (nth 2 (nth 2 expr))
812 (and (setq t1 (math-integral
813 (math-div (nth 1 expr)
814 (nth 1 (nth 2 expr)))))
815 (math-div t1 (nth 2 (nth 2 expr)))))
816 ((eq (car-safe (nth 2 expr)) 'calcFunc-exp)
818 (math-mul (nth 1 expr)
820 (math-neg (nth 1 (nth 2 expr)))))))))
822 (cond ((not (math-expr-contains (nth 1 expr) math-integ-var))
823 (or (and (setq t1 (math-is-polynomial (nth 2 expr)
832 (math-mul (nth 2 expr)
837 ((not (math-expr-contains (nth 2 expr) math-integ-var))
838 (if (and (integerp (nth 2 expr)) (< (nth 2 expr) 0))
840 (list '/ 1 (math-pow (nth 1 expr) (- (nth 2 expr))))
842 (or (and (setq t1 (math-is-polynomial (nth 1 expr)
845 (setq t2 (math-add (nth 2 expr) 1))
846 (math-div (math-pow (nth 1 expr) t2)
847 (math-mul t2 (nth 1 t1))))
848 (and (Math-negp (nth 2 expr))
851 (math-pow (nth 1 expr)
857 ;; Integral of a polynomial.
858 (and (setq t1 (math-is-polynomial expr math-integ-var 20))
862 (if (setq accum (math-add accum
863 (math-div (math-mul (car t1)
872 ;; Try looking it up!
873 (cond ((= (length expr) 2)
874 (and (symbolp (car expr))
875 (setq t1 (get (car expr) 'math-integral))
878 (not (setq t2 (funcall (car t1)
881 (and t2 (math-normalize t2)))))
883 (and (symbolp (car expr))
884 (setq t1 (get (car expr) 'math-integral-2))
887 (not (setq t2 (funcall (car t1)
891 (and t2 (math-normalize t2))))))
893 ;; Integral of a rational function.
894 (and (math-ratpoly-p expr math-integ-var)
895 (setq t1 (calcFunc-apart expr math-integ-var))
896 (not (equal t1 expr))
899 ;; Try user-defined integration rules.
901 (let ((math-old-integ (symbol-function 'calcFunc-integ))
902 (input (list 'calcFunc-integtry expr math-integ-var))
906 (fset 'calcFunc-integ 'math-sub-integration)
907 (setq res (math-rewrite input
908 '(var IntegRules var-IntegRules)
910 (fset 'calcFunc-integ math-old-integ)
911 (and (not (equal res input))
912 (if (setq part (math-expr-calls
913 res '(calcFunc-integsubst)))
914 (and (memq (length part) '(3 4 5))
922 (math-integrate-by-substitution
925 (list 'calcFunc-integfailed
928 (if (not (math-expr-calls res
930 calcFunc-integfailed)))
932 (fset 'calcFunc-integ math-old-integ))))
934 ;; See if the function is a symbolic derivative.
935 (and (string-match "'" (symbol-name (car expr)))
936 (let ((name (symbol-name (car expr)))
937 (p expr) (n 0) (which nil) (bad nil))
938 (while (setq n (1+ n) p (cdr p))
939 (if (equal (car p) math-integ-var)
940 (if which (setq bad t) (setq which n))
941 (if (math-expr-contains (car p) math-integ-var)
944 (let ((prime (if (= which 1) "'" (format "'%d" which))))
945 (and (string-match (concat prime "\\('['0-9]*\\|$\\)")
949 (substring name 0 (match-beginning 0))
950 (substring name (+ (match-beginning 0)
954 ;; Try transformation methods (parts, substitutions).
955 (and (> math-integ-level 0)
956 (math-do-integral-methods expr))
958 ;; Try expanding the function's definition.
959 (let ((res (math-expand-formula expr)))
961 (math-integral res))))))
963 (defun math-sub-integration (expr &rest rest)
964 (or (if (or (not rest)
965 (and (< math-integ-level math-integral-limit)
966 (eq (car rest) math-integ-var)))
968 (let ((res (apply math-old-integ expr rest)))
969 (and (or (= math-integ-level math-integral-limit)
970 (not (math-expr-calls res 'calcFunc-integ)))
972 (list 'calcFunc-integfailed expr)))
974 (defun math-do-integral-methods (expr)
975 (let ((so-far math-integ-var-list-list)
978 ;; Integration by substitution, for various likely sub-expressions.
979 ;; (In first pass, we look only for sub-exprs that are linear in X.)
980 (or (if math-enable-subst
981 (math-integ-try-substitutions expr)
982 (math-integ-try-linear-substitutions expr))
984 ;; If function has sines and cosines, try tan(x/2) substitution.
985 (and (let ((p (setq rat-in (math-expr-rational-in expr))))
987 (memq (car (car p)) '(calcFunc-sin
990 (equal (nth 1 (car p)) math-integ-var))
993 (or (and (math-integ-parts-easy expr)
994 (math-integ-try-parts expr t))
995 (math-integrate-by-good-substitution
996 expr (list 'calcFunc-tan (math-div math-integ-var 2)))))
998 ;; If function has sinh and cosh, try tanh(x/2) substitution.
999 (and (let ((p rat-in))
1001 (memq (car (car p)) '(calcFunc-sinh
1005 (equal (nth 1 (car p)) math-integ-var))
1008 (or (and (math-integ-parts-easy expr)
1009 (math-integ-try-parts expr t))
1010 (math-integrate-by-good-substitution
1011 expr (list 'calcFunc-tanh (math-div math-integ-var 2)))))
1013 ;; If function has square roots, try sin, tan, or sec substitution.
1014 (and (let ((p rat-in))
1017 (or (equal (car p) math-integ-var)
1018 (and (eq (car (car p)) 'calcFunc-sqrt)
1019 (setq t1 (math-is-polynomial
1020 (nth 1 (setq t2 (car p)))
1021 math-integ-var 2)))))
1025 (if (math-guess-if-neg (nth 2 t1))
1026 (let* ((c (math-sqrt (math-neg (nth 2 t1))))
1027 (d (math-div (nth 1 t1) (math-mul -2 c)))
1028 (a (math-sqrt (math-add (car t1) (math-sqr d)))))
1029 (math-integrate-by-good-substitution
1030 expr (list 'calcFunc-arcsin
1032 (math-add (math-mul c math-integ-var) d)
1034 (let* ((c (math-sqrt (nth 2 t1)))
1035 (d (math-div (nth 1 t1) (math-mul 2 c)))
1036 (aa (math-sub (car t1) (math-sqr d))))
1037 (if (and nil (not (and (eq d 0) (eq c 1))))
1038 (math-integrate-by-good-substitution
1039 expr (math-add (math-mul c math-integ-var) d))
1040 (if (math-guess-if-neg aa)
1041 (math-integrate-by-good-substitution
1042 expr (list 'calcFunc-arccosh
1044 (math-add (math-mul c math-integ-var)
1046 (math-sqrt (math-neg aa)))))
1047 (math-integrate-by-good-substitution
1048 expr (list 'calcFunc-arcsinh
1050 (math-add (math-mul c math-integ-var)
1052 (math-sqrt aa))))))))
1053 (math-integrate-by-good-substitution expr t2)) )
1055 ;; Try integration by parts.
1056 (math-integ-try-parts expr)
1061 (defun math-integ-parts-easy (expr)
1062 (cond ((Math-primp expr) t)
1063 ((memq (car expr) '(+ - *))
1064 (and (math-integ-parts-easy (nth 1 expr))
1065 (math-integ-parts-easy (nth 2 expr))))
1067 (and (math-integ-parts-easy (nth 1 expr))
1068 (math-atomic-factorp (nth 2 expr))))
1070 (and (natnump (nth 2 expr))
1071 (math-integ-parts-easy (nth 1 expr))))
1072 ((eq (car expr) 'neg)
1073 (math-integ-parts-easy (nth 1 expr)))
1076 (defun math-integ-try-parts (expr &optional math-good-parts)
1077 ;; Integration by parts:
1078 ;; integ(f(x) g(x),x) = f(x) h(x) - integ(h(x) f'(x),x)
1079 ;; where h(x) = integ(g(x),x).
1080 (or (let ((exp (calcFunc-expand expr)))
1081 (and (not (equal exp expr))
1082 (math-integral exp)))
1083 (and (eq (car expr) '*)
1084 (let ((first-bad (or (math-polynomial-p (nth 1 expr)
1086 (equal (nth 2 expr) math-prev-parts-v))))
1087 (or (and first-bad ; so try this one first
1088 (math-integrate-by-parts (nth 1 expr) (nth 2 expr)))
1089 (math-integrate-by-parts (nth 2 expr) (nth 1 expr))
1090 (and (not first-bad)
1091 (math-integrate-by-parts (nth 1 expr) (nth 2 expr))))))
1092 (and (eq (car expr) '/)
1093 (math-expr-contains (nth 1 expr) math-integ-var)
1094 (let ((recip (math-div 1 (nth 2 expr))))
1095 (or (math-integrate-by-parts (nth 1 expr) recip)
1096 (math-integrate-by-parts recip (nth 1 expr)))))
1097 (and (eq (car expr) '^)
1098 (math-integrate-by-parts (math-pow (nth 1 expr)
1099 (math-sub (nth 2 expr) 1))
1102 (defun math-integrate-by-parts (u vprime)
1103 (let ((math-integ-level (if (or math-good-parts
1104 (math-polynomial-p u math-integ-var))
1106 (1- math-integ-level)))
1107 (math-doing-parts t)
1109 (and (>= math-integ-level 0)
1112 (setcar (cdr cur-record) 'parts)
1113 (math-tracing-integral "Integrating by parts, u = "
1114 (math-format-value u 1000)
1116 (math-format-value vprime 1000)
1118 (and (setq v (math-integral vprime))
1119 (setq temp (calcFunc-deriv u math-integ-var nil t))
1120 (setq temp (let ((math-prev-parts-v v))
1121 (math-integral (math-mul v temp) 'yes)))
1122 (setq temp (math-sub (math-mul u v) temp))
1123 (if (eq (nth 1 cur-record) 'parts)
1124 (calcFunc-expand temp)
1125 (setq v (list 'var 'PARTS cur-record)
1126 var-thing (list 'vec (math-sub v temp) v)
1127 temp (let (calc-next-why)
1128 (math-solve-for (math-sub v temp) 0 v nil)))
1129 (and temp (not (integerp temp))
1130 (math-simplify-extended temp)))))
1131 (setcar (cdr cur-record) 'busy)))))
1133 ;;; This tries two different formulations, hoping the algebraic simplifier
1134 ;;; will be strong enough to handle at least one.
1135 (defun math-integrate-by-substitution (expr u &optional user uinv uinvprime)
1136 (and (> math-integ-level 0)
1137 (let ((math-integ-level (max (- math-integ-level 2) 0)))
1138 (math-integrate-by-good-substitution expr u user uinv uinvprime))))
1140 (defun math-integrate-by-good-substitution (expr u &optional user
1142 (let ((math-living-dangerously t)
1144 (and (setq uinv (if uinv
1145 (math-expr-subst uinv math-integ-var
1147 (let (calc-next-why)
1150 math-integ-var nil))))
1152 (math-tracing-integral "Integrating by substitution, u = "
1153 (math-format-value u 1000)
1155 (or (and (setq deriv (calcFunc-deriv u
1158 (setq temp (math-integral (math-expr-subst
1161 (math-div expr deriv)
1169 (and (setq deriv (or uinvprime
1170 (calcFunc-deriv uinv
1174 (setq temp (math-integral (math-mul
1187 (math-simplify-extended
1188 (math-expr-subst temp math-integ-var u)))))
1190 ;;; Look for substitutions of the form u = a x + b.
1191 (defun math-integ-try-linear-substitutions (sub-expr)
1192 (and (not (Math-primp sub-expr))
1193 (or (and (not (memq (car sub-expr) '(+ - * / neg)))
1194 (not (and (eq (car sub-expr) '^)
1195 (integerp (nth 2 sub-expr))))
1196 (math-expr-contains sub-expr math-integ-var)
1198 (while (and (setq sub-expr (cdr sub-expr))
1199 (or (not (math-linear-in (car sub-expr)
1201 (assoc (car sub-expr) so-far)
1203 (setq so-far (cons (list (car sub-expr))
1206 (math-integrate-by-substitution
1207 expr (car sub-expr))))))))
1210 (while (and (setq sub-expr (cdr sub-expr))
1211 (not (setq res (math-integ-try-linear-substitutions
1215 ;;; Recursively try different substitutions based on various sub-expressions.
1216 (defun math-integ-try-substitutions (sub-expr &optional allow-rat)
1217 (and (not (Math-primp sub-expr))
1218 (not (assoc sub-expr so-far))
1219 (math-expr-contains sub-expr math-integ-var)
1220 (or (and (if (and (not (memq (car sub-expr) '(+ - * / neg)))
1221 (not (and (eq (car sub-expr) '^)
1222 (integerp (nth 2 sub-expr)))))
1224 (prog1 allow-rat (setq allow-rat nil)))
1225 (not (eq sub-expr expr))
1226 (or (math-integrate-by-substitution expr sub-expr)
1227 (and (eq (car sub-expr) '^)
1228 (integerp (nth 2 sub-expr))
1229 (< (nth 2 sub-expr) 0)
1230 (math-integ-try-substitutions
1231 (math-pow (nth 1 sub-expr) (- (nth 2 sub-expr)))
1234 (setq so-far (cons (list sub-expr) so-far))
1235 (while (and (setq sub-expr (cdr sub-expr))
1236 (not (setq res (math-integ-try-substitutions
1237 (car sub-expr) allow-rat)))))
1240 (defun math-expr-rational-in (expr)
1242 (math-expr-rational-in-rec expr)
1243 (mapcar 'car parts)))
1245 (defun math-expr-rational-in-rec (expr)
1246 (cond ((Math-primp expr)
1247 (and (equal expr math-integ-var)
1248 (not (assoc expr parts))
1249 (setq parts (cons (list expr) parts))))
1250 ((or (memq (car expr) '(+ - * / neg))
1251 (and (eq (car expr) '^) (integerp (nth 2 expr))))
1252 (math-expr-rational-in-rec (nth 1 expr))
1253 (and (nth 2 expr) (math-expr-rational-in-rec (nth 2 expr))))
1254 ((and (eq (car expr) '^)
1255 (eq (math-quarter-integer (nth 2 expr)) 2))
1256 (math-expr-rational-in-rec (list 'calcFunc-sqrt (nth 1 expr))))
1258 (and (not (assoc expr parts))
1259 (math-expr-contains expr math-integ-var)
1260 (setq parts (cons (list expr) parts))))))
1262 (defun math-expr-calls (expr funcs &optional arg-contains)
1264 (if (or (memq (car expr) funcs)
1265 (and (eq (car expr) '^) (eq (car funcs) 'calcFunc-sqrt)
1266 (eq (math-quarter-integer (nth 2 expr)) 2)))
1267 (and (or (not arg-contains)
1268 (math-expr-contains expr arg-contains))
1270 (and (not (Math-primp expr))
1272 (while (and (setq expr (cdr expr))
1273 (not (setq res (math-expr-calls
1274 (car expr) funcs arg-contains)))))
1277 (defun math-fix-const-terms (expr except-vars)
1278 (cond ((not (math-expr-depends expr except-vars)) 0)
1279 ((Math-primp expr) expr)
1281 (math-add (math-fix-const-terms (nth 1 expr) except-vars)
1282 (math-fix-const-terms (nth 2 expr) except-vars)))
1284 (math-sub (math-fix-const-terms (nth 1 expr) except-vars)
1285 (math-fix-const-terms (nth 2 expr) except-vars)))
1288 ;; Command for debugging the Calculator's symbolic integrator.
1289 (defun calc-dump-integral-cache (&optional arg)
1291 (let ((buf (current-buffer)))
1293 (let ((p math-integral-cache)
1295 (display-buffer (get-buffer-create "*Integral Cache*"))
1296 (set-buffer (get-buffer "*Integral Cache*"))
1299 (setq cur-record (car p))
1300 (or arg (math-replace-integral-parts cur-record))
1301 (insert (math-format-flat-expr (car cur-record) 0)
1303 (if (symbolp (nth 1 cur-record))
1304 (concat "(" (symbol-name (nth 1 cur-record)) ")")
1305 (math-format-flat-expr (nth 1 cur-record) 0))
1308 (goto-char (point-min)))
1311 (defun math-try-integral (expr)
1312 (let ((math-integ-level math-integral-limit)
1313 (math-integ-depth 0)
1314 (math-integ-msg "Working...done")
1315 (cur-record nil) ; a technicality
1316 (math-integrating t)
1317 (calc-prefer-frac t)
1318 (calc-symbolic-mode t)
1319 (has-rules (calc-has-rules 'var-IntegRules)))
1320 (or (math-integral expr 'yes)
1321 (and math-any-substs
1322 (setq math-enable-subst t)
1323 (math-integral expr 'yes))
1324 (and (> math-max-integral-limit math-integral-limit)
1325 (setq math-integral-limit math-max-integral-limit
1326 math-integ-level math-integral-limit)
1327 (math-integral expr 'yes)))))
1329 (defun calcFunc-integ (expr var &optional low high)
1331 ;; Do these even if the parts turn out not to be integrable.
1332 ((eq (car-safe expr) '+)
1333 (math-add (calcFunc-integ (nth 1 expr) var low high)
1334 (calcFunc-integ (nth 2 expr) var low high)))
1335 ((eq (car-safe expr) '-)
1336 (math-sub (calcFunc-integ (nth 1 expr) var low high)
1337 (calcFunc-integ (nth 2 expr) var low high)))
1338 ((eq (car-safe expr) 'neg)
1339 (math-neg (calcFunc-integ (nth 1 expr) var low high)))
1340 ((and (eq (car-safe expr) '*)
1341 (not (math-expr-contains (nth 1 expr) var)))
1342 (math-mul (nth 1 expr) (calcFunc-integ (nth 2 expr) var low high)))
1343 ((and (eq (car-safe expr) '*)
1344 (not (math-expr-contains (nth 2 expr) var)))
1345 (math-mul (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr)))
1346 ((and (eq (car-safe expr) '/)
1347 (not (math-expr-contains (nth 1 expr) var))
1348 (not (math-equal-int (nth 1 expr) 1)))
1349 (math-mul (nth 1 expr)
1350 (calcFunc-integ (math-div 1 (nth 2 expr)) var low high)))
1351 ((and (eq (car-safe expr) '/)
1352 (not (math-expr-contains (nth 2 expr) var)))
1353 (math-div (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr)))
1354 ((and (eq (car-safe expr) '/)
1355 (eq (car-safe (nth 1 expr)) '*)
1356 (not (math-expr-contains (nth 1 (nth 1 expr)) var)))
1357 (math-mul (nth 1 (nth 1 expr))
1358 (calcFunc-integ (math-div (nth 2 (nth 1 expr)) (nth 2 expr))
1360 ((and (eq (car-safe expr) '/)
1361 (eq (car-safe (nth 1 expr)) '*)
1362 (not (math-expr-contains (nth 2 (nth 1 expr)) var)))
1363 (math-mul (nth 2 (nth 1 expr))
1364 (calcFunc-integ (math-div (nth 1 (nth 1 expr)) (nth 2 expr))
1366 ((and (eq (car-safe expr) '/)
1367 (eq (car-safe (nth 2 expr)) '*)
1368 (not (math-expr-contains (nth 1 (nth 2 expr)) var)))
1369 (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 2 (nth 2 expr)))
1371 (nth 1 (nth 2 expr))))
1372 ((and (eq (car-safe expr) '/)
1373 (eq (car-safe (nth 2 expr)) '*)
1374 (not (math-expr-contains (nth 2 (nth 2 expr)) var)))
1375 (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 1 (nth 2 expr)))
1377 (nth 2 (nth 2 expr))))
1378 ((eq (car-safe expr) 'vec)
1379 (cons 'vec (mapcar (function (lambda (x) (calcFunc-integ x var low high)))
1382 (let ((state (list calc-angle-mode
1383 ;;calc-symbolic-mode
1386 (calc-var-value 'var-IntegRules)
1387 (calc-var-value 'var-IntegSimpRules))))
1388 (or (equal state math-integral-cache-state)
1389 (setq math-integral-cache-state state
1390 math-integral-cache nil)))
1391 (let* ((math-max-integral-limit (or (and (boundp 'var-IntegLimit)
1392 (natnump var-IntegLimit)
1395 (math-integral-limit 1)
1396 (sexpr (math-expr-subst expr var math-integ-var))
1397 (trace-buffer (get-buffer "*Trace*"))
1398 (calc-language (if (eq calc-language 'big) nil calc-language))
1400 (math-enable-subst nil)
1401 (math-prev-parts-v nil)
1402 (math-doing-parts nil)
1403 (math-good-parts nil)
1406 (let ((calcbuf (current-buffer))
1407 (calcwin (selected-window)))
1410 (if (get-buffer-window trace-buffer)
1411 (select-window (get-buffer-window trace-buffer)))
1412 (set-buffer trace-buffer)
1413 (goto-char (point-max))
1414 (or (assq 'scroll-stop (buffer-local-variables))
1416 (make-local-variable 'scroll-step)
1417 (setq scroll-step 3)))
1419 (set-buffer calcbuf)
1420 (math-try-integral sexpr))
1421 (select-window calcwin)
1422 (set-buffer calcbuf)))
1423 (math-try-integral sexpr))))
1426 (if (calc-has-rules 'var-IntegAfterRules)
1427 (setq res (math-rewrite res '(var IntegAfterRules
1428 var-IntegAfterRules))))
1431 (math-sub (math-expr-subst res math-integ-var high)
1432 (math-expr-subst res math-integ-var low))
1433 (setq res (math-fix-const-terms res math-integ-vars))
1435 (math-expr-subst res math-integ-var low)
1436 (math-expr-subst res math-integ-var var)))))
1437 (append (list 'calcFunc-integ expr var)
1438 (and low (list low))
1439 (and high (list high))))))))
1442 (math-defintegral calcFunc-inv
1443 (math-integral (math-div 1 u)))
1445 (math-defintegral calcFunc-conj
1446 (let ((int (math-integral u)))
1448 (list 'calcFunc-conj int))))
1450 (math-defintegral calcFunc-deg
1451 (let ((int (math-integral u)))
1453 (list 'calcFunc-deg int))))
1455 (math-defintegral calcFunc-rad
1456 (let ((int (math-integral u)))
1458 (list 'calcFunc-rad int))))
1460 (math-defintegral calcFunc-re
1461 (let ((int (math-integral u)))
1463 (list 'calcFunc-re int))))
1465 (math-defintegral calcFunc-im
1466 (let ((int (math-integral u)))
1468 (list 'calcFunc-im int))))
1470 (math-defintegral calcFunc-sqrt
1471 (and (equal u math-integ-var)
1472 (math-mul '(frac 2 3)
1473 (list 'calcFunc-sqrt (math-pow u 3)))))
1475 (math-defintegral calcFunc-exp
1476 (or (and (equal u math-integ-var)
1477 (list 'calcFunc-exp u))
1478 (let ((p (math-is-polynomial u math-integ-var 2)))
1480 (let ((sqa (math-sqrt (math-neg (nth 2 p)))))
1483 (math-mul (math-div (list 'calcFunc-sqrt '(var pi var-pi))
1487 (math-div (math-sub (math-mul (car p)
1490 (math-sqr (nth 1 p))
1494 (math-sub (math-mul sqa math-integ-var)
1495 (math-div (nth 1 p) (math-mul 2 sqa)))))
1498 (math-defintegral calcFunc-ln
1499 (or (and (equal u math-integ-var)
1500 (math-sub (math-mul u (list 'calcFunc-ln u)) u))
1501 (and (eq (car u) '*)
1502 (math-integral (math-add (list 'calcFunc-ln (nth 1 u))
1503 (list 'calcFunc-ln (nth 2 u)))))
1504 (and (eq (car u) '/)
1505 (math-integral (math-sub (list 'calcFunc-ln (nth 1 u))
1506 (list 'calcFunc-ln (nth 2 u)))))
1507 (and (eq (car u) '^)
1508 (math-integral (math-mul (nth 2 u)
1509 (list 'calcFunc-ln (nth 1 u)))))))
1511 (math-defintegral calcFunc-log10
1512 (and (equal u math-integ-var)
1513 (math-sub (math-mul u (list 'calcFunc-ln u))
1514 (math-div u (list 'calcFunc-ln 10)))))
1516 (math-defintegral-2 calcFunc-log
1517 (math-integral (math-div (list 'calcFunc-ln u)
1518 (list 'calcFunc-ln v))))
1520 (math-defintegral calcFunc-sin
1521 (or (and (equal u math-integ-var)
1522 (math-neg (math-from-radians-2 (list 'calcFunc-cos u))))
1523 (and (nth 2 (math-is-polynomial u math-integ-var 2))
1524 (math-integral (math-to-exponentials (list 'calcFunc-sin u))))))
1526 (math-defintegral calcFunc-cos
1527 (or (and (equal u math-integ-var)
1528 (math-from-radians-2 (list 'calcFunc-sin u)))
1529 (and (nth 2 (math-is-polynomial u math-integ-var 2))
1530 (math-integral (math-to-exponentials (list 'calcFunc-cos u))))))
1532 (math-defintegral calcFunc-tan
1533 (and (equal u math-integ-var)
1534 (math-neg (math-from-radians-2
1535 (list 'calcFunc-ln (list 'calcFunc-cos u))))))
1537 (math-defintegral calcFunc-arcsin
1538 (and (equal u math-integ-var)
1539 (math-add (math-mul u (list 'calcFunc-arcsin u))
1540 (math-from-radians-2
1541 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))))
1543 (math-defintegral calcFunc-arccos
1544 (and (equal u math-integ-var)
1545 (math-sub (math-mul u (list 'calcFunc-arccos u))
1546 (math-from-radians-2
1547 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))))
1549 (math-defintegral calcFunc-arctan
1550 (and (equal u math-integ-var)
1551 (math-sub (math-mul u (list 'calcFunc-arctan u))
1552 (math-from-radians-2
1553 (math-div (list 'calcFunc-ln (math-add 1 (math-sqr u)))
1556 (math-defintegral calcFunc-sinh
1557 (and (equal u math-integ-var)
1558 (list 'calcFunc-cosh u)))
1560 (math-defintegral calcFunc-cosh
1561 (and (equal u math-integ-var)
1562 (list 'calcFunc-sinh u)))
1564 (math-defintegral calcFunc-tanh
1565 (and (equal u math-integ-var)
1566 (list 'calcFunc-ln (list 'calcFunc-cosh u))))
1568 (math-defintegral calcFunc-arcsinh
1569 (and (equal u math-integ-var)
1570 (math-sub (math-mul u (list 'calcFunc-arcsinh u))
1571 (list 'calcFunc-sqrt (math-add (math-sqr u) 1)))))
1573 (math-defintegral calcFunc-arccosh
1574 (and (equal u math-integ-var)
1575 (math-sub (math-mul u (list 'calcFunc-arccosh u))
1576 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))
1578 (math-defintegral calcFunc-arctanh
1579 (and (equal u math-integ-var)
1580 (math-sub (math-mul u (list 'calcFunc-arctan u))
1581 (math-div (list 'calcFunc-ln
1582 (math-add 1 (math-sqr u)))
1585 ;;; (Ax + B) / (ax^2 + bx + c)^n forms.
1586 (math-defintegral-2 /
1587 (math-integral-rational-funcs u v))
1589 (defun math-integral-rational-funcs (u v)
1590 (let ((pu (math-is-polynomial u math-integ-var 1))
1594 (if (and (eq (car-safe v) '^) (natnump (nth 2 v)))
1595 (setq vpow (nth 2 v)
1597 (and (setq pv (math-is-polynomial v math-integ-var 2))
1598 (let ((int (math-mul-thru
1600 (math-integral-q02 (car pv) (nth 1 pv)
1601 (nth 2 pv) v vpow))))
1603 (setq int (math-add int
1608 (nth 2 pv) v vpow)))))
1611 (defun math-integral-q12 (a b c v vpow)
1615 (math-sub (math-div math-integ-var b)
1616 (math-mul (math-div a (math-sqr b))
1617 (list 'calcFunc-ln v))))
1619 (math-div (math-add (list 'calcFunc-ln v)
1623 (let ((nm1 (math-sub vpow 1))
1624 (nm2 (math-sub vpow 2)))
1626 (math-div a (math-mul nm1 (math-pow v nm1)))
1627 (math-div 1 (math-mul nm2 (math-pow v nm2))))
1630 (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b))))
1631 (let ((part (math-div b (math-mul 2 c))))
1632 (math-mul-thru (math-pow c vpow)
1633 (math-integral-q12 part 1 nil
1634 (math-add math-integ-var part)
1637 (and (math-ratp q) (math-negp q)
1638 (let ((calc-symbolic-mode t))
1639 (math-ratp (math-sqrt (math-neg q))))
1640 (throw 'int-rat nil)) ; should have used calcFunc-apart first
1641 (math-sub (math-div (list 'calcFunc-ln v) (math-mul 2 c))
1642 (math-mul-thru (math-div b (math-mul 2 c))
1643 (math-integral-q02 a b c v 1))))
1645 (let ((n (1- vpow)))
1646 (math-sub (math-neg (math-div
1647 (math-add (math-mul b math-integ-var)
1649 (math-mul n (math-mul q (math-pow v n)))))
1650 (math-mul-thru (math-div (math-mul b (1- (* 2 n)))
1652 (math-integral-q02 a b c v n))))))))
1654 (defun math-integral-q02 (a b c v vpow)
1658 (math-div (list 'calcFunc-ln v) b))
1660 (math-div (math-pow v (- 1 vpow))
1661 (math-mul (- 1 vpow) b)))))
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-q02 part 1 nil
1667 (math-add math-integ-var part)
1670 (setq part (math-add (math-mul 2 (math-mul c math-integ-var)) b))
1672 (let ((n (1- vpow)))
1673 (math-add (math-div part (math-mul n (math-mul q (math-pow v n))))
1674 (math-mul-thru (math-div (math-mul (- (* 4 n) 2) c)
1676 (math-integral-q02 a b c v n)))))
1677 ((math-guess-if-neg q)
1678 (setq rq (list 'calcFunc-sqrt (math-neg q)))
1679 ;;(math-div-thru (list 'calcFunc-ln
1680 ;; (math-div (math-sub part rq)
1681 ;; (math-add part rq)))
1683 (math-div (math-mul -2 (list 'calcFunc-arctanh
1684 (math-div part rq)))
1687 (setq rq (list 'calcFunc-sqrt q))
1688 (math-div (math-mul 2 (math-to-radians-2
1689 (list 'calcFunc-arctan
1690 (math-div part rq))))
1694 (math-defintegral calcFunc-erf
1695 (and (equal u math-integ-var)
1696 (math-add (math-mul u (list 'calcFunc-erf u))
1697 (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u))
1698 (list 'calcFunc-sqrt
1699 '(var pi var-pi)))))))
1701 (math-defintegral calcFunc-erfc
1702 (and (equal u math-integ-var)
1703 (math-sub (math-mul u (list 'calcFunc-erfc u))
1704 (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u))
1705 (list 'calcFunc-sqrt
1706 '(var pi var-pi)))))))
1711 (defvar math-tabulate-initial nil)
1712 (defvar math-tabulate-function nil)
1713 (defun calcFunc-table (expr var &optional low high step)
1714 (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
1715 (or high (setq high low low 1))
1716 (and (or (math-infinitep low) (math-infinitep high))
1718 (math-scan-for-limits expr))
1719 (and step (math-zerop step) (math-reject-arg step 'nonzerop))
1720 (let ((known (+ (if (Math-objectp low) 1 0)
1721 (if (Math-objectp high) 1 0)
1722 (if (or (null step) (Math-objectp step)) 1 0)))
1723 (count '(var inf var-inf))
1725 (or (= known 2) ; handy optimization
1726 (equal high '(var inf var-inf))
1728 (setq count (math-div (math-sub high low) (or step 1)))
1729 (or (Math-objectp count)
1730 (setq count (math-simplify count)))
1731 (if (Math-messy-integerp count)
1732 (setq count (math-trunc count)))))
1733 (if (Math-negp count)
1735 (if (integerp count)
1736 (let ((var-DUMMY nil)
1737 (vec math-tabulate-initial)
1738 (math-working-step-2 (1+ count))
1739 (math-working-step 0))
1740 (setq expr (math-evaluate-expr
1741 (math-expr-subst expr var '(var DUMMY var-DUMMY))))
1743 (setq math-working-step (1+ math-working-step)
1745 vec (cond ((eq math-tabulate-function 'calcFunc-sum)
1746 (math-add vec (math-evaluate-expr expr)))
1747 ((eq math-tabulate-function 'calcFunc-prod)
1748 (math-mul vec (math-evaluate-expr expr)))
1750 (cons (math-evaluate-expr expr) vec)))
1751 low (math-add low (or step 1))
1753 (if math-tabulate-function
1755 (cons 'vec (nreverse vec))))
1756 (if (Math-integerp count)
1757 (calc-record-why 'fixnump high)
1758 (if (Math-num-integerp low)
1759 (if (Math-num-integerp high)
1760 (calc-record-why 'integerp step)
1761 (calc-record-why 'integerp high))
1762 (calc-record-why 'integerp low)))
1763 (append (list (or math-tabulate-function 'calcFunc-table)
1765 (and (not (and (equal low '(neg (var inf var-inf)))
1766 (equal high '(var inf var-inf))))
1768 (and step (list step))))))
1770 (defun math-scan-for-limits (x)
1771 (cond ((Math-primp x))
1772 ((and (eq (car x) 'calcFunc-subscr)
1773 (Math-vectorp (nth 1 x))
1774 (math-expr-contains (nth 2 x) var))
1775 (let* ((calc-next-why nil)
1776 (low-val (math-solve-for (nth 2 x) 1 var nil))
1777 (high-val (math-solve-for (nth 2 x) (1- (length (nth 1 x)))
1780 (and low-val (math-realp low-val)
1781 high-val (math-realp high-val))
1782 (and (Math-lessp high-val low-val)
1783 (setq temp low-val low-val high-val high-val temp))
1784 (setq low (math-max low (math-ceiling low-val))
1785 high (math-min high (math-floor high-val)))))
1787 (while (setq x (cdr x))
1788 (math-scan-for-limits (car x))))))
1791 (defvar math-disable-sums nil)
1792 (defun calcFunc-sum (expr var &optional low high step)
1793 (if math-disable-sums (math-reject-arg))
1794 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2)))
1795 (math-sum-rec expr var low high step)))
1796 (math-disable-sums t))
1797 (math-normalize res)))
1799 (defun math-sum-rec (expr var &optional low high step)
1800 (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
1801 (and low (not high) (setq high low low 1))
1805 ((not (math-expr-contains expr var))
1806 (math-mul expr (math-add (math-div (math-sub high low) (or step 1))
1808 ((and step (not (math-equal-int step 1)))
1809 (if (math-negp step)
1810 (math-sum-rec expr var high low (math-neg step))
1811 (let ((lo (math-simplify (math-div low step))))
1812 (if (math-known-num-integerp lo)
1813 (math-sum-rec (math-normalize
1814 (math-expr-subst expr var
1815 (math-mul step var)))
1816 var lo (math-simplify (math-div high step)))
1817 (math-sum-rec (math-normalize
1818 (math-expr-subst expr var
1819 (math-add (math-mul step var)
1822 (math-simplify (math-div (math-sub high low)
1824 ((memq (setq t1 (math-compare low high)) '(0 1))
1826 (math-expr-subst expr var low)
1828 ((setq t1 (math-is-polynomial expr var 20))
1832 (setq poly (math-poly-mix poly 1
1833 (math-sum-integer-power n) (car t1))
1836 (setq n (math-build-polynomial-expr poly high))
1837 (if (memq low '(0 1))
1839 (math-sub n (math-build-polynomial-expr poly
1840 (math-sub low 1))))))
1841 ((and (memq (car expr) '(+ -))
1842 (setq t1 (math-sum-rec (nth 1 expr) var low high)
1843 t2 (math-sum-rec (nth 2 expr) var low high))
1844 (not (and (math-expr-calls t1 '(calcFunc-sum))
1845 (math-expr-calls t2 '(calcFunc-sum)))))
1846 (list (car expr) t1 t2))
1847 ((and (eq (car expr) '*)
1848 (setq t1 (math-sum-const-factors expr var)))
1849 (math-mul (car t1) (math-sum-rec (cdr t1) var low high)))
1850 ((and (eq (car expr) '*) (memq (car-safe (nth 1 expr)) '(+ -)))
1851 (math-sum-rec (math-add-or-sub (math-mul (nth 1 (nth 1 expr))
1853 (math-mul (nth 2 (nth 1 expr))
1855 nil (eq (car (nth 1 expr)) '-))
1857 ((and (eq (car expr) '*) (memq (car-safe (nth 2 expr)) '(+ -)))
1858 (math-sum-rec (math-add-or-sub (math-mul (nth 1 expr)
1859 (nth 1 (nth 2 expr)))
1860 (math-mul (nth 1 expr)
1861 (nth 2 (nth 2 expr)))
1862 nil (eq (car (nth 2 expr)) '-))
1864 ((and (eq (car expr) '/)
1865 (not (math-primp (nth 1 expr)))
1866 (setq t1 (math-sum-const-factors (nth 1 expr) var)))
1868 (math-sum-rec (math-div (cdr t1) (nth 2 expr))
1870 ((and (eq (car expr) '/)
1871 (setq t1 (math-sum-const-factors (nth 2 expr) var)))
1872 (math-div (math-sum-rec (math-div (nth 1 expr) (cdr t1))
1875 ((eq (car expr) 'neg)
1876 (math-neg (math-sum-rec (nth 1 expr) var low high)))
1877 ((and (eq (car expr) '^)
1878 (not (math-expr-contains (nth 1 expr) var))
1879 (setq t1 (math-is-polynomial (nth 2 expr) var 1)))
1880 (let ((x (math-pow (nth 1 expr) (nth 1 t1))))
1881 (math-div (math-mul (math-sub (math-pow x (math-add 1 high))
1883 (math-pow (nth 1 expr) (car t1)))
1885 ((and (setq t1 (math-to-exponentials expr))
1886 (setq t1 (math-sum-rec t1 var low high))
1887 (not (math-expr-calls t1 '(calcFunc-sum))))
1889 ((memq (car expr) '(calcFunc-ln calcFunc-log10))
1890 (list (car expr) (calcFunc-prod (nth 1 expr) var low high)))
1891 ((and (eq (car expr) 'calcFunc-log)
1893 (not (math-expr-contains (nth 2 expr) var)))
1895 (calcFunc-prod (nth 1 expr) var low high)
1897 (if (equal val '(var nan var-nan)) (setq val nil))
1899 (let* ((math-tabulate-initial 0)
1900 (math-tabulate-function 'calcFunc-sum))
1901 (calcFunc-table expr var low high)))))
1903 (defun calcFunc-asum (expr var low &optional high step no-mul-flag)
1904 (or high (setq high low low 1))
1905 (if (and step (not (math-equal-int step 1)))
1906 (if (math-negp step)
1907 (math-mul (math-pow -1 low)
1908 (calcFunc-asum expr var high low (math-neg step) t))
1909 (let ((lo (math-simplify (math-div low step))))
1910 (if (math-num-integerp lo)
1911 (calcFunc-asum (math-normalize
1912 (math-expr-subst expr var
1913 (math-mul step var)))
1914 var lo (math-simplify (math-div high step)))
1915 (calcFunc-asum (math-normalize
1916 (math-expr-subst expr var
1917 (math-add (math-mul step var)
1920 (math-simplify (math-div (math-sub high low)
1922 (math-mul (if no-mul-flag 1 (math-pow -1 low))
1923 (calcFunc-sum (math-mul (math-pow -1 var) expr) var low high))))
1925 (defun math-sum-const-factors (expr var)
1929 (while (eq (car-safe p) '*)
1930 (if (math-expr-contains (nth 1 p) var)
1931 (setq not-const (cons (nth 1 p) not-const))
1932 (setq const (cons (nth 1 p) const)))
1934 (if (math-expr-contains p var)
1935 (setq not-const (cons p not-const))
1936 (setq const (cons p const)))
1938 (cons (let ((temp (car const)))
1939 (while (setq const (cdr const))
1940 (setq temp (list '* (car const) temp)))
1942 (let ((temp (or (car not-const) 1)))
1943 (while (setq not-const (cdr not-const))
1944 (setq temp (list '* (car not-const) temp)))
1947 (defvar math-sum-int-pow-cache (list '(0 1)))
1948 ;; Following is from CRC Math Tables, 27th ed, pp. 52-53.
1949 (defun math-sum-integer-power (pow)
1950 (let ((calc-prefer-frac t)
1951 (n (length math-sum-int-pow-cache)))
1953 (let* ((new (list 0 0))
1955 (pp (cdr (nth (1- n) math-sum-int-pow-cache)))
1960 (setq q (math-div (car pp) p)
1961 new (cons (math-mul q n) new)
1962 sum (math-add sum q)
1965 (setcar lin (math-sub 1 (math-mul n sum)))
1966 (setq math-sum-int-pow-cache
1967 (nconc math-sum-int-pow-cache (list (nreverse new)))
1969 (nth pow math-sum-int-pow-cache)))
1971 (defun math-to-exponentials (expr)
1974 (let ((x (nth 1 expr))
1975 (pi (if calc-symbolic-mode '(var pi var-pi) (math-pi)))
1976 (i (if calc-symbolic-mode '(var i var-i) '(cplx 0 1))))
1977 (cond ((eq (car expr) 'calcFunc-exp)
1978 (list '^ '(var e var-e) x))
1979 ((eq (car expr) 'calcFunc-sin)
1980 (or (eq calc-angle-mode 'rad)
1981 (setq x (list '/ (list '* x pi) 180)))
1983 (list '^ '(var e var-e) (list '* x i))
1984 (list '^ '(var e var-e)
1985 (list 'neg (list '* x i))))
1987 ((eq (car expr) 'calcFunc-cos)
1988 (or (eq calc-angle-mode 'rad)
1989 (setq x (list '/ (list '* x pi) 180)))
1991 (list '^ '(var e var-e)
1993 (list '^ '(var e var-e)
1994 (list 'neg (list '* x i))))
1996 ((eq (car expr) 'calcFunc-sinh)
1998 (list '^ '(var e var-e) x)
1999 (list '^ '(var e var-e) (list 'neg x)))
2001 ((eq (car expr) 'calcFunc-cosh)
2003 (list '^ '(var e var-e) x)
2004 (list '^ '(var e var-e) (list 'neg x)))
2008 (defun math-to-exps (expr)
2009 (cond (calc-symbolic-mode expr)
2011 (if (equal expr '(var e var-e)) (math-e) expr))
2012 ((and (eq (car expr) '^)
2013 (equal (nth 1 expr) '(var e var-e)))
2014 (list 'calcFunc-exp (nth 2 expr)))
2016 (cons (car expr) (mapcar 'math-to-exps (cdr expr))))))
2019 (defvar math-disable-prods nil)
2020 (defun calcFunc-prod (expr var &optional low high step)
2021 (if math-disable-prods (math-reject-arg))
2022 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2)))
2023 (math-prod-rec expr var low high step)))
2024 (math-disable-prods t))
2025 (math-normalize res)))
2027 (defun math-prod-rec (expr var &optional low high step)
2028 (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
2029 (and low (not high) (setq high '(var inf var-inf)))
2033 ((not (math-expr-contains expr var))
2034 (math-pow expr (math-add (math-div (math-sub high low) (or step 1))
2036 ((and step (not (math-equal-int step 1)))
2037 (if (math-negp step)
2038 (math-prod-rec expr var high low (math-neg step))
2039 (let ((lo (math-simplify (math-div low step))))
2040 (if (math-known-num-integerp lo)
2041 (math-prod-rec (math-normalize
2042 (math-expr-subst expr var
2043 (math-mul step var)))
2044 var lo (math-simplify (math-div high step)))
2045 (math-prod-rec (math-normalize
2046 (math-expr-subst expr var
2047 (math-add (math-mul step
2051 (math-simplify (math-div (math-sub high low)
2053 ((and (memq (car expr) '(* /))
2054 (setq t1 (math-prod-rec (nth 1 expr) var low high)
2055 t2 (math-prod-rec (nth 2 expr) var low high))
2056 (not (and (math-expr-calls t1 '(calcFunc-prod))
2057 (math-expr-calls t2 '(calcFunc-prod)))))
2058 (list (car expr) t1 t2))
2059 ((and (eq (car expr) '^)
2060 (not (math-expr-contains (nth 2 expr) var)))
2061 (math-pow (math-prod-rec (nth 1 expr) var low high)
2063 ((and (eq (car expr) '^)
2064 (not (math-expr-contains (nth 1 expr) var)))
2065 (math-pow (nth 1 expr)
2066 (calcFunc-sum (nth 2 expr) var low high)))
2067 ((eq (car expr) 'sqrt)
2068 (math-normalize (list 'calcFunc-sqrt
2069 (list 'calcFunc-prod (nth 1 expr)
2071 ((eq (car expr) 'neg)
2072 (math-mul (math-pow -1 (math-add (math-sub high low) 1))
2073 (math-prod-rec (nth 1 expr) var low high)))
2074 ((eq (car expr) 'calcFunc-exp)
2075 (list 'calcFunc-exp (calcFunc-sum (nth 1 expr) var low high)))
2076 ((and (setq t1 (math-is-polynomial expr var 1))
2079 ((or (and (math-equal-int (nth 1 t1) 1)
2080 (setq low (math-simplify
2081 (math-add low (car t1)))
2083 (math-add high (car t1)))))
2084 (and (math-equal-int (nth 1 t1) -1)
2087 (math-sub (car t1) high))
2089 (math-sub (car t1) t2)))))
2090 (if (or (math-zerop low) (math-zerop high))
2092 (if (and (or (math-negp low) (math-negp high))
2093 (or (math-num-integerp low)
2094 (math-num-integerp high)))
2095 (if (math-posp high)
2097 (math-mul (math-pow -1
2099 (math-add low high) 1))
2101 (list 'calcFunc-fact
2103 (list 'calcFunc-fact
2104 (math-sub -1 high)))))
2106 (list 'calcFunc-fact high)
2107 (list 'calcFunc-fact (math-sub low 1))))))
2108 ((and (or (and (math-equal-int (nth 1 t1) 2)
2109 (setq t2 (math-simplify
2110 (math-add (math-mul low 2)
2113 (math-add (math-mul high 2)
2115 (and (math-equal-int (nth 1 t1) -2)
2116 (setq t2 (math-simplify
2123 (or (math-integerp t2)
2124 (and (math-messy-integerp t2)
2125 (setq t2 (math-trunc t2)))
2127 (and (math-messy-integerp t3)
2128 (setq t3 (math-trunc t3)))))
2129 (if (or (math-zerop t2) (math-zerop t3))
2131 (if (or (math-evenp t2) (math-evenp t3))
2132 (if (or (math-negp t2) (math-negp t3))
2133 (if (math-posp high)
2136 (list 'calcFunc-dfact
2138 (list 'calcFunc-dfact
2141 (list 'calcFunc-dfact t3)
2142 (list 'calcFunc-dfact
2147 (list '/ (list '- (list '- t2 t3)
2151 (list 'calcFunc-dfact
2153 (list 'calcFunc-dfact
2157 (list 'calcFunc-dfact t3)
2158 (list 'calcFunc-dfact
2162 (if (equal val '(var nan var-nan)) (setq val nil))
2164 (let* ((math-tabulate-initial 1)
2165 (math-tabulate-function 'calcFunc-prod))
2166 (calcFunc-table expr var low high)))))
2171 (defvar math-solve-ranges nil)
2172 ;;; Attempt to reduce lhs = rhs to solve-var = rhs', where solve-var appears
2173 ;;; in lhs but not in rhs or rhs'; return rhs'.
2174 ;;; Uses global values: solve-*.
2175 (defun math-try-solve-for (lhs rhs &optional sign no-poly)
2177 (cond ((equal lhs solve-var)
2178 (setq math-solve-sign sign)
2179 (if (eq solve-full 'all)
2180 (let ((vec (list 'vec (math-evaluate-expr rhs)))
2182 (while math-solve-ranges
2183 (setq p (car math-solve-ranges)
2186 (while (setq p (cdr p))
2187 (setq newvec (nconc newvec
2188 (cdr (math-expr-subst
2189 vec var (car p))))))
2191 math-solve-ranges (cdr math-solve-ranges)))
2192 (math-normalize vec))
2196 ((and (eq (car lhs) '-)
2197 (eq (car-safe (nth 1 lhs)) (car-safe (nth 2 lhs)))
2199 (= (length (nth 1 lhs)) 2)
2200 (= (length (nth 2 lhs)) 2)
2201 (setq t1 (get (car (nth 1 lhs)) 'math-inverse))
2202 (setq t2 (funcall t1 '(var SOLVEDUM SOLVEDUM)))
2203 (eq (math-expr-contains-count t2 '(var SOLVEDUM SOLVEDUM)) 1)
2204 (setq t3 (math-solve-above-dummy t2))
2205 (setq t1 (math-try-solve-for (math-sub (nth 1 (nth 1 lhs))
2208 (nth 1 (nth 2 lhs))))
2211 ((eq (car lhs) 'neg)
2212 (math-try-solve-for (nth 1 lhs) (math-neg rhs)
2213 (and sign (- sign))))
2214 ((and (not (eq solve-full 't)) (math-try-solve-prod)))
2216 (setq t2 (math-decompose-poly lhs solve-var 15 rhs)))
2217 (setq t1 (cdr (nth 1 t2))
2218 t1 (let ((math-solve-ranges math-solve-ranges))
2219 (cond ((= (length t1) 5)
2220 (apply 'math-solve-quartic (car t2) t1))
2222 (apply 'math-solve-cubic (car t2) t1))
2224 (apply 'math-solve-quadratic (car t2) t1))
2226 (apply 'math-solve-linear (car t2) sign t1))
2228 (math-poly-all-roots (car t2) t1))
2229 (calc-symbolic-mode nil)
2233 (math-poly-any-root (reverse t1) 0 t)
2236 (if (eq (nth 2 t2) 1)
2238 (math-solve-prod t1 (math-try-solve-for (nth 2 t2) 0 nil t)))
2239 (calc-record-why "*Unable to find a symbolic solution")
2241 ((and (math-solve-find-root-term lhs nil)
2242 (eq (math-expr-contains-count lhs t1) 1)) ; just in case
2243 (math-try-solve-for (math-simplify
2244 (math-sub (if (or t3 (math-evenp t2))
2246 (math-neg (math-pow t1 t2)))
2248 (math-sub (math-normalize
2255 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
2256 (math-try-solve-for (nth 2 lhs)
2257 (math-sub rhs (nth 1 lhs))
2259 ((not (math-expr-contains (nth 2 lhs) solve-var))
2260 (math-try-solve-for (nth 1 lhs)
2261 (math-sub rhs (nth 2 lhs))
2263 ((eq (car lhs) 'calcFunc-eq)
2264 (math-try-solve-for (math-sub (nth 1 lhs) (nth 2 lhs))
2267 (cond ((or (and (eq (car-safe (nth 1 lhs)) 'calcFunc-sin)
2268 (eq (car-safe (nth 2 lhs)) 'calcFunc-cos))
2269 (and (eq (car-safe (nth 1 lhs)) 'calcFunc-cos)
2270 (eq (car-safe (nth 2 lhs)) 'calcFunc-sin)))
2271 (math-try-solve-for (math-sub (nth 1 lhs)
2272 (list (car (nth 1 lhs))
2274 (math-quarter-circle t)
2275 (nth 1 (nth 2 lhs)))))
2277 ((not (math-expr-contains (nth 1 lhs) solve-var))
2278 (math-try-solve-for (nth 2 lhs)
2279 (math-sub (nth 1 lhs) rhs)
2280 (and sign (- sign))))
2281 ((not (math-expr-contains (nth 2 lhs) solve-var))
2282 (math-try-solve-for (nth 1 lhs)
2283 (math-add rhs (nth 2 lhs))
2285 ((and (eq solve-full 't) (math-try-solve-prod)))
2286 ((and (eq (car lhs) '%)
2287 (not (math-expr-contains (nth 2 lhs) solve-var)))
2288 (math-try-solve-for (nth 1 lhs) (math-add rhs
2291 ((eq (car lhs) 'calcFunc-log)
2292 (cond ((not (math-expr-contains (nth 2 lhs) solve-var))
2293 (math-try-solve-for (nth 1 lhs) (math-pow (nth 2 lhs) rhs)))
2294 ((not (math-expr-contains (nth 1 lhs) solve-var))
2295 (math-try-solve-for (nth 2 lhs) (math-pow
2297 (math-div 1 rhs))))))
2298 ((and (= (length lhs) 2)
2300 (setq t1 (get (car lhs) 'math-inverse))
2301 (setq t2 (funcall t1 rhs)))
2302 (setq t1 (get (car lhs) 'math-inverse-sign))
2303 (math-try-solve-for (nth 1 lhs) (math-normalize t2)
2307 (funcall t1 lhs sign)))))
2308 ((and (symbolp (car lhs))
2309 (setq t1 (get (car lhs) 'math-inverse-n))
2310 (setq t2 (funcall t1 lhs rhs)))
2312 ((setq t1 (math-expand-formula lhs))
2313 (math-try-solve-for t1 rhs sign))
2315 (calc-record-why "*No inverse known" lhs)
2319 (defun math-try-solve-prod ()
2320 (cond ((eq (car lhs) '*)
2321 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
2322 (math-try-solve-for (nth 2 lhs)
2323 (math-div rhs (nth 1 lhs))
2324 (math-solve-sign sign (nth 1 lhs))))
2325 ((not (math-expr-contains (nth 2 lhs) solve-var))
2326 (math-try-solve-for (nth 1 lhs)
2327 (math-div rhs (nth 2 lhs))
2328 (math-solve-sign sign (nth 2 lhs))))
2330 (math-solve-prod (let ((math-solve-ranges math-solve-ranges))
2331 (math-try-solve-for (nth 2 lhs) 0))
2332 (math-try-solve-for (nth 1 lhs) 0)))))
2334 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
2335 (math-try-solve-for (nth 2 lhs)
2336 (math-div (nth 1 lhs) rhs)
2337 (math-solve-sign sign (nth 1 lhs))))
2338 ((not (math-expr-contains (nth 2 lhs) solve-var))
2339 (math-try-solve-for (nth 1 lhs)
2340 (math-mul rhs (nth 2 lhs))
2341 (math-solve-sign sign (nth 2 lhs))))
2342 ((setq t1 (math-try-solve-for (math-sub (nth 1 lhs)
2343 (math-mul (nth 2 lhs)
2348 (cond ((not (math-expr-contains (nth 1 lhs) solve-var))
2351 (math-add (math-normalize
2352 (list 'calcFunc-log rhs (nth 1 lhs)))
2355 (math-mul '(var pi var-pi)
2359 (list 'calcFunc-ln (nth 1 lhs)))))))
2360 ((not (math-expr-contains (nth 2 lhs) solve-var))
2361 (cond ((and (integerp (nth 2 lhs))
2363 (setq t1 (math-integer-log2 (nth 2 lhs))))
2365 (if (and (eq solve-full t)
2366 (math-known-realp (nth 1 lhs)))
2368 (while (>= (setq t1 (1- t1)) 0)
2369 (setq t2 (list 'calcFunc-sqrt t2)))
2370 (setq t2 (math-solve-get-sign t2)))
2371 (while (>= (setq t1 (1- t1)) 0)
2372 (setq t2 (math-solve-get-sign
2374 (list 'calcFunc-sqrt t2))))))
2377 (math-normalize t2)))
2378 ((math-looks-negp (nth 2 lhs))
2380 (list '^ (nth 1 lhs) (math-neg (nth 2 lhs)))
2382 ((and (eq solve-full t)
2383 (Math-integerp (nth 2 lhs))
2384 (math-known-realp (nth 1 lhs)))
2385 (setq t1 (math-normalize
2386 (list 'calcFunc-nroot rhs (nth 2 lhs))))
2387 (if (math-evenp (nth 2 lhs))
2388 (setq t1 (math-solve-get-sign t1)))
2392 (math-oddp (nth 2 lhs))
2393 (math-solve-sign sign (nth 2 lhs)))))
2394 (t (math-try-solve-for
2399 (if (Math-realp (nth 2 lhs))
2404 (and (integerp (nth 2 lhs))
2407 (math-div (nth 2 lhs) 2))
2414 (and (integerp (nth 2 lhs))
2419 (list 'calcFunc-nroot
2423 (math-oddp (nth 2 lhs))
2424 (math-solve-sign sign (nth 2 lhs)))))))))
2427 (defun math-solve-prod (lsoln rsoln)
2432 ((eq solve-full 'all)
2433 (cons 'vec (append (cdr lsoln) (cdr rsoln))))
2436 (list 'calcFunc-gt (math-solve-get-sign 1) 0)
2441 ;;; This deals with negative, fractional, and symbolic powers of "x".
2442 (defun math-solve-poly-funny-powers (sub-rhs) ; uses "t1", "t2"
2444 (let ((pp math-poly-neg-powers)
2447 (setq fac (math-pow (car pp) (or math-poly-mult-powers 1))
2448 t1 (math-mul t1 fac)
2449 rhs (math-mul rhs fac)
2451 (if sub-rhs (setq t1 (math-sub t1 rhs)))
2452 (let ((math-poly-neg-powers nil))
2453 (setq t2 (math-mul (or math-poly-mult-powers 1)
2454 (let ((calc-prefer-frac t))
2455 (math-div 1 math-poly-frac-powers)))
2456 t1 (math-is-polynomial (math-simplify (calcFunc-expand t1)) b 50))))
2458 ;;; This converts "a x^8 + b x^5 + c x^2" to "(a (x^3)^2 + b (x^3) + c) * x^2".
2459 (defun math-solve-crunch-poly (max-degree) ; uses "t1", "t3"
2461 (while (and t1 (Math-zerop (car t1)))
2465 (let* ((degree (1- (length t1)))
2467 (while (and (> scale 1) (= (car t3) 1))
2468 (and (= (% degree scale) 0)
2474 (if (= (% n scale) 0)
2475 (setq new-t1 (nconc new-t1 (list (car p))))
2476 (or (Math-zerop (car p))
2481 (setq t3 (cons scale (cdr t3))
2483 (setq scale (1- scale)))
2484 (setq t3 (list (math-mul (car t3) t2) (math-mul count t2)))
2485 (<= (1- (length t1)) max-degree)))))
2487 (defun calcFunc-poly (expr var &optional degree)
2489 (or (natnump degree) (math-reject-arg degree 'fixnatnump))
2491 (let ((p (math-is-polynomial expr var degree 'gen)))
2496 (math-reject-arg expr "Expected a polynomial"))))
2498 (defun calcFunc-gpoly (expr var &optional degree)
2500 (or (natnump degree) (math-reject-arg degree 'fixnatnump))
2502 (let* ((math-poly-base-variable var)
2503 (d (math-decompose-poly expr var degree nil)))
2506 (math-reject-arg expr "Expected a polynomial"))))
2508 (defun math-decompose-poly (lhs solve-var degree sub-rhs)
2509 (let ((rhs (or sub-rhs 1))
2511 (setq t2 (math-polynomial-base
2515 (let ((math-poly-neg-powers '(1))
2516 (math-poly-mult-powers nil)
2517 (math-poly-frac-powers 1)
2518 (math-poly-exp-base t))
2519 (and (not (equal b lhs))
2520 (or (not (memq (car-safe b) '(+ -))) sub-rhs)
2521 (setq t3 '(1 0) t2 1
2522 t1 (math-is-polynomial lhs b 50))
2523 (if (and (equal math-poly-neg-powers '(1))
2524 (memq math-poly-mult-powers '(nil 1))
2525 (eq math-poly-frac-powers 1)
2527 (setq t1 (cons (math-sub (car t1) rhs)
2529 (math-solve-poly-funny-powers sub-rhs))
2530 (math-solve-crunch-poly degree)
2531 (or (math-expr-contains b solve-var)
2532 (math-expr-contains (car t3) solve-var))))))))
2534 (list (math-pow t2 (car t3))
2537 (math-pow t2 (nth 1 t3))
2538 (math-div (math-pow t2 (nth 1 t3)) rhs))))))
2540 (defun math-solve-linear (var sign b a)
2541 (math-try-solve-for var
2542 (math-div (math-neg b) a)
2543 (math-solve-sign sign a)
2546 (defun math-solve-quadratic (var c b a)
2549 (if (math-looks-evenp b)
2550 (let ((halfb (math-div b 2)))
2554 (math-solve-get-sign
2556 (list 'calcFunc-sqrt
2557 (math-add (math-sqr halfb)
2558 (math-mul (math-neg c) a))))))
2563 (math-solve-get-sign
2565 (list 'calcFunc-sqrt
2566 (math-add (math-sqr b)
2567 (math-mul 4 (math-mul (math-neg c) a)))))))
2571 (defun math-solve-cubic (var d c b a)
2572 (let* ((p (math-div b a))
2576 (aa (math-sub q (math-div psqr 3)))
2578 (math-div (math-sub (math-mul 2 (math-mul psqr p))
2579 (math-mul 9 (math-mul p q)))
2583 (math-try-solve-for (math-pow (math-add var (math-div p 3)) 3)
2584 (math-neg bb) nil t)
2587 (math-mul (math-add var (math-div p 3))
2588 (math-add (math-sqr (math-add var (math-div p 3)))
2591 (setq m (math-mul 2 (list 'calcFunc-sqrt (math-div aa -3))))
2600 (math-sub (list 'calcFunc-arccos
2601 (math-div (math-mul 3 bb)
2605 (math-add 1 (math-solve-get-int
2608 calc-symbolic-mode))))
2613 (defun math-solve-quartic (var d c b a aa)
2614 (setq a (math-div a aa))
2615 (setq b (math-div b aa))
2616 (setq c (math-div c aa))
2617 (setq d (math-div d aa))
2620 (let* ((asqr (math-sqr a))
2621 (asqr4 (math-div asqr 4))
2622 (y (let ((solve-full nil)
2624 (math-solve-cubic solve-var
2626 (math-mul 4 (math-mul b d))
2629 (math-sub (math-mul a c)
2633 (rsqr (math-add (math-sub asqr4 b) y))
2634 (r (list 'calcFunc-sqrt rsqr))
2635 (sign1 (math-solve-get-sign 1))
2636 (de (list 'calcFunc-sqrt
2638 (math-sub (math-mul 3 asqr4)
2640 (if (Math-zerop rsqr)
2644 (list 'calcFunc-sqrt
2645 (math-sub (math-sqr y)
2651 (math-mul 4 (math-mul a b))
2657 (math-sub (math-add (math-mul sign1 (math-div r 2))
2658 (math-solve-get-sign (math-div de 2)))
2662 (defvar math-symbolic-solve nil)
2663 (defvar math-int-coefs nil)
2664 (defun math-poly-all-roots (var p &optional math-factoring)
2666 (let* ((math-symbolic-solve calc-symbolic-mode)
2668 (deg (1- (length p)))
2669 (orig-p (reverse p))
2670 (math-int-coefs nil)
2671 (math-int-scale nil)
2672 (math-double-roots nil)
2673 (math-int-factors nil)
2674 (math-int-threshold nil)
2676 ;; If rational coefficients, look for exact rational factors.
2677 (while (and pp (Math-ratp (car pp)))
2680 (if (or math-factoring math-symbolic-solve)
2682 (let ((lead (car orig-p))
2683 (calc-prefer-frac t)
2684 (scale (apply 'math-lcm-denoms p)))
2685 (setq math-int-scale (math-abs (math-mul scale lead))
2686 math-int-threshold (math-div '(float 5 -2) math-int-scale)
2687 math-int-coefs (cdr (math-div (cons 'vec orig-p) lead)))))
2689 (let ((calc-prefer-frac nil)
2690 (calc-symbolic-mode nil)
2692 (def-p (copy-sequence orig-p)))
2694 (if (Math-numberp (car pp))
2697 (while (> deg (if math-symbolic-solve 2 4))
2698 (let* ((x (math-poly-any-root def-p '(float 0 0) nil))
2700 (if (and (eq (car-safe x) 'cplx)
2701 (math-nearly-zerop (nth 2 x) (nth 1 x)))
2702 (setq x (calcFunc-re x)))
2704 (setq roots (cons x roots)))
2705 (or (math-numberp x)
2706 (setq x (math-evaluate-expr x)))
2709 (while (setq pp (cdr pp))
2712 (setq b (math-add (math-mul x b) c)))
2713 (setq def-p (cdr def-p)
2715 (setq p (reverse def-p))))
2717 (let ((solve-var '(var DUMMY var-DUMMY))
2718 (math-solve-sign nil)
2719 (math-solve-ranges nil)
2721 (if (= (length p) (length math-int-coefs))
2722 (setq p (reverse math-int-coefs)))
2723 (setq roots (append (cdr (apply (cond ((= deg 2)
2724 'math-solve-quadratic)
2728 'math-solve-quartic))
2732 (setq roots (cons (math-div (math-neg (car p)) (nth 1 p))
2737 (math-poly-integer-root (car roots))
2738 (setq roots (cdr roots)))
2739 (list math-int-factors (nreverse math-int-coefs) math-int-scale))
2740 (let ((vec nil) res)
2742 (let ((root (car roots))
2743 (solve-full (and solve-full 'all)))
2744 (if (math-floatp root)
2745 (setq root (math-poly-any-root orig-p root t)))
2746 (setq vec (append vec
2747 (cdr (or (math-try-solve-for var root nil t)
2748 (throw 'ouch nil))))))
2749 (setq roots (cdr roots)))
2750 (setq vec (cons 'vec (nreverse vec)))
2751 (if math-symbolic-solve
2752 (setq vec (math-normalize vec)))
2753 (if (eq solve-full t)
2754 (list 'calcFunc-subscr
2756 (math-solve-get-int 1 (1- (length orig-p)) 1))
2759 (defun math-lcm-denoms (&rest fracs)
2762 (if (eq (car-safe (car fracs)) 'frac)
2763 (setq den (calcFunc-lcm den (nth 2 (car fracs)))))
2764 (setq fracs (cdr fracs)))
2767 (defun math-poly-any-root (p x polish) ; p is a reverse poly coeff list
2768 (let* ((newt (if (math-zerop x)
2769 (math-poly-newton-root
2770 p '(cplx (float 123 -6) (float 1 -4)) 4)
2771 (math-poly-newton-root p x 4)))
2772 (res (if (math-zerop (cdr newt))
2774 (if (and (math-lessp (cdr newt) '(float 1 -3)) (not polish))
2775 (setq newt (math-poly-newton-root p (car newt) 30)))
2776 (if (math-zerop (cdr newt))
2778 (math-poly-laguerre-root p x polish)))))
2779 (and math-symbolic-solve (math-floatp res)
2783 (defun math-poly-newton-root (p x iters)
2784 (let* ((calc-prefer-frac nil)
2785 (calc-symbolic-mode nil)
2786 (try-integer math-int-coefs)
2788 (while (and (> (setq iters (1- iters)) 0)
2790 (math-working "newton" x)
2793 (while (setq pp (cdr pp))
2794 (setq d (math-add (math-mul x d) b)
2795 b (math-add (math-mul x b) (car pp))))
2796 (not (math-zerop d)))
2798 (setq dx (math-div b d)
2801 (let ((adx (math-abs-approx dx)))
2802 (and (math-lessp adx math-int-threshold)
2803 (let ((iroot (math-poly-integer-root x)))
2806 (setq try-integer nil))))))
2807 (or (not (or (eq dx 0)
2808 (math-nearly-zerop dx (math-abs-approx x))))
2809 (progn (setq dx 0) nil)))))
2810 (cons x (if (math-zerop x)
2811 1 (math-div (math-abs-approx dx) (math-abs-approx x))))))
2813 (defun math-poly-integer-root (x)
2814 (and (math-lessp (calcFunc-xpon (math-abs-approx x)) calc-internal-prec)
2816 (let* ((calc-prefer-frac t)
2817 (xre (calcFunc-re x))
2818 (xim (calcFunc-im x))
2819 (xresq (math-sqr xre))
2820 (ximsq (math-sqr xim)))
2821 (if (math-lessp ximsq (calcFunc-scf xresq -1))
2822 ;; Look for linear factor
2823 (let* ((rnd (math-div (math-round (math-mul xre math-int-scale))
2825 (icp math-int-coefs)
2828 (while (setq icp (cdr icp))
2829 (setq newcoef (cons rem newcoef)
2830 rem (math-add (car icp)
2831 (math-mul rem rnd))))
2832 (and (math-zerop rem)
2834 (setq math-int-coefs (nreverse newcoef)
2835 math-int-factors (cons (list (math-neg rnd))
2838 ;; Look for irreducible quadratic factor
2839 (let* ((rnd1 (math-div (math-round
2840 (math-mul xre (math-mul -2 math-int-scale)))
2842 (sqscale (math-sqr math-int-scale))
2843 (rnd0 (math-div (math-round (math-mul (math-add xresq ximsq)
2846 (rem1 (car math-int-coefs))
2847 (icp (cdr math-int-coefs))
2850 (found (assoc (list rnd0 rnd1 (math-posp xim))
2854 (setq math-double-roots (delq found math-double-roots)
2856 (while (setq icp (cdr icp))
2858 newcoef (cons rem1 newcoef)
2859 rem1 (math-sub rem0 (math-mul this rnd1))
2860 rem0 (math-sub (car icp) (math-mul this rnd0)))))
2861 (and (math-zerop rem0)
2863 (let ((aa (math-div rnd1 -2)))
2864 (or found (setq math-int-coefs (reverse newcoef)
2865 math-double-roots (cons (list
2870 math-int-factors (cons (cons rnd0 rnd1)
2873 (let ((calc-symbolic-mode math-symbolic-solve))
2874 (math-mul (math-sqrt (math-sub (math-sqr aa)
2876 (if (math-negp xim) -1 1)))))))))))
2878 ;;; The following routine is from Numerical Recipes, section 9.5.
2879 (defun math-poly-laguerre-root (p x polish)
2880 (let* ((calc-prefer-frac nil)
2881 (calc-symbolic-mode nil)
2884 (try-newt (not polish))
2888 (and (or (< (setq iters (1+ iters)) 50)
2889 (math-reject-arg x "*Laguerre's method failed to converge"))
2890 (let ((err (math-abs-approx (car p)))
2891 (abx (math-abs-approx x))
2895 (while (setq pp (cdr pp))
2896 (setq f (math-add (math-mul x f) d)
2897 d (math-add (math-mul x d) b)
2898 b (math-add (math-mul x b) (car pp))
2899 err (math-add (math-abs-approx b) (math-mul abx err))))
2900 (math-lessp (calcFunc-scf err (- -2 calc-internal-prec))
2901 (math-abs-approx b)))
2902 (or (not (math-zerop d))
2903 (not (math-zerop f))
2905 (setq x (math-pow (math-neg b) (list 'frac 1 m)))
2907 (let* ((g (math-div d b))
2909 (h (math-sub g2 (math-mul 2 (math-div f b))))
2911 (math-mul (1- m) (math-sub (math-mul m h) g2))))
2912 (gp (math-add g sq))
2913 (gm (math-sub g sq)))
2914 (if (math-lessp (calcFunc-abssqr gp) (calcFunc-abssqr gm))
2916 (setq dx (math-div m gp)
2919 (math-lessp (math-abs-approx dx)
2920 (calcFunc-scf (math-abs-approx x) -3)))
2921 (let ((newt (math-poly-newton-root p x1 7)))
2924 (if (math-zerop (cdr newt))
2925 (setq x (car newt) x1 x)
2926 (if (math-lessp (cdr newt) '(float 1 -6))
2927 (let ((newt2 (math-poly-newton-root
2929 (if (math-zerop (cdr newt2))
2930 (setq x (car newt2) x1 x)
2931 (setq x (car newt))))))))
2933 (math-nearly-equal x x1))))
2934 (let ((cdx (math-abs-approx dx)))
2939 (math-lessp cdx dxold)
2942 (let ((digs (calcFunc-xpon
2943 (math-div (math-abs-approx x) cdx))))
2945 "*Could not attain full precision")
2947 (let ((calc-internal-prec (max 3 digs)))
2948 (setq x (math-normalize x))))))
2952 (math-lessp (calcFunc-scf (math-abs-approx x)
2953 (- calc-internal-prec))
2955 (or (and (math-floatp x)
2956 (math-poly-integer-root x))
2959 (defun math-solve-above-dummy (x)
2960 (and (not (Math-primp x))
2961 (if (and (equal (nth 1 x) '(var SOLVEDUM SOLVEDUM))
2965 (while (and (setq x (cdr x))
2966 (not (setq res (math-solve-above-dummy (car x))))))
2969 (defun math-solve-find-root-term (x neg) ; sets "t2", "t3"
2970 (if (math-solve-find-root-in-prod x)
2973 (and (memq (car-safe x) '(+ -))
2974 (or (math-solve-find-root-term (nth 1 x) neg)
2975 (math-solve-find-root-term (nth 2 x)
2976 (if (eq (car x) '-) (not neg) neg))))))
2978 (defun math-solve-find-root-in-prod (x)
2980 (math-expr-contains x solve-var)
2981 (or (and (eq (car x) 'calcFunc-sqrt)
2983 (and (eq (car x) '^)
2984 (or (and (memq (math-quarter-integer (nth 2 x)) '(1 2 3))
2986 (and (eq (car-safe (nth 2 x)) 'frac)
2987 (eq (nth 2 (nth 2 x)) 3)
2989 (and (memq (car x) '(* /))
2990 (or (and (not (math-expr-contains (nth 1 x) solve-var))
2991 (math-solve-find-root-in-prod (nth 2 x)))
2992 (and (not (math-expr-contains (nth 2 x) solve-var))
2993 (math-solve-find-root-in-prod (nth 1 x))))))))
2996 (defun math-solve-system (exprs solve-vars solve-full)
2997 (setq exprs (mapcar 'list (if (Math-vectorp exprs)
3000 solve-vars (if (Math-vectorp solve-vars)
3003 (or (let ((math-solve-simplifying nil))
3004 (math-solve-system-rec exprs solve-vars nil))
3005 (let ((math-solve-simplifying t))
3006 (math-solve-system-rec exprs solve-vars nil))))
3008 ;;; The following backtracking solver works by choosing a variable
3009 ;;; and equation, and trying to solve the equation for the variable.
3010 ;;; If it succeeds it calls itself recursively with that variable and
3011 ;;; equation removed from their respective lists, and with the solution
3012 ;;; added to solns as well as being substituted into all existing
3013 ;;; equations. The algorithm terminates when any solution path
3014 ;;; manages to remove all the variables from var-list.
3016 ;;; To support calcFunc-roots, entries in eqn-list and solns are
3017 ;;; actually lists of equations.
3019 (defun math-solve-system-rec (eqn-list var-list solns)
3024 ;; Try each variable in turn.
3030 (elim (eq (car-safe vv) 'calcFunc-elim)))
3032 (setq vv (nth 1 vv)))
3034 ;; Try each equation in turn.
3043 ;; Try to solve for vv the list of equations e2.
3045 (setq res2 (or (and (eq (car e2) eprev)
3047 (math-solve-for (car e2) 0 vv
3049 (setq eprev (car e2)
3050 res (cons (if (eq solve-full 'all)
3058 ;; Found a solution. Now try other variables.
3059 (setq res (nreverse res)
3060 res (math-solve-system-rec
3062 'math-solve-system-subst
3064 (copy-sequence eqn-list)))
3065 (delq (car v) (copy-sequence var-list))
3066 (let ((math-solve-simplifying nil)
3072 (math-solve-system-subst
3077 (cons (cons vv (apply 'append res))
3085 ;; Eliminated all variables, so now put solution into the proper format.
3086 (setq solns (sort solns
3089 (not (memq (car x) (memq (car y) solve-vars)))))))
3090 (if (eq solve-full 'all)
3095 (mapcar (function (lambda (x) (cons 'vec (cdr x)))) solns)
3096 (mapcar (function (lambda (x) (cons 'vec x))) eqn-list)))))
3100 (mapcar (function (lambda (x) (cons 'calcFunc-eq x))) solns)
3101 (mapcar 'car eqn-list)))))))
3103 (defun math-solve-system-subst (x) ; uses "res" and "v"
3107 (setq accum (nconc accum
3110 (if math-solve-simplifying
3112 (math-expr-subst (car x) vv r))
3113 (math-expr-subst (car x) vv r))))
3120 (defun math-get-from-counter (name)
3121 (let ((ctr (assq name calc-command-flags)))
3123 (setcdr ctr (1+ (cdr ctr)))
3124 (setq ctr (cons name 1)
3125 calc-command-flags (cons ctr calc-command-flags)))
3128 (defun math-solve-get-sign (val)
3129 (setq val (math-simplify val))
3130 (if (and (eq (car-safe val) '*)
3131 (Math-numberp (nth 1 val)))
3132 (list '* (nth 1 val) (math-solve-get-sign (nth 2 val)))
3133 (and (eq (car-safe val) 'calcFunc-sqrt)
3134 (eq (car-safe (nth 1 val)) '^)
3135 (setq val (math-normalize (list '^
3137 (math-div (nth 2 (nth 1 val)) 2)))))
3139 (if (and (calc-var-value 'var-GenCount)
3140 (Math-natnump var-GenCount)
3141 (not (eq solve-full 'all)))
3143 (math-mul (list 'calcFunc-as var-GenCount) val)
3144 (setq var-GenCount (math-add var-GenCount 1))
3145 (calc-refresh-evaltos 'var-GenCount))
3146 (let* ((var (concat "s" (int-to-string (math-get-from-counter 'solve-sign))))
3147 (var2 (list 'var (intern var) (intern (concat "var-" var)))))
3148 (if (eq solve-full 'all)
3149 (setq math-solve-ranges (cons (list var2 1 -1)
3150 math-solve-ranges)))
3151 (math-mul var2 val)))
3152 (calc-record-why "*Choosing positive solution")
3155 (defun math-solve-get-int (val &optional range first)
3157 (if (and (calc-var-value 'var-GenCount)
3158 (Math-natnump var-GenCount)
3159 (not (eq solve-full 'all)))
3161 (math-mul val (list 'calcFunc-an var-GenCount))
3162 (setq var-GenCount (math-add var-GenCount 1))
3163 (calc-refresh-evaltos 'var-GenCount))
3164 (let* ((var (concat "n" (int-to-string
3165 (math-get-from-counter 'solve-int))))
3166 (var2 (list 'var (intern var) (intern (concat "var-" var)))))
3167 (if (and range (eq solve-full 'all))
3168 (setq math-solve-ranges (cons (cons var2
3169 (cdr (calcFunc-index
3170 range (or first 0))))
3171 math-solve-ranges)))
3172 (math-mul val var2)))
3173 (calc-record-why "*Choosing 0 for arbitrary integer in solution")
3176 (defun math-solve-sign (sign expr)
3178 (let ((s1 (math-possible-signs expr)))
3179 (cond ((memq s1 '(4 6))
3184 (defun math-looks-evenp (expr)
3185 (if (Math-integerp expr)
3187 (if (memq (car expr) '(* /))
3188 (math-looks-evenp (nth 1 expr)))))
3190 (defun math-solve-for (lhs rhs solve-var solve-full &optional sign)
3191 (if (math-expr-contains rhs solve-var)
3192 (math-solve-for (math-sub lhs rhs) 0 solve-var solve-full)
3193 (and (math-expr-contains lhs solve-var)
3194 (math-with-extra-prec 1
3195 (let* ((math-poly-base-variable solve-var)
3196 (res (math-try-solve-for lhs rhs sign)))
3197 (if (and (eq solve-full 'all)
3198 (math-known-realp solve-var))
3199 (let ((old-len (length res))
3204 (and (not (memq (car-safe x)
3208 new-len (length res))
3209 (if (< new-len old-len)
3210 (calc-record-why (if (= new-len 1)
3211 "*All solutions were complex"
3213 "*Omitted %d complex solutions"
3214 (- old-len new-len)))))))
3217 (defun math-solve-eqn (expr var full)
3218 (if (memq (car-safe expr) '(calcFunc-neq calcFunc-lt calcFunc-gt
3219 calcFunc-leq calcFunc-geq))
3220 (let ((res (math-solve-for (cons '- (cdr expr))
3222 (if (eq (car expr) 'calcFunc-neq) nil 1))))
3224 (if (eq math-solve-sign 1)
3225 (list (car expr) var res)
3226 (if (eq math-solve-sign -1)
3227 (list (car expr) res var)
3228 (or (eq (car expr) 'calcFunc-neq)
3230 "*Can't determine direction of inequality"))
3231 (and (memq (car expr) '(calcFunc-neq calcFunc-lt calcFunc-gt))
3232 (list 'calcFunc-neq var res))))))
3233 (let ((res (math-solve-for expr 0 var full)))
3235 (list 'calcFunc-eq var res)))))
3237 (defun math-reject-solution (expr var func)
3238 (if (math-expr-contains expr var)
3239 (or (equal (car calc-next-why) '(* "Unable to find a symbolic solution"))
3240 (calc-record-why "*Unable to find a solution")))
3241 (list func expr var))
3243 (defun calcFunc-solve (expr var)
3244 (or (if (or (Math-vectorp expr) (Math-vectorp var))
3245 (math-solve-system expr var nil)
3246 (math-solve-eqn expr var nil))
3247 (math-reject-solution expr var 'calcFunc-solve)))
3249 (defun calcFunc-fsolve (expr var)
3250 (or (if (or (Math-vectorp expr) (Math-vectorp var))
3251 (math-solve-system expr var t)
3252 (math-solve-eqn expr var t))
3253 (math-reject-solution expr var 'calcFunc-fsolve)))
3255 (defun calcFunc-roots (expr var)
3256 (let ((math-solve-ranges nil))
3257 (or (if (or (Math-vectorp expr) (Math-vectorp var))
3258 (math-solve-system expr var 'all)
3259 (math-solve-for expr 0 var 'all))
3260 (math-reject-solution expr var 'calcFunc-roots))))
3262 (defun calcFunc-finv (expr var)
3263 (let ((res (math-solve-for expr math-integ-var var nil)))
3265 (math-normalize (math-expr-subst res math-integ-var var))
3266 (math-reject-solution expr var 'calcFunc-finv))))
3268 (defun calcFunc-ffinv (expr var)
3269 (let ((res (math-solve-for expr math-integ-var var t)))
3271 (math-normalize (math-expr-subst res math-integ-var var))
3272 (math-reject-solution expr var 'calcFunc-finv))))
3275 (put 'calcFunc-inv 'math-inverse
3276 (function (lambda (x) (math-div 1 x))))
3277 (put 'calcFunc-inv 'math-inverse-sign -1)
3279 (put 'calcFunc-sqrt 'math-inverse
3280 (function (lambda (x) (math-sqr x))))
3282 (put 'calcFunc-conj 'math-inverse
3283 (function (lambda (x) (list 'calcFunc-conj x))))
3285 (put 'calcFunc-abs 'math-inverse
3286 (function (lambda (x) (math-solve-get-sign x))))
3288 (put 'calcFunc-deg 'math-inverse
3289 (function (lambda (x) (list 'calcFunc-rad x))))
3290 (put 'calcFunc-deg 'math-inverse-sign 1)
3292 (put 'calcFunc-rad 'math-inverse
3293 (function (lambda (x) (list 'calcFunc-deg x))))
3294 (put 'calcFunc-rad 'math-inverse-sign 1)
3296 (put 'calcFunc-ln 'math-inverse
3297 (function (lambda (x) (list 'calcFunc-exp x))))
3298 (put 'calcFunc-ln 'math-inverse-sign 1)
3300 (put 'calcFunc-log10 'math-inverse
3301 (function (lambda (x) (list 'calcFunc-exp10 x))))
3302 (put 'calcFunc-log10 'math-inverse-sign 1)
3304 (put 'calcFunc-lnp1 'math-inverse
3305 (function (lambda (x) (list 'calcFunc-expm1 x))))
3306 (put 'calcFunc-lnp1 'math-inverse-sign 1)
3308 (put 'calcFunc-exp 'math-inverse
3309 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-ln x))
3311 (math-mul '(var pi var-pi)
3313 '(var i var-i))))))))
3314 (put 'calcFunc-exp 'math-inverse-sign 1)
3316 (put 'calcFunc-expm1 'math-inverse
3317 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-lnp1 x))
3319 (math-mul '(var pi var-pi)
3321 '(var i var-i))))))))
3322 (put 'calcFunc-expm1 'math-inverse-sign 1)
3324 (put 'calcFunc-sin 'math-inverse
3325 (function (lambda (x) (let ((n (math-solve-get-int 1)))
3326 (math-add (math-mul (math-normalize
3327 (list 'calcFunc-arcsin x))
3329 (math-mul (math-half-circle t)
3332 (put 'calcFunc-cos 'math-inverse
3333 (function (lambda (x) (math-add (math-solve-get-sign
3335 (list 'calcFunc-arccos x)))
3337 (math-full-circle t))))))
3339 (put 'calcFunc-tan 'math-inverse
3340 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-arctan x))
3342 (math-half-circle t))))))
3344 (put 'calcFunc-arcsin 'math-inverse
3345 (function (lambda (x) (math-normalize (list 'calcFunc-sin x)))))
3347 (put 'calcFunc-arccos 'math-inverse
3348 (function (lambda (x) (math-normalize (list 'calcFunc-cos x)))))
3350 (put 'calcFunc-arctan 'math-inverse
3351 (function (lambda (x) (math-normalize (list 'calcFunc-tan x)))))
3353 (put 'calcFunc-sinh 'math-inverse
3354 (function (lambda (x) (let ((n (math-solve-get-int 1)))
3355 (math-add (math-mul (math-normalize
3356 (list 'calcFunc-arcsinh x))
3358 (math-mul (math-half-circle t)
3362 (put 'calcFunc-sinh 'math-inverse-sign 1)
3364 (put 'calcFunc-cosh 'math-inverse
3365 (function (lambda (x) (math-add (math-solve-get-sign
3367 (list 'calcFunc-arccosh x)))
3368 (math-mul (math-full-circle t)
3370 '(var i var-i)))))))
3372 (put 'calcFunc-tanh 'math-inverse
3373 (function (lambda (x) (math-add (math-normalize
3374 (list 'calcFunc-arctanh x))
3375 (math-mul (math-half-circle t)
3377 '(var i var-i)))))))
3378 (put 'calcFunc-tanh 'math-inverse-sign 1)
3380 (put 'calcFunc-arcsinh 'math-inverse
3381 (function (lambda (x) (math-normalize (list 'calcFunc-sinh x)))))
3382 (put 'calcFunc-arcsinh 'math-inverse-sign 1)
3384 (put 'calcFunc-arccosh 'math-inverse
3385 (function (lambda (x) (math-normalize (list 'calcFunc-cosh x)))))
3387 (put 'calcFunc-arctanh 'math-inverse
3388 (function (lambda (x) (math-normalize (list 'calcFunc-tanh x)))))
3389 (put 'calcFunc-arctanh 'math-inverse-sign 1)
3393 (defun calcFunc-taylor (expr var num)
3394 (let ((x0 0) (v var))
3395 (if (memq (car-safe var) '(+ - calcFunc-eq))
3396 (setq x0 (if (eq (car var) '+) (math-neg (nth 2 var)) (nth 2 var))
3398 (or (and (eq (car-safe v) 'var)
3399 (math-expr-contains expr v)
3401 (let ((accum (math-expr-subst expr v x0))
3402 (var2 (if (eq (car var) 'calcFunc-eq)
3408 (while (and (<= (setq n (1+ n)) num)
3409 (setq fprime (calcFunc-deriv fprime v nil t)))
3410 (setq fprime (math-simplify fprime)
3411 nfac (math-mul nfac n)
3412 accum (math-add accum
3413 (math-div (math-mul (math-pow var2 n)
3418 (math-normalize accum))))
3419 (list 'calcFunc-taylor expr var num))))
3421 ;;; calcalg2.el ends here