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 (defvar math-linear-subst-tried t
742 "Non-nil means that a linear substitution has been tried.")
744 (defun math-do-integral (expr)
745 (let ((math-linear-subst-tried nil
)
747 (or (cond ((not (math-expr-contains expr math-integ-var
))
748 (math-mul expr math-integ-var
))
749 ((equal expr math-integ-var
)
750 (math-div (math-sqr expr
) 2))
752 (and (setq t1
(math-integral (nth 1 expr
)))
753 (setq t2
(math-integral (nth 2 expr
)))
756 (and (setq t1
(math-integral (nth 1 expr
)))
757 (setq t2
(math-integral (nth 2 expr
)))
759 ((eq (car expr
) 'neg
)
760 (and (setq t1
(math-integral (nth 1 expr
)))
763 (cond ((not (math-expr-contains (nth 1 expr
) math-integ-var
))
764 (and (setq t1
(math-integral (nth 2 expr
)))
765 (math-mul (nth 1 expr
) t1
)))
766 ((not (math-expr-contains (nth 2 expr
) math-integ-var
))
767 (and (setq t1
(math-integral (nth 1 expr
)))
768 (math-mul t1
(nth 2 expr
))))
769 ((memq (car-safe (nth 1 expr
)) '(+ -
))
770 (math-integral (list (car (nth 1 expr
))
771 (math-mul (nth 1 (nth 1 expr
))
773 (math-mul (nth 2 (nth 1 expr
))
776 ((memq (car-safe (nth 2 expr
)) '(+ -
))
777 (math-integral (list (car (nth 2 expr
))
778 (math-mul (nth 1 (nth 2 expr
))
780 (math-mul (nth 2 (nth 2 expr
))
784 (cond ((and (not (math-expr-contains (nth 1 expr
)
786 (not (math-equal-int (nth 1 expr
) 1)))
787 (and (setq t1
(math-integral (math-div 1 (nth 2 expr
))))
788 (math-mul (nth 1 expr
) t1
)))
789 ((not (math-expr-contains (nth 2 expr
) math-integ-var
))
790 (and (setq t1
(math-integral (nth 1 expr
)))
791 (math-div t1
(nth 2 expr
))))
792 ((and (eq (car-safe (nth 1 expr
)) '*)
793 (not (math-expr-contains (nth 1 (nth 1 expr
))
795 (and (setq t1
(math-integral
796 (math-div (nth 2 (nth 1 expr
))
798 (math-mul t1
(nth 1 (nth 1 expr
)))))
799 ((and (eq (car-safe (nth 1 expr
)) '*)
800 (not (math-expr-contains (nth 2 (nth 1 expr
))
802 (and (setq t1
(math-integral
803 (math-div (nth 1 (nth 1 expr
))
805 (math-mul t1
(nth 2 (nth 1 expr
)))))
806 ((and (eq (car-safe (nth 2 expr
)) '*)
807 (not (math-expr-contains (nth 1 (nth 2 expr
))
809 (and (setq t1
(math-integral
810 (math-div (nth 1 expr
)
811 (nth 2 (nth 2 expr
)))))
812 (math-div t1
(nth 1 (nth 2 expr
)))))
813 ((and (eq (car-safe (nth 2 expr
)) '*)
814 (not (math-expr-contains (nth 2 (nth 2 expr
))
816 (and (setq t1
(math-integral
817 (math-div (nth 1 expr
)
818 (nth 1 (nth 2 expr
)))))
819 (math-div t1
(nth 2 (nth 2 expr
)))))
820 ((eq (car-safe (nth 2 expr
)) 'calcFunc-exp
)
822 (math-mul (nth 1 expr
)
824 (math-neg (nth 1 (nth 2 expr
)))))))))
826 (cond ((not (math-expr-contains (nth 1 expr
) math-integ-var
))
827 (or (and (setq t1
(math-is-polynomial (nth 2 expr
)
836 (math-mul (nth 2 expr
)
841 ((not (math-expr-contains (nth 2 expr
) math-integ-var
))
842 (if (and (integerp (nth 2 expr
)) (< (nth 2 expr
) 0))
844 (list '/ 1 (math-pow (nth 1 expr
) (- (nth 2 expr
))))
846 (or (and (setq t1
(math-is-polynomial (nth 1 expr
)
849 (setq t2
(math-add (nth 2 expr
) 1))
850 (math-div (math-pow (nth 1 expr
) t2
)
851 (math-mul t2
(nth 1 t1
))))
852 (and (Math-negp (nth 2 expr
))
855 (math-pow (nth 1 expr
)
861 ;; Integral of a polynomial.
862 (and (setq t1
(math-is-polynomial expr math-integ-var
20))
866 (if (setq accum
(math-add accum
867 (math-div (math-mul (car t1
)
876 ;; Try looking it up!
877 (cond ((= (length expr
) 2)
878 (and (symbolp (car expr
))
879 (setq t1
(get (car expr
) 'math-integral
))
882 (not (setq t2
(funcall (car t1
)
885 (and t2
(math-normalize t2
)))))
887 (and (symbolp (car expr
))
888 (setq t1
(get (car expr
) 'math-integral-2
))
891 (not (setq t2
(funcall (car t1
)
895 (and t2
(math-normalize t2
))))))
897 ;; Integral of a rational function.
898 (and (math-ratpoly-p expr math-integ-var
)
899 (setq t1
(calcFunc-apart expr math-integ-var
))
900 (not (equal t1 expr
))
903 ;; Try user-defined integration rules.
905 (let ((math-old-integ (symbol-function 'calcFunc-integ
))
906 (input (list 'calcFunc-integtry expr math-integ-var
))
910 (fset 'calcFunc-integ
'math-sub-integration
)
911 (setq res
(math-rewrite input
912 '(var IntegRules var-IntegRules
)
914 (fset 'calcFunc-integ math-old-integ
)
915 (and (not (equal res input
))
916 (if (setq part
(math-expr-calls
917 res
'(calcFunc-integsubst)))
918 (and (memq (length part
) '(3 4 5))
926 (math-integrate-by-substitution
929 (list 'calcFunc-integfailed
932 (if (not (math-expr-calls res
934 calcFunc-integfailed
)))
936 (fset 'calcFunc-integ math-old-integ
))))
938 ;; See if the function is a symbolic derivative.
939 (and (string-match "'" (symbol-name (car expr
)))
940 (let ((name (symbol-name (car expr
)))
941 (p expr
) (n 0) (which nil
) (bad nil
))
942 (while (setq n
(1+ n
) p
(cdr p
))
943 (if (equal (car p
) math-integ-var
)
944 (if which
(setq bad t
) (setq which n
))
945 (if (math-expr-contains (car p
) math-integ-var
)
948 (let ((prime (if (= which
1) "'" (format "'%d" which
))))
949 (and (string-match (concat prime
"\\('['0-9]*\\|$\\)")
953 (substring name
0 (match-beginning 0))
954 (substring name
(+ (match-beginning 0)
958 ;; Try transformation methods (parts, substitutions).
959 (and (> math-integ-level
0)
960 (math-do-integral-methods expr
))
962 ;; Try expanding the function's definition.
963 (let ((res (math-expand-formula expr
)))
965 (math-integral res
))))))
967 (defun math-sub-integration (expr &rest rest
)
968 (or (if (or (not rest
)
969 (and (< math-integ-level math-integral-limit
)
970 (eq (car rest
) math-integ-var
)))
972 (let ((res (apply math-old-integ expr rest
)))
973 (and (or (= math-integ-level math-integral-limit
)
974 (not (math-expr-calls res
'calcFunc-integ
)))
976 (list 'calcFunc-integfailed expr
)))
978 (defun math-do-integral-methods (expr)
979 (let ((so-far math-integ-var-list-list
)
982 ;; Integration by substitution, for various likely sub-expressions.
983 ;; (In first pass, we look only for sub-exprs that are linear in X.)
984 (or (if math-linear-subst-tried
985 (math-integ-try-substitutions expr
)
986 (math-integ-try-linear-substitutions expr
))
988 ;; If function has sines and cosines, try tan(x/2) substitution.
989 (and (let ((p (setq rat-in
(math-expr-rational-in expr
))))
991 (memq (car (car p
)) '(calcFunc-sin
994 (equal (nth 1 (car p
)) math-integ-var
))
997 (or (and (math-integ-parts-easy expr
)
998 (math-integ-try-parts expr t
))
999 (math-integrate-by-good-substitution
1000 expr
(list 'calcFunc-tan
(math-div math-integ-var
2)))))
1002 ;; If function has sinh and cosh, try tanh(x/2) substitution.
1003 (and (let ((p rat-in
))
1005 (memq (car (car p
)) '(calcFunc-sinh
1009 (equal (nth 1 (car p
)) math-integ-var
))
1012 (or (and (math-integ-parts-easy expr
)
1013 (math-integ-try-parts expr t
))
1014 (math-integrate-by-good-substitution
1015 expr
(list 'calcFunc-tanh
(math-div math-integ-var
2)))))
1017 ;; If function has square roots, try sin, tan, or sec substitution.
1018 (and (let ((p rat-in
))
1021 (or (equal (car p
) math-integ-var
)
1022 (and (eq (car (car p
)) 'calcFunc-sqrt
)
1023 (setq t1
(math-is-polynomial
1024 (nth 1 (setq t2
(car p
)))
1025 math-integ-var
2)))))
1029 (if (math-guess-if-neg (nth 2 t1
))
1030 (let* ((c (math-sqrt (math-neg (nth 2 t1
))))
1031 (d (math-div (nth 1 t1
) (math-mul -
2 c
)))
1032 (a (math-sqrt (math-add (car t1
) (math-sqr d
)))))
1033 (math-integrate-by-good-substitution
1034 expr
(list 'calcFunc-arcsin
1036 (math-add (math-mul c math-integ-var
) d
)
1038 (let* ((c (math-sqrt (nth 2 t1
)))
1039 (d (math-div (nth 1 t1
) (math-mul 2 c
)))
1040 (aa (math-sub (car t1
) (math-sqr d
))))
1041 (if (and nil
(not (and (eq d
0) (eq c
1))))
1042 (math-integrate-by-good-substitution
1043 expr
(math-add (math-mul c math-integ-var
) d
))
1044 (if (math-guess-if-neg aa
)
1045 (math-integrate-by-good-substitution
1046 expr
(list 'calcFunc-arccosh
1048 (math-add (math-mul c math-integ-var
)
1050 (math-sqrt (math-neg aa
)))))
1051 (math-integrate-by-good-substitution
1052 expr
(list 'calcFunc-arcsinh
1054 (math-add (math-mul c math-integ-var
)
1056 (math-sqrt aa
))))))))
1057 (math-integrate-by-good-substitution expr t2
)) )
1059 ;; Try integration by parts.
1060 (math-integ-try-parts expr
)
1065 (defun math-integ-parts-easy (expr)
1066 (cond ((Math-primp expr
) t
)
1067 ((memq (car expr
) '(+ -
*))
1068 (and (math-integ-parts-easy (nth 1 expr
))
1069 (math-integ-parts-easy (nth 2 expr
))))
1071 (and (math-integ-parts-easy (nth 1 expr
))
1072 (math-atomic-factorp (nth 2 expr
))))
1074 (and (natnump (nth 2 expr
))
1075 (math-integ-parts-easy (nth 1 expr
))))
1076 ((eq (car expr
) 'neg
)
1077 (math-integ-parts-easy (nth 1 expr
)))
1080 (defun math-integ-try-parts (expr &optional math-good-parts
)
1081 ;; Integration by parts:
1082 ;; integ(f(x) g(x),x) = f(x) h(x) - integ(h(x) f'(x),x)
1083 ;; where h(x) = integ(g(x),x).
1084 (or (let ((exp (calcFunc-expand expr
)))
1085 (and (not (equal exp expr
))
1086 (math-integral exp
)))
1087 (and (eq (car expr
) '*)
1088 (let ((first-bad (or (math-polynomial-p (nth 1 expr
)
1090 (equal (nth 2 expr
) math-prev-parts-v
))))
1091 (or (and first-bad
; so try this one first
1092 (math-integrate-by-parts (nth 1 expr
) (nth 2 expr
)))
1093 (math-integrate-by-parts (nth 2 expr
) (nth 1 expr
))
1094 (and (not first-bad
)
1095 (math-integrate-by-parts (nth 1 expr
) (nth 2 expr
))))))
1096 (and (eq (car expr
) '/)
1097 (math-expr-contains (nth 1 expr
) math-integ-var
)
1098 (let ((recip (math-div 1 (nth 2 expr
))))
1099 (or (math-integrate-by-parts (nth 1 expr
) recip
)
1100 (math-integrate-by-parts recip
(nth 1 expr
)))))
1101 (and (eq (car expr
) '^
)
1102 (math-integrate-by-parts (math-pow (nth 1 expr
)
1103 (math-sub (nth 2 expr
) 1))
1106 (defun math-integrate-by-parts (u vprime
)
1107 (let ((math-integ-level (if (or math-good-parts
1108 (math-polynomial-p u math-integ-var
))
1110 (1- math-integ-level
)))
1111 (math-doing-parts t
)
1113 (and (>= math-integ-level
0)
1116 (setcar (cdr cur-record
) 'parts
)
1117 (math-tracing-integral "Integrating by parts, u = "
1118 (math-format-value u
1000)
1120 (math-format-value vprime
1000)
1122 (and (setq v
(math-integral vprime
))
1123 (setq temp
(calcFunc-deriv u math-integ-var nil t
))
1124 (setq temp
(let ((math-prev-parts-v v
))
1125 (math-integral (math-mul v temp
) 'yes
)))
1126 (setq temp
(math-sub (math-mul u v
) temp
))
1127 (if (eq (nth 1 cur-record
) 'parts
)
1128 (calcFunc-expand temp
)
1129 (setq v
(list 'var
'PARTS cur-record
)
1130 var-thing
(list 'vec
(math-sub v temp
) v
)
1131 temp
(let (calc-next-why)
1132 (math-solve-for (math-sub v temp
) 0 v nil
)))
1133 (and temp
(not (integerp temp
))
1134 (math-simplify-extended temp
)))))
1135 (setcar (cdr cur-record
) 'busy
)))))
1137 ;;; This tries two different formulations, hoping the algebraic simplifier
1138 ;;; will be strong enough to handle at least one.
1139 (defun math-integrate-by-substitution (expr u
&optional user uinv uinvprime
)
1140 (and (> math-integ-level
0)
1141 (let ((math-integ-level (max (- math-integ-level
2) 0)))
1142 (math-integrate-by-good-substitution expr u user uinv uinvprime
))))
1144 (defun math-integrate-by-good-substitution (expr u
&optional user
1146 (let ((math-living-dangerously t
)
1148 (and (setq uinv
(if uinv
1149 (math-expr-subst uinv math-integ-var
1151 (let (calc-next-why)
1154 math-integ-var nil
))))
1156 (math-tracing-integral "Integrating by substitution, u = "
1157 (math-format-value u
1000)
1159 (or (and (setq deriv
(calcFunc-deriv u
1162 (setq temp
(math-integral (math-expr-subst
1165 (math-div expr deriv
)
1173 (and (setq deriv
(or uinvprime
1174 (calcFunc-deriv uinv
1178 (setq temp
(math-integral (math-mul
1191 (math-simplify-extended
1192 (math-expr-subst temp math-integ-var u
)))))
1194 ;;; Look for substitutions of the form u = a x + b.
1195 (defun math-integ-try-linear-substitutions (sub-expr)
1196 (setq math-linear-subst-tried t
)
1197 (and (not (Math-primp sub-expr
))
1198 (or (and (not (memq (car sub-expr
) '(+ -
* / neg
)))
1199 (not (and (eq (car sub-expr
) '^
)
1200 (integerp (nth 2 sub-expr
))))
1201 (math-expr-contains sub-expr math-integ-var
)
1203 (while (and (setq sub-expr
(cdr sub-expr
))
1204 (or (not (math-linear-in (car sub-expr
)
1206 (assoc (car sub-expr
) so-far
)
1208 (setq so-far
(cons (list (car sub-expr
))
1211 (math-integrate-by-substitution
1212 expr
(car sub-expr
))))))))
1215 (while (and (setq sub-expr
(cdr sub-expr
))
1216 (not (setq res
(math-integ-try-linear-substitutions
1220 ;;; Recursively try different substitutions based on various sub-expressions.
1221 (defun math-integ-try-substitutions (sub-expr &optional allow-rat
)
1222 (and (not (Math-primp sub-expr
))
1223 (not (assoc sub-expr so-far
))
1224 (math-expr-contains sub-expr math-integ-var
)
1225 (or (and (if (and (not (memq (car sub-expr
) '(+ -
* / neg
)))
1226 (not (and (eq (car sub-expr
) '^
)
1227 (integerp (nth 2 sub-expr
)))))
1229 (prog1 allow-rat
(setq allow-rat nil
)))
1230 (not (eq sub-expr expr
))
1231 (or (math-integrate-by-substitution expr sub-expr
)
1232 (and (eq (car sub-expr
) '^
)
1233 (integerp (nth 2 sub-expr
))
1234 (< (nth 2 sub-expr
) 0)
1235 (math-integ-try-substitutions
1236 (math-pow (nth 1 sub-expr
) (- (nth 2 sub-expr
)))
1239 (setq so-far
(cons (list sub-expr
) so-far
))
1240 (while (and (setq sub-expr
(cdr sub-expr
))
1241 (not (setq res
(math-integ-try-substitutions
1242 (car sub-expr
) allow-rat
)))))
1245 (defun math-expr-rational-in (expr)
1247 (math-expr-rational-in-rec expr
)
1248 (mapcar 'car parts
)))
1250 (defun math-expr-rational-in-rec (expr)
1251 (cond ((Math-primp expr
)
1252 (and (equal expr math-integ-var
)
1253 (not (assoc expr parts
))
1254 (setq parts
(cons (list expr
) parts
))))
1255 ((or (memq (car expr
) '(+ -
* / neg
))
1256 (and (eq (car expr
) '^
) (integerp (nth 2 expr
))))
1257 (math-expr-rational-in-rec (nth 1 expr
))
1258 (and (nth 2 expr
) (math-expr-rational-in-rec (nth 2 expr
))))
1259 ((and (eq (car expr
) '^
)
1260 (eq (math-quarter-integer (nth 2 expr
)) 2))
1261 (math-expr-rational-in-rec (list 'calcFunc-sqrt
(nth 1 expr
))))
1263 (and (not (assoc expr parts
))
1264 (math-expr-contains expr math-integ-var
)
1265 (setq parts
(cons (list expr
) parts
))))))
1267 (defun math-expr-calls (expr funcs
&optional arg-contains
)
1269 (if (or (memq (car expr
) funcs
)
1270 (and (eq (car expr
) '^
) (eq (car funcs
) 'calcFunc-sqrt
)
1271 (eq (math-quarter-integer (nth 2 expr
)) 2)))
1272 (and (or (not arg-contains
)
1273 (math-expr-contains expr arg-contains
))
1275 (and (not (Math-primp expr
))
1277 (while (and (setq expr
(cdr expr
))
1278 (not (setq res
(math-expr-calls
1279 (car expr
) funcs arg-contains
)))))
1282 (defun math-fix-const-terms (expr except-vars
)
1283 (cond ((not (math-expr-depends expr except-vars
)) 0)
1284 ((Math-primp expr
) expr
)
1286 (math-add (math-fix-const-terms (nth 1 expr
) except-vars
)
1287 (math-fix-const-terms (nth 2 expr
) except-vars
)))
1289 (math-sub (math-fix-const-terms (nth 1 expr
) except-vars
)
1290 (math-fix-const-terms (nth 2 expr
) except-vars
)))
1293 ;; Command for debugging the Calculator's symbolic integrator.
1294 (defun calc-dump-integral-cache (&optional arg
)
1296 (let ((buf (current-buffer)))
1298 (let ((p math-integral-cache
)
1300 (display-buffer (get-buffer-create "*Integral Cache*"))
1301 (set-buffer (get-buffer "*Integral Cache*"))
1304 (setq cur-record
(car p
))
1305 (or arg
(math-replace-integral-parts cur-record
))
1306 (insert (math-format-flat-expr (car cur-record
) 0)
1308 (if (symbolp (nth 1 cur-record
))
1309 (concat "(" (symbol-name (nth 1 cur-record
)) ")")
1310 (math-format-flat-expr (nth 1 cur-record
) 0))
1313 (goto-char (point-min)))
1316 (defun math-try-integral (expr)
1317 (let ((math-integ-level math-integral-limit
)
1318 (math-integ-depth 0)
1319 (math-integ-msg "Working...done")
1320 (cur-record nil
) ; a technicality
1321 (math-integrating t
)
1322 (calc-prefer-frac t
)
1323 (calc-symbolic-mode t
)
1324 (has-rules (calc-has-rules 'var-IntegRules
)))
1325 (or (math-integral expr
'yes
)
1326 (and math-any-substs
1327 (setq math-enable-subst t
)
1328 (math-integral expr
'yes
))
1329 (and (> math-max-integral-limit math-integral-limit
)
1330 (setq math-integral-limit math-max-integral-limit
1331 math-integ-level math-integral-limit
)
1332 (math-integral expr
'yes
)))))
1334 (defun calcFunc-integ (expr var
&optional low high
)
1336 ;; Do these even if the parts turn out not to be integrable.
1337 ((eq (car-safe expr
) '+)
1338 (math-add (calcFunc-integ (nth 1 expr
) var low high
)
1339 (calcFunc-integ (nth 2 expr
) var low high
)))
1340 ((eq (car-safe expr
) '-
)
1341 (math-sub (calcFunc-integ (nth 1 expr
) var low high
)
1342 (calcFunc-integ (nth 2 expr
) var low high
)))
1343 ((eq (car-safe expr
) 'neg
)
1344 (math-neg (calcFunc-integ (nth 1 expr
) var low high
)))
1345 ((and (eq (car-safe expr
) '*)
1346 (not (math-expr-contains (nth 1 expr
) var
)))
1347 (math-mul (nth 1 expr
) (calcFunc-integ (nth 2 expr
) var low high
)))
1348 ((and (eq (car-safe expr
) '*)
1349 (not (math-expr-contains (nth 2 expr
) var
)))
1350 (math-mul (calcFunc-integ (nth 1 expr
) var low high
) (nth 2 expr
)))
1351 ((and (eq (car-safe expr
) '/)
1352 (not (math-expr-contains (nth 1 expr
) var
))
1353 (not (math-equal-int (nth 1 expr
) 1)))
1354 (math-mul (nth 1 expr
)
1355 (calcFunc-integ (math-div 1 (nth 2 expr
)) var low high
)))
1356 ((and (eq (car-safe expr
) '/)
1357 (not (math-expr-contains (nth 2 expr
) var
)))
1358 (math-div (calcFunc-integ (nth 1 expr
) var low high
) (nth 2 expr
)))
1359 ((and (eq (car-safe expr
) '/)
1360 (eq (car-safe (nth 1 expr
)) '*)
1361 (not (math-expr-contains (nth 1 (nth 1 expr
)) var
)))
1362 (math-mul (nth 1 (nth 1 expr
))
1363 (calcFunc-integ (math-div (nth 2 (nth 1 expr
)) (nth 2 expr
))
1365 ((and (eq (car-safe expr
) '/)
1366 (eq (car-safe (nth 1 expr
)) '*)
1367 (not (math-expr-contains (nth 2 (nth 1 expr
)) var
)))
1368 (math-mul (nth 2 (nth 1 expr
))
1369 (calcFunc-integ (math-div (nth 1 (nth 1 expr
)) (nth 2 expr
))
1371 ((and (eq (car-safe expr
) '/)
1372 (eq (car-safe (nth 2 expr
)) '*)
1373 (not (math-expr-contains (nth 1 (nth 2 expr
)) var
)))
1374 (math-div (calcFunc-integ (math-div (nth 1 expr
) (nth 2 (nth 2 expr
)))
1376 (nth 1 (nth 2 expr
))))
1377 ((and (eq (car-safe expr
) '/)
1378 (eq (car-safe (nth 2 expr
)) '*)
1379 (not (math-expr-contains (nth 2 (nth 2 expr
)) var
)))
1380 (math-div (calcFunc-integ (math-div (nth 1 expr
) (nth 1 (nth 2 expr
)))
1382 (nth 2 (nth 2 expr
))))
1383 ((eq (car-safe expr
) 'vec
)
1384 (cons 'vec
(mapcar (function (lambda (x) (calcFunc-integ x var low high
)))
1387 (let ((state (list calc-angle-mode
1388 ;;calc-symbolic-mode
1391 (calc-var-value 'var-IntegRules
)
1392 (calc-var-value 'var-IntegSimpRules
))))
1393 (or (equal state math-integral-cache-state
)
1394 (setq math-integral-cache-state state
1395 math-integral-cache nil
)))
1396 (let* ((math-max-integral-limit (or (and (boundp 'var-IntegLimit
)
1397 (natnump var-IntegLimit
)
1400 (math-integral-limit 1)
1401 (sexpr (math-expr-subst expr var math-integ-var
))
1402 (trace-buffer (get-buffer "*Trace*"))
1403 (calc-language (if (eq calc-language
'big
) nil calc-language
))
1405 (math-enable-subst nil
)
1406 (math-prev-parts-v nil
)
1407 (math-doing-parts nil
)
1408 (math-good-parts nil
)
1411 (let ((calcbuf (current-buffer))
1412 (calcwin (selected-window)))
1415 (if (get-buffer-window trace-buffer
)
1416 (select-window (get-buffer-window trace-buffer
)))
1417 (set-buffer trace-buffer
)
1418 (goto-char (point-max))
1419 (or (assq 'scroll-stop
(buffer-local-variables))
1421 (make-local-variable 'scroll-step
)
1422 (setq scroll-step
3)))
1424 (set-buffer calcbuf
)
1425 (math-try-integral sexpr
))
1426 (select-window calcwin
)
1427 (set-buffer calcbuf
)))
1428 (math-try-integral sexpr
))))
1431 (if (calc-has-rules 'var-IntegAfterRules
)
1432 (setq res
(math-rewrite res
'(var IntegAfterRules
1433 var-IntegAfterRules
))))
1436 (math-sub (math-expr-subst res math-integ-var high
)
1437 (math-expr-subst res math-integ-var low
))
1438 (setq res
(math-fix-const-terms res math-integ-vars
))
1440 (math-expr-subst res math-integ-var low
)
1441 (math-expr-subst res math-integ-var var
)))))
1442 (append (list 'calcFunc-integ expr var
)
1443 (and low
(list low
))
1444 (and high
(list high
))))))))
1447 (math-defintegral calcFunc-inv
1448 (math-integral (math-div 1 u
)))
1450 (math-defintegral calcFunc-conj
1451 (let ((int (math-integral u
)))
1453 (list 'calcFunc-conj int
))))
1455 (math-defintegral calcFunc-deg
1456 (let ((int (math-integral u
)))
1458 (list 'calcFunc-deg int
))))
1460 (math-defintegral calcFunc-rad
1461 (let ((int (math-integral u
)))
1463 (list 'calcFunc-rad int
))))
1465 (math-defintegral calcFunc-re
1466 (let ((int (math-integral u
)))
1468 (list 'calcFunc-re int
))))
1470 (math-defintegral calcFunc-im
1471 (let ((int (math-integral u
)))
1473 (list 'calcFunc-im int
))))
1475 (math-defintegral calcFunc-sqrt
1476 (and (equal u math-integ-var
)
1477 (math-mul '(frac 2 3)
1478 (list 'calcFunc-sqrt
(math-pow u
3)))))
1480 (math-defintegral calcFunc-exp
1481 (or (and (equal u math-integ-var
)
1482 (list 'calcFunc-exp u
))
1483 (let ((p (math-is-polynomial u math-integ-var
2)))
1485 (let ((sqa (math-sqrt (math-neg (nth 2 p
)))))
1488 (math-mul (math-div (list 'calcFunc-sqrt
'(var pi var-pi
))
1492 (math-div (math-sub (math-mul (car p
)
1495 (math-sqr (nth 1 p
))
1499 (math-sub (math-mul sqa math-integ-var
)
1500 (math-div (nth 1 p
) (math-mul 2 sqa
)))))
1503 (math-defintegral calcFunc-ln
1504 (or (and (equal u math-integ-var
)
1505 (math-sub (math-mul u
(list 'calcFunc-ln u
)) u
))
1506 (and (eq (car u
) '*)
1507 (math-integral (math-add (list 'calcFunc-ln
(nth 1 u
))
1508 (list 'calcFunc-ln
(nth 2 u
)))))
1509 (and (eq (car u
) '/)
1510 (math-integral (math-sub (list 'calcFunc-ln
(nth 1 u
))
1511 (list 'calcFunc-ln
(nth 2 u
)))))
1512 (and (eq (car u
) '^
)
1513 (math-integral (math-mul (nth 2 u
)
1514 (list 'calcFunc-ln
(nth 1 u
)))))))
1516 (math-defintegral calcFunc-log10
1517 (and (equal u math-integ-var
)
1518 (math-sub (math-mul u
(list 'calcFunc-ln u
))
1519 (math-div u
(list 'calcFunc-ln
10)))))
1521 (math-defintegral-2 calcFunc-log
1522 (math-integral (math-div (list 'calcFunc-ln u
)
1523 (list 'calcFunc-ln v
))))
1525 (math-defintegral calcFunc-sin
1526 (or (and (equal u math-integ-var
)
1527 (math-neg (math-from-radians-2 (list 'calcFunc-cos u
))))
1528 (and (nth 2 (math-is-polynomial u math-integ-var
2))
1529 (math-integral (math-to-exponentials (list 'calcFunc-sin u
))))))
1531 (math-defintegral calcFunc-cos
1532 (or (and (equal u math-integ-var
)
1533 (math-from-radians-2 (list 'calcFunc-sin u
)))
1534 (and (nth 2 (math-is-polynomial u math-integ-var
2))
1535 (math-integral (math-to-exponentials (list 'calcFunc-cos u
))))))
1537 (math-defintegral calcFunc-tan
1538 (and (equal u math-integ-var
)
1539 (math-neg (math-from-radians-2
1540 (list 'calcFunc-ln
(list 'calcFunc-cos u
))))))
1542 (math-defintegral calcFunc-arcsin
1543 (and (equal u math-integ-var
)
1544 (math-add (math-mul u
(list 'calcFunc-arcsin u
))
1545 (math-from-radians-2
1546 (list 'calcFunc-sqrt
(math-sub 1 (math-sqr u
)))))))
1548 (math-defintegral calcFunc-arccos
1549 (and (equal u math-integ-var
)
1550 (math-sub (math-mul u
(list 'calcFunc-arccos u
))
1551 (math-from-radians-2
1552 (list 'calcFunc-sqrt
(math-sub 1 (math-sqr u
)))))))
1554 (math-defintegral calcFunc-arctan
1555 (and (equal u math-integ-var
)
1556 (math-sub (math-mul u
(list 'calcFunc-arctan u
))
1557 (math-from-radians-2
1558 (math-div (list 'calcFunc-ln
(math-add 1 (math-sqr u
)))
1561 (math-defintegral calcFunc-sinh
1562 (and (equal u math-integ-var
)
1563 (list 'calcFunc-cosh u
)))
1565 (math-defintegral calcFunc-cosh
1566 (and (equal u math-integ-var
)
1567 (list 'calcFunc-sinh u
)))
1569 (math-defintegral calcFunc-tanh
1570 (and (equal u math-integ-var
)
1571 (list 'calcFunc-ln
(list 'calcFunc-cosh u
))))
1573 (math-defintegral calcFunc-arcsinh
1574 (and (equal u math-integ-var
)
1575 (math-sub (math-mul u
(list 'calcFunc-arcsinh u
))
1576 (list 'calcFunc-sqrt
(math-add (math-sqr u
) 1)))))
1578 (math-defintegral calcFunc-arccosh
1579 (and (equal u math-integ-var
)
1580 (math-sub (math-mul u
(list 'calcFunc-arccosh u
))
1581 (list 'calcFunc-sqrt
(math-sub 1 (math-sqr u
))))))
1583 (math-defintegral calcFunc-arctanh
1584 (and (equal u math-integ-var
)
1585 (math-sub (math-mul u
(list 'calcFunc-arctan u
))
1586 (math-div (list 'calcFunc-ln
1587 (math-add 1 (math-sqr u
)))
1590 ;;; (Ax + B) / (ax^2 + bx + c)^n forms.
1591 (math-defintegral-2 /
1592 (math-integral-rational-funcs u v
))
1594 (defun math-integral-rational-funcs (u v
)
1595 (let ((pu (math-is-polynomial u math-integ-var
1))
1599 (if (and (eq (car-safe v
) '^
) (natnump (nth 2 v
)))
1600 (setq vpow
(nth 2 v
)
1602 (and (setq pv
(math-is-polynomial v math-integ-var
2))
1603 (let ((int (math-mul-thru
1605 (math-integral-q02 (car pv
) (nth 1 pv
)
1606 (nth 2 pv
) v vpow
))))
1608 (setq int
(math-add int
1613 (nth 2 pv
) v vpow
)))))
1616 (defun math-integral-q12 (a b c v vpow
)
1620 (math-sub (math-div math-integ-var b
)
1621 (math-mul (math-div a
(math-sqr b
))
1622 (list 'calcFunc-ln v
))))
1624 (math-div (math-add (list 'calcFunc-ln v
)
1628 (let ((nm1 (math-sub vpow
1))
1629 (nm2 (math-sub vpow
2)))
1631 (math-div a
(math-mul nm1
(math-pow v nm1
)))
1632 (math-div 1 (math-mul nm2
(math-pow v nm2
))))
1635 (setq q
(math-sub (math-mul 4 (math-mul a c
)) (math-sqr b
))))
1636 (let ((part (math-div b
(math-mul 2 c
))))
1637 (math-mul-thru (math-pow c vpow
)
1638 (math-integral-q12 part
1 nil
1639 (math-add math-integ-var part
)
1642 (and (math-ratp q
) (math-negp q
)
1643 (let ((calc-symbolic-mode t
))
1644 (math-ratp (math-sqrt (math-neg q
))))
1645 (throw 'int-rat nil
)) ; should have used calcFunc-apart first
1646 (math-sub (math-div (list 'calcFunc-ln v
) (math-mul 2 c
))
1647 (math-mul-thru (math-div b
(math-mul 2 c
))
1648 (math-integral-q02 a b c v
1))))
1650 (let ((n (1- vpow
)))
1651 (math-sub (math-neg (math-div
1652 (math-add (math-mul b math-integ-var
)
1654 (math-mul n
(math-mul q
(math-pow v n
)))))
1655 (math-mul-thru (math-div (math-mul b
(1- (* 2 n
)))
1657 (math-integral-q02 a b c v n
))))))))
1659 (defun math-integral-q02 (a b c v vpow
)
1663 (math-div (list 'calcFunc-ln v
) b
))
1665 (math-div (math-pow v
(- 1 vpow
))
1666 (math-mul (- 1 vpow
) b
)))))
1668 (setq q
(math-sub (math-mul 4 (math-mul a c
)) (math-sqr b
))))
1669 (let ((part (math-div b
(math-mul 2 c
))))
1670 (math-mul-thru (math-pow c vpow
)
1671 (math-integral-q02 part
1 nil
1672 (math-add math-integ-var part
)
1675 (setq part
(math-add (math-mul 2 (math-mul c math-integ-var
)) b
))
1677 (let ((n (1- vpow
)))
1678 (math-add (math-div part
(math-mul n
(math-mul q
(math-pow v n
))))
1679 (math-mul-thru (math-div (math-mul (- (* 4 n
) 2) c
)
1681 (math-integral-q02 a b c v n
)))))
1682 ((math-guess-if-neg q
)
1683 (setq rq
(list 'calcFunc-sqrt
(math-neg q
)))
1684 ;;(math-div-thru (list 'calcFunc-ln
1685 ;; (math-div (math-sub part rq)
1686 ;; (math-add part rq)))
1688 (math-div (math-mul -
2 (list 'calcFunc-arctanh
1689 (math-div part rq
)))
1692 (setq rq
(list 'calcFunc-sqrt q
))
1693 (math-div (math-mul 2 (math-to-radians-2
1694 (list 'calcFunc-arctan
1695 (math-div part rq
))))
1699 (math-defintegral calcFunc-erf
1700 (and (equal u math-integ-var
)
1701 (math-add (math-mul u
(list 'calcFunc-erf u
))
1702 (math-div 1 (math-mul (list 'calcFunc-exp
(math-sqr u
))
1703 (list 'calcFunc-sqrt
1704 '(var pi var-pi
)))))))
1706 (math-defintegral calcFunc-erfc
1707 (and (equal u math-integ-var
)
1708 (math-sub (math-mul u
(list 'calcFunc-erfc u
))
1709 (math-div 1 (math-mul (list 'calcFunc-exp
(math-sqr u
))
1710 (list 'calcFunc-sqrt
1711 '(var pi var-pi
)))))))
1716 (defvar math-tabulate-initial nil
)
1717 (defvar math-tabulate-function nil
)
1718 (defun calcFunc-table (expr var
&optional low high step
)
1719 (or low
(setq low
'(neg (var inf var-inf
)) high
'(var inf var-inf
)))
1720 (or high
(setq high low low
1))
1721 (and (or (math-infinitep low
) (math-infinitep high
))
1723 (math-scan-for-limits expr
))
1724 (and step
(math-zerop step
) (math-reject-arg step
'nonzerop
))
1725 (let ((known (+ (if (Math-objectp low
) 1 0)
1726 (if (Math-objectp high
) 1 0)
1727 (if (or (null step
) (Math-objectp step
)) 1 0)))
1728 (count '(var inf var-inf
))
1730 (or (= known
2) ; handy optimization
1731 (equal high
'(var inf var-inf
))
1733 (setq count
(math-div (math-sub high low
) (or step
1)))
1734 (or (Math-objectp count
)
1735 (setq count
(math-simplify count
)))
1736 (if (Math-messy-integerp count
)
1737 (setq count
(math-trunc count
)))))
1738 (if (Math-negp count
)
1740 (if (integerp count
)
1741 (let ((var-DUMMY nil
)
1742 (vec math-tabulate-initial
)
1743 (math-working-step-2 (1+ count
))
1744 (math-working-step 0))
1745 (setq expr
(math-evaluate-expr
1746 (math-expr-subst expr var
'(var DUMMY var-DUMMY
))))
1748 (setq math-working-step
(1+ math-working-step
)
1750 vec
(cond ((eq math-tabulate-function
'calcFunc-sum
)
1751 (math-add vec
(math-evaluate-expr expr
)))
1752 ((eq math-tabulate-function
'calcFunc-prod
)
1753 (math-mul vec
(math-evaluate-expr expr
)))
1755 (cons (math-evaluate-expr expr
) vec
)))
1756 low
(math-add low
(or step
1))
1758 (if math-tabulate-function
1760 (cons 'vec
(nreverse vec
))))
1761 (if (Math-integerp count
)
1762 (calc-record-why 'fixnump high
)
1763 (if (Math-num-integerp low
)
1764 (if (Math-num-integerp high
)
1765 (calc-record-why 'integerp step
)
1766 (calc-record-why 'integerp high
))
1767 (calc-record-why 'integerp low
)))
1768 (append (list (or math-tabulate-function
'calcFunc-table
)
1770 (and (not (and (equal low
'(neg (var inf var-inf
)))
1771 (equal high
'(var inf var-inf
))))
1773 (and step
(list step
))))))
1775 (defun math-scan-for-limits (x)
1776 (cond ((Math-primp x
))
1777 ((and (eq (car x
) 'calcFunc-subscr
)
1778 (Math-vectorp (nth 1 x
))
1779 (math-expr-contains (nth 2 x
) var
))
1780 (let* ((calc-next-why nil
)
1781 (low-val (math-solve-for (nth 2 x
) 1 var nil
))
1782 (high-val (math-solve-for (nth 2 x
) (1- (length (nth 1 x
)))
1785 (and low-val
(math-realp low-val
)
1786 high-val
(math-realp high-val
))
1787 (and (Math-lessp high-val low-val
)
1788 (setq temp low-val low-val high-val high-val temp
))
1789 (setq low
(math-max low
(math-ceiling low-val
))
1790 high
(math-min high
(math-floor high-val
)))))
1792 (while (setq x
(cdr x
))
1793 (math-scan-for-limits (car x
))))))
1796 (defvar math-disable-sums nil
)
1797 (defun calcFunc-sum (expr var
&optional low high step
)
1798 (if math-disable-sums
(math-reject-arg))
1799 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec
2)))
1800 (math-sum-rec expr var low high step
)))
1801 (math-disable-sums t
))
1802 (math-normalize res
)))
1804 (defun math-sum-rec (expr var
&optional low high step
)
1805 (or low
(setq low
'(neg (var inf var-inf
)) high
'(var inf var-inf
)))
1806 (and low
(not high
) (setq high low low
1))
1810 ((not (math-expr-contains expr var
))
1811 (math-mul expr
(math-add (math-div (math-sub high low
) (or step
1))
1813 ((and step
(not (math-equal-int step
1)))
1814 (if (math-negp step
)
1815 (math-sum-rec expr var high low
(math-neg step
))
1816 (let ((lo (math-simplify (math-div low step
))))
1817 (if (math-known-num-integerp lo
)
1818 (math-sum-rec (math-normalize
1819 (math-expr-subst expr var
1820 (math-mul step var
)))
1821 var lo
(math-simplify (math-div high step
)))
1822 (math-sum-rec (math-normalize
1823 (math-expr-subst expr var
1824 (math-add (math-mul step var
)
1827 (math-simplify (math-div (math-sub high low
)
1829 ((memq (setq t1
(math-compare low high
)) '(0 1))
1831 (math-expr-subst expr var low
)
1833 ((setq t1
(math-is-polynomial expr var
20))
1837 (setq poly
(math-poly-mix poly
1
1838 (math-sum-integer-power n
) (car t1
))
1841 (setq n
(math-build-polynomial-expr poly high
))
1842 (if (memq low
'(0 1))
1844 (math-sub n
(math-build-polynomial-expr poly
1845 (math-sub low
1))))))
1846 ((and (memq (car expr
) '(+ -
))
1847 (setq t1
(math-sum-rec (nth 1 expr
) var low high
)
1848 t2
(math-sum-rec (nth 2 expr
) var low high
))
1849 (not (and (math-expr-calls t1
'(calcFunc-sum))
1850 (math-expr-calls t2
'(calcFunc-sum)))))
1851 (list (car expr
) t1 t2
))
1852 ((and (eq (car expr
) '*)
1853 (setq t1
(math-sum-const-factors expr var
)))
1854 (math-mul (car t1
) (math-sum-rec (cdr t1
) var low high
)))
1855 ((and (eq (car expr
) '*) (memq (car-safe (nth 1 expr
)) '(+ -
)))
1856 (math-sum-rec (math-add-or-sub (math-mul (nth 1 (nth 1 expr
))
1858 (math-mul (nth 2 (nth 1 expr
))
1860 nil
(eq (car (nth 1 expr
)) '-
))
1862 ((and (eq (car expr
) '*) (memq (car-safe (nth 2 expr
)) '(+ -
)))
1863 (math-sum-rec (math-add-or-sub (math-mul (nth 1 expr
)
1864 (nth 1 (nth 2 expr
)))
1865 (math-mul (nth 1 expr
)
1866 (nth 2 (nth 2 expr
)))
1867 nil
(eq (car (nth 2 expr
)) '-
))
1869 ((and (eq (car expr
) '/)
1870 (not (math-primp (nth 1 expr
)))
1871 (setq t1
(math-sum-const-factors (nth 1 expr
) var
)))
1873 (math-sum-rec (math-div (cdr t1
) (nth 2 expr
))
1875 ((and (eq (car expr
) '/)
1876 (setq t1
(math-sum-const-factors (nth 2 expr
) var
)))
1877 (math-div (math-sum-rec (math-div (nth 1 expr
) (cdr t1
))
1880 ((eq (car expr
) 'neg
)
1881 (math-neg (math-sum-rec (nth 1 expr
) var low high
)))
1882 ((and (eq (car expr
) '^
)
1883 (not (math-expr-contains (nth 1 expr
) var
))
1884 (setq t1
(math-is-polynomial (nth 2 expr
) var
1)))
1885 (let ((x (math-pow (nth 1 expr
) (nth 1 t1
))))
1886 (math-div (math-mul (math-sub (math-pow x
(math-add 1 high
))
1888 (math-pow (nth 1 expr
) (car t1
)))
1890 ((and (setq t1
(math-to-exponentials expr
))
1891 (setq t1
(math-sum-rec t1 var low high
))
1892 (not (math-expr-calls t1
'(calcFunc-sum))))
1894 ((memq (car expr
) '(calcFunc-ln calcFunc-log10
))
1895 (list (car expr
) (calcFunc-prod (nth 1 expr
) var low high
)))
1896 ((and (eq (car expr
) 'calcFunc-log
)
1898 (not (math-expr-contains (nth 2 expr
) var
)))
1900 (calcFunc-prod (nth 1 expr
) var low high
)
1902 (if (equal val
'(var nan var-nan
)) (setq val nil
))
1904 (let* ((math-tabulate-initial 0)
1905 (math-tabulate-function 'calcFunc-sum
))
1906 (calcFunc-table expr var low high
)))))
1908 (defun calcFunc-asum (expr var low
&optional high step no-mul-flag
)
1909 (or high
(setq high low low
1))
1910 (if (and step
(not (math-equal-int step
1)))
1911 (if (math-negp step
)
1912 (math-mul (math-pow -
1 low
)
1913 (calcFunc-asum expr var high low
(math-neg step
) t
))
1914 (let ((lo (math-simplify (math-div low step
))))
1915 (if (math-num-integerp lo
)
1916 (calcFunc-asum (math-normalize
1917 (math-expr-subst expr var
1918 (math-mul step var
)))
1919 var lo
(math-simplify (math-div high step
)))
1920 (calcFunc-asum (math-normalize
1921 (math-expr-subst expr var
1922 (math-add (math-mul step var
)
1925 (math-simplify (math-div (math-sub high low
)
1927 (math-mul (if no-mul-flag
1 (math-pow -
1 low
))
1928 (calcFunc-sum (math-mul (math-pow -
1 var
) expr
) var low high
))))
1930 (defun math-sum-const-factors (expr var
)
1934 (while (eq (car-safe p
) '*)
1935 (if (math-expr-contains (nth 1 p
) var
)
1936 (setq not-const
(cons (nth 1 p
) not-const
))
1937 (setq const
(cons (nth 1 p
) const
)))
1939 (if (math-expr-contains p var
)
1940 (setq not-const
(cons p not-const
))
1941 (setq const
(cons p const
)))
1943 (cons (let ((temp (car const
)))
1944 (while (setq const
(cdr const
))
1945 (setq temp
(list '* (car const
) temp
)))
1947 (let ((temp (or (car not-const
) 1)))
1948 (while (setq not-const
(cdr not-const
))
1949 (setq temp
(list '* (car not-const
) temp
)))
1952 (defvar math-sum-int-pow-cache
(list '(0 1)))
1953 ;; Following is from CRC Math Tables, 27th ed, pp. 52-53.
1954 (defun math-sum-integer-power (pow)
1955 (let ((calc-prefer-frac t
)
1956 (n (length math-sum-int-pow-cache
)))
1958 (let* ((new (list 0 0))
1960 (pp (cdr (nth (1- n
) math-sum-int-pow-cache
)))
1965 (setq q
(math-div (car pp
) p
)
1966 new
(cons (math-mul q n
) new
)
1967 sum
(math-add sum q
)
1970 (setcar lin
(math-sub 1 (math-mul n sum
)))
1971 (setq math-sum-int-pow-cache
1972 (nconc math-sum-int-pow-cache
(list (nreverse new
)))
1974 (nth pow math-sum-int-pow-cache
)))
1976 (defun math-to-exponentials (expr)
1979 (let ((x (nth 1 expr
))
1980 (pi (if calc-symbolic-mode
'(var pi var-pi
) (math-pi)))
1981 (i (if calc-symbolic-mode
'(var i var-i
) '(cplx 0 1))))
1982 (cond ((eq (car expr
) 'calcFunc-exp
)
1983 (list '^
'(var e var-e
) x
))
1984 ((eq (car expr
) 'calcFunc-sin
)
1985 (or (eq calc-angle-mode
'rad
)
1986 (setq x
(list '/ (list '* x pi
) 180)))
1988 (list '^
'(var e var-e
) (list '* x i
))
1989 (list '^
'(var e var-e
)
1990 (list 'neg
(list '* x i
))))
1992 ((eq (car expr
) 'calcFunc-cos
)
1993 (or (eq calc-angle-mode
'rad
)
1994 (setq x
(list '/ (list '* x pi
) 180)))
1996 (list '^
'(var e var-e
)
1998 (list '^
'(var e var-e
)
1999 (list 'neg
(list '* x i
))))
2001 ((eq (car expr
) 'calcFunc-sinh
)
2003 (list '^
'(var e var-e
) x
)
2004 (list '^
'(var e var-e
) (list 'neg x
)))
2006 ((eq (car expr
) 'calcFunc-cosh
)
2008 (list '^
'(var e var-e
) x
)
2009 (list '^
'(var e var-e
) (list 'neg x
)))
2013 (defun math-to-exps (expr)
2014 (cond (calc-symbolic-mode expr
)
2016 (if (equal expr
'(var e var-e
)) (math-e) expr
))
2017 ((and (eq (car expr
) '^
)
2018 (equal (nth 1 expr
) '(var e var-e
)))
2019 (list 'calcFunc-exp
(nth 2 expr
)))
2021 (cons (car expr
) (mapcar 'math-to-exps
(cdr expr
))))))
2024 (defvar math-disable-prods nil
)
2025 (defun calcFunc-prod (expr var
&optional low high step
)
2026 (if math-disable-prods
(math-reject-arg))
2027 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec
2)))
2028 (math-prod-rec expr var low high step
)))
2029 (math-disable-prods t
))
2030 (math-normalize res
)))
2032 (defun math-prod-rec (expr var
&optional low high step
)
2033 (or low
(setq low
'(neg (var inf var-inf
)) high
'(var inf var-inf
)))
2034 (and low
(not high
) (setq high
'(var inf var-inf
)))
2038 ((not (math-expr-contains expr var
))
2039 (math-pow expr
(math-add (math-div (math-sub high low
) (or step
1))
2041 ((and step
(not (math-equal-int step
1)))
2042 (if (math-negp step
)
2043 (math-prod-rec expr var high low
(math-neg step
))
2044 (let ((lo (math-simplify (math-div low step
))))
2045 (if (math-known-num-integerp lo
)
2046 (math-prod-rec (math-normalize
2047 (math-expr-subst expr var
2048 (math-mul step var
)))
2049 var lo
(math-simplify (math-div high step
)))
2050 (math-prod-rec (math-normalize
2051 (math-expr-subst expr var
2052 (math-add (math-mul step
2056 (math-simplify (math-div (math-sub high low
)
2058 ((and (memq (car expr
) '(* /))
2059 (setq t1
(math-prod-rec (nth 1 expr
) var low high
)
2060 t2
(math-prod-rec (nth 2 expr
) var low high
))
2061 (not (and (math-expr-calls t1
'(calcFunc-prod))
2062 (math-expr-calls t2
'(calcFunc-prod)))))
2063 (list (car expr
) t1 t2
))
2064 ((and (eq (car expr
) '^
)
2065 (not (math-expr-contains (nth 2 expr
) var
)))
2066 (math-pow (math-prod-rec (nth 1 expr
) var low high
)
2068 ((and (eq (car expr
) '^
)
2069 (not (math-expr-contains (nth 1 expr
) var
)))
2070 (math-pow (nth 1 expr
)
2071 (calcFunc-sum (nth 2 expr
) var low high
)))
2072 ((eq (car expr
) 'sqrt
)
2073 (math-normalize (list 'calcFunc-sqrt
2074 (list 'calcFunc-prod
(nth 1 expr
)
2076 ((eq (car expr
) 'neg
)
2077 (math-mul (math-pow -
1 (math-add (math-sub high low
) 1))
2078 (math-prod-rec (nth 1 expr
) var low high
)))
2079 ((eq (car expr
) 'calcFunc-exp
)
2080 (list 'calcFunc-exp
(calcFunc-sum (nth 1 expr
) var low high
)))
2081 ((and (setq t1
(math-is-polynomial expr var
1))
2084 ((or (and (math-equal-int (nth 1 t1
) 1)
2085 (setq low
(math-simplify
2086 (math-add low
(car t1
)))
2088 (math-add high
(car t1
)))))
2089 (and (math-equal-int (nth 1 t1
) -
1)
2092 (math-sub (car t1
) high
))
2094 (math-sub (car t1
) t2
)))))
2095 (if (or (math-zerop low
) (math-zerop high
))
2097 (if (and (or (math-negp low
) (math-negp high
))
2098 (or (math-num-integerp low
)
2099 (math-num-integerp high
)))
2100 (if (math-posp high
)
2102 (math-mul (math-pow -
1
2104 (math-add low high
) 1))
2106 (list 'calcFunc-fact
2108 (list 'calcFunc-fact
2109 (math-sub -
1 high
)))))
2111 (list 'calcFunc-fact high
)
2112 (list 'calcFunc-fact
(math-sub low
1))))))
2113 ((and (or (and (math-equal-int (nth 1 t1
) 2)
2114 (setq t2
(math-simplify
2115 (math-add (math-mul low
2)
2118 (math-add (math-mul high
2)
2120 (and (math-equal-int (nth 1 t1
) -
2)
2121 (setq t2
(math-simplify
2128 (or (math-integerp t2
)
2129 (and (math-messy-integerp t2
)
2130 (setq t2
(math-trunc t2
)))
2132 (and (math-messy-integerp t3
)
2133 (setq t3
(math-trunc t3
)))))
2134 (if (or (math-zerop t2
) (math-zerop t3
))
2136 (if (or (math-evenp t2
) (math-evenp t3
))
2137 (if (or (math-negp t2
) (math-negp t3
))
2138 (if (math-posp high
)
2141 (list 'calcFunc-dfact
2143 (list 'calcFunc-dfact
2146 (list 'calcFunc-dfact t3
)
2147 (list 'calcFunc-dfact
2152 (list '/ (list '-
(list '- t2 t3
)
2156 (list 'calcFunc-dfact
2158 (list 'calcFunc-dfact
2162 (list 'calcFunc-dfact t3
)
2163 (list 'calcFunc-dfact
2167 (if (equal val
'(var nan var-nan
)) (setq val nil
))
2169 (let* ((math-tabulate-initial 1)
2170 (math-tabulate-function 'calcFunc-prod
))
2171 (calcFunc-table expr var low high
)))))
2176 (defvar math-solve-ranges nil
)
2177 ;;; Attempt to reduce lhs = rhs to solve-var = rhs', where solve-var appears
2178 ;;; in lhs but not in rhs or rhs'; return rhs'.
2179 ;;; Uses global values: solve-*.
2180 (defun math-try-solve-for (lhs rhs
&optional sign no-poly
)
2182 (cond ((equal lhs solve-var
)
2183 (setq math-solve-sign sign
)
2184 (if (eq solve-full
'all
)
2185 (let ((vec (list 'vec
(math-evaluate-expr rhs
)))
2187 (while math-solve-ranges
2188 (setq p
(car math-solve-ranges
)
2191 (while (setq p
(cdr p
))
2192 (setq newvec
(nconc newvec
2193 (cdr (math-expr-subst
2194 vec var
(car p
))))))
2196 math-solve-ranges
(cdr math-solve-ranges
)))
2197 (math-normalize vec
))
2201 ((and (eq (car lhs
) '-
)
2202 (eq (car-safe (nth 1 lhs
)) (car-safe (nth 2 lhs
)))
2204 (= (length (nth 1 lhs
)) 2)
2205 (= (length (nth 2 lhs
)) 2)
2206 (setq t1
(get (car (nth 1 lhs
)) 'math-inverse
))
2207 (setq t2
(funcall t1
'(var SOLVEDUM SOLVEDUM
)))
2208 (eq (math-expr-contains-count t2
'(var SOLVEDUM SOLVEDUM
)) 1)
2209 (setq t3
(math-solve-above-dummy t2
))
2210 (setq t1
(math-try-solve-for (math-sub (nth 1 (nth 1 lhs
))
2213 (nth 1 (nth 2 lhs
))))
2216 ((eq (car lhs
) 'neg
)
2217 (math-try-solve-for (nth 1 lhs
) (math-neg rhs
)
2218 (and sign
(- sign
))))
2219 ((and (not (eq solve-full
't
)) (math-try-solve-prod)))
2221 (setq t2
(math-decompose-poly lhs solve-var
15 rhs
)))
2222 (setq t1
(cdr (nth 1 t2
))
2223 t1
(let ((math-solve-ranges math-solve-ranges
))
2224 (cond ((= (length t1
) 5)
2225 (apply 'math-solve-quartic
(car t2
) t1
))
2227 (apply 'math-solve-cubic
(car t2
) t1
))
2229 (apply 'math-solve-quadratic
(car t2
) t1
))
2231 (apply 'math-solve-linear
(car t2
) sign t1
))
2233 (math-poly-all-roots (car t2
) t1
))
2234 (calc-symbolic-mode nil
)
2238 (math-poly-any-root (reverse t1
) 0 t
)
2241 (if (eq (nth 2 t2
) 1)
2243 (math-solve-prod t1
(math-try-solve-for (nth 2 t2
) 0 nil t
)))
2244 (calc-record-why "*Unable to find a symbolic solution")
2246 ((and (math-solve-find-root-term lhs nil
)
2247 (eq (math-expr-contains-count lhs t1
) 1)) ; just in case
2248 (math-try-solve-for (math-simplify
2249 (math-sub (if (or t3
(math-evenp t2
))
2251 (math-neg (math-pow t1 t2
)))
2253 (math-sub (math-normalize
2260 (cond ((not (math-expr-contains (nth 1 lhs
) solve-var
))
2261 (math-try-solve-for (nth 2 lhs
)
2262 (math-sub rhs
(nth 1 lhs
))
2264 ((not (math-expr-contains (nth 2 lhs
) solve-var
))
2265 (math-try-solve-for (nth 1 lhs
)
2266 (math-sub rhs
(nth 2 lhs
))
2268 ((eq (car lhs
) 'calcFunc-eq
)
2269 (math-try-solve-for (math-sub (nth 1 lhs
) (nth 2 lhs
))
2272 (cond ((or (and (eq (car-safe (nth 1 lhs
)) 'calcFunc-sin
)
2273 (eq (car-safe (nth 2 lhs
)) 'calcFunc-cos
))
2274 (and (eq (car-safe (nth 1 lhs
)) 'calcFunc-cos
)
2275 (eq (car-safe (nth 2 lhs
)) 'calcFunc-sin
)))
2276 (math-try-solve-for (math-sub (nth 1 lhs
)
2277 (list (car (nth 1 lhs
))
2279 (math-quarter-circle t
)
2280 (nth 1 (nth 2 lhs
)))))
2282 ((not (math-expr-contains (nth 1 lhs
) solve-var
))
2283 (math-try-solve-for (nth 2 lhs
)
2284 (math-sub (nth 1 lhs
) rhs
)
2285 (and sign
(- sign
))))
2286 ((not (math-expr-contains (nth 2 lhs
) solve-var
))
2287 (math-try-solve-for (nth 1 lhs
)
2288 (math-add rhs
(nth 2 lhs
))
2290 ((and (eq solve-full
't
) (math-try-solve-prod)))
2291 ((and (eq (car lhs
) '%
)
2292 (not (math-expr-contains (nth 2 lhs
) solve-var
)))
2293 (math-try-solve-for (nth 1 lhs
) (math-add rhs
2296 ((eq (car lhs
) 'calcFunc-log
)
2297 (cond ((not (math-expr-contains (nth 2 lhs
) solve-var
))
2298 (math-try-solve-for (nth 1 lhs
) (math-pow (nth 2 lhs
) rhs
)))
2299 ((not (math-expr-contains (nth 1 lhs
) solve-var
))
2300 (math-try-solve-for (nth 2 lhs
) (math-pow
2302 (math-div 1 rhs
))))))
2303 ((and (= (length lhs
) 2)
2305 (setq t1
(get (car lhs
) 'math-inverse
))
2306 (setq t2
(funcall t1 rhs
)))
2307 (setq t1
(get (car lhs
) 'math-inverse-sign
))
2308 (math-try-solve-for (nth 1 lhs
) (math-normalize t2
)
2312 (funcall t1 lhs sign
)))))
2313 ((and (symbolp (car lhs
))
2314 (setq t1
(get (car lhs
) 'math-inverse-n
))
2315 (setq t2
(funcall t1 lhs rhs
)))
2317 ((setq t1
(math-expand-formula lhs
))
2318 (math-try-solve-for t1 rhs sign
))
2320 (calc-record-why "*No inverse known" lhs
)
2324 (defun math-try-solve-prod ()
2325 (cond ((eq (car lhs
) '*)
2326 (cond ((not (math-expr-contains (nth 1 lhs
) solve-var
))
2327 (math-try-solve-for (nth 2 lhs
)
2328 (math-div rhs
(nth 1 lhs
))
2329 (math-solve-sign sign
(nth 1 lhs
))))
2330 ((not (math-expr-contains (nth 2 lhs
) solve-var
))
2331 (math-try-solve-for (nth 1 lhs
)
2332 (math-div rhs
(nth 2 lhs
))
2333 (math-solve-sign sign
(nth 2 lhs
))))
2335 (math-solve-prod (let ((math-solve-ranges math-solve-ranges
))
2336 (math-try-solve-for (nth 2 lhs
) 0))
2337 (math-try-solve-for (nth 1 lhs
) 0)))))
2339 (cond ((not (math-expr-contains (nth 1 lhs
) solve-var
))
2340 (math-try-solve-for (nth 2 lhs
)
2341 (math-div (nth 1 lhs
) rhs
)
2342 (math-solve-sign sign
(nth 1 lhs
))))
2343 ((not (math-expr-contains (nth 2 lhs
) solve-var
))
2344 (math-try-solve-for (nth 1 lhs
)
2345 (math-mul rhs
(nth 2 lhs
))
2346 (math-solve-sign sign
(nth 2 lhs
))))
2347 ((setq t1
(math-try-solve-for (math-sub (nth 1 lhs
)
2348 (math-mul (nth 2 lhs
)
2353 (cond ((not (math-expr-contains (nth 1 lhs
) solve-var
))
2356 (math-add (math-normalize
2357 (list 'calcFunc-log rhs
(nth 1 lhs
)))
2360 (math-mul '(var pi var-pi
)
2364 (list 'calcFunc-ln
(nth 1 lhs
)))))))
2365 ((not (math-expr-contains (nth 2 lhs
) solve-var
))
2366 (cond ((and (integerp (nth 2 lhs
))
2368 (setq t1
(math-integer-log2 (nth 2 lhs
))))
2370 (if (and (eq solve-full t
)
2371 (math-known-realp (nth 1 lhs
)))
2373 (while (>= (setq t1
(1- t1
)) 0)
2374 (setq t2
(list 'calcFunc-sqrt t2
)))
2375 (setq t2
(math-solve-get-sign t2
)))
2376 (while (>= (setq t1
(1- t1
)) 0)
2377 (setq t2
(math-solve-get-sign
2379 (list 'calcFunc-sqrt t2
))))))
2382 (math-normalize t2
)))
2383 ((math-looks-negp (nth 2 lhs
))
2385 (list '^
(nth 1 lhs
) (math-neg (nth 2 lhs
)))
2387 ((and (eq solve-full t
)
2388 (Math-integerp (nth 2 lhs
))
2389 (math-known-realp (nth 1 lhs
)))
2390 (setq t1
(math-normalize
2391 (list 'calcFunc-nroot rhs
(nth 2 lhs
))))
2392 (if (math-evenp (nth 2 lhs
))
2393 (setq t1
(math-solve-get-sign t1
)))
2397 (math-oddp (nth 2 lhs
))
2398 (math-solve-sign sign
(nth 2 lhs
)))))
2399 (t (math-try-solve-for
2404 (if (Math-realp (nth 2 lhs
))
2409 (and (integerp (nth 2 lhs
))
2412 (math-div (nth 2 lhs
) 2))
2419 (and (integerp (nth 2 lhs
))
2424 (list 'calcFunc-nroot
2428 (math-oddp (nth 2 lhs
))
2429 (math-solve-sign sign
(nth 2 lhs
)))))))))
2432 (defun math-solve-prod (lsoln rsoln
)
2437 ((eq solve-full
'all
)
2438 (cons 'vec
(append (cdr lsoln
) (cdr rsoln
))))
2441 (list 'calcFunc-gt
(math-solve-get-sign 1) 0)
2446 ;;; This deals with negative, fractional, and symbolic powers of "x".
2447 (defun math-solve-poly-funny-powers (sub-rhs) ; uses "t1", "t2"
2449 (let ((pp math-poly-neg-powers
)
2452 (setq fac
(math-pow (car pp
) (or math-poly-mult-powers
1))
2453 t1
(math-mul t1 fac
)
2454 rhs
(math-mul rhs fac
)
2456 (if sub-rhs
(setq t1
(math-sub t1 rhs
)))
2457 (let ((math-poly-neg-powers nil
))
2458 (setq t2
(math-mul (or math-poly-mult-powers
1)
2459 (let ((calc-prefer-frac t
))
2460 (math-div 1 math-poly-frac-powers
)))
2461 t1
(math-is-polynomial (math-simplify (calcFunc-expand t1
)) b
50))))
2463 ;;; This converts "a x^8 + b x^5 + c x^2" to "(a (x^3)^2 + b (x^3) + c) * x^2".
2464 (defun math-solve-crunch-poly (max-degree) ; uses "t1", "t3"
2466 (while (and t1
(Math-zerop (car t1
)))
2470 (let* ((degree (1- (length t1
)))
2472 (while (and (> scale
1) (= (car t3
) 1))
2473 (and (= (% degree scale
) 0)
2479 (if (= (% n scale
) 0)
2480 (setq new-t1
(nconc new-t1
(list (car p
))))
2481 (or (Math-zerop (car p
))
2486 (setq t3
(cons scale
(cdr t3
))
2488 (setq scale
(1- scale
)))
2489 (setq t3
(list (math-mul (car t3
) t2
) (math-mul count t2
)))
2490 (<= (1- (length t1
)) max-degree
)))))
2492 (defun calcFunc-poly (expr var
&optional degree
)
2494 (or (natnump degree
) (math-reject-arg degree
'fixnatnump
))
2496 (let ((p (math-is-polynomial expr var degree
'gen
)))
2501 (math-reject-arg expr
"Expected a polynomial"))))
2503 (defun calcFunc-gpoly (expr var
&optional degree
)
2505 (or (natnump degree
) (math-reject-arg degree
'fixnatnump
))
2507 (let* ((math-poly-base-variable var
)
2508 (d (math-decompose-poly expr var degree nil
)))
2511 (math-reject-arg expr
"Expected a polynomial"))))
2513 (defun math-decompose-poly (lhs solve-var degree sub-rhs
)
2514 (let ((rhs (or sub-rhs
1))
2516 (setq t2
(math-polynomial-base
2520 (let ((math-poly-neg-powers '(1))
2521 (math-poly-mult-powers nil
)
2522 (math-poly-frac-powers 1)
2523 (math-poly-exp-base t
))
2524 (and (not (equal b lhs
))
2525 (or (not (memq (car-safe b
) '(+ -
))) sub-rhs
)
2526 (setq t3
'(1 0) t2
1
2527 t1
(math-is-polynomial lhs b
50))
2528 (if (and (equal math-poly-neg-powers
'(1))
2529 (memq math-poly-mult-powers
'(nil 1))
2530 (eq math-poly-frac-powers
1)
2532 (setq t1
(cons (math-sub (car t1
) rhs
)
2534 (math-solve-poly-funny-powers sub-rhs
))
2535 (math-solve-crunch-poly degree
)
2536 (or (math-expr-contains b solve-var
)
2537 (math-expr-contains (car t3
) solve-var
))))))))
2539 (list (math-pow t2
(car t3
))
2542 (math-pow t2
(nth 1 t3
))
2543 (math-div (math-pow t2
(nth 1 t3
)) rhs
))))))
2545 (defun math-solve-linear (var sign b a
)
2546 (math-try-solve-for var
2547 (math-div (math-neg b
) a
)
2548 (math-solve-sign sign a
)
2551 (defun math-solve-quadratic (var c b a
)
2554 (if (math-looks-evenp b
)
2555 (let ((halfb (math-div b
2)))
2559 (math-solve-get-sign
2561 (list 'calcFunc-sqrt
2562 (math-add (math-sqr halfb
)
2563 (math-mul (math-neg c
) a
))))))
2568 (math-solve-get-sign
2570 (list 'calcFunc-sqrt
2571 (math-add (math-sqr b
)
2572 (math-mul 4 (math-mul (math-neg c
) a
)))))))
2576 (defun math-solve-cubic (var d c b a
)
2577 (let* ((p (math-div b a
))
2581 (aa (math-sub q
(math-div psqr
3)))
2583 (math-div (math-sub (math-mul 2 (math-mul psqr p
))
2584 (math-mul 9 (math-mul p q
)))
2588 (math-try-solve-for (math-pow (math-add var
(math-div p
3)) 3)
2589 (math-neg bb
) nil t
)
2592 (math-mul (math-add var
(math-div p
3))
2593 (math-add (math-sqr (math-add var
(math-div p
3)))
2596 (setq m
(math-mul 2 (list 'calcFunc-sqrt
(math-div aa -
3))))
2605 (math-sub (list 'calcFunc-arccos
2606 (math-div (math-mul 3 bb
)
2610 (math-add 1 (math-solve-get-int
2613 calc-symbolic-mode
))))
2618 (defun math-solve-quartic (var d c b a aa
)
2619 (setq a
(math-div a aa
))
2620 (setq b
(math-div b aa
))
2621 (setq c
(math-div c aa
))
2622 (setq d
(math-div d aa
))
2625 (let* ((asqr (math-sqr a
))
2626 (asqr4 (math-div asqr
4))
2627 (y (let ((solve-full nil
)
2629 (math-solve-cubic solve-var
2631 (math-mul 4 (math-mul b d
))
2634 (math-sub (math-mul a c
)
2638 (rsqr (math-add (math-sub asqr4 b
) y
))
2639 (r (list 'calcFunc-sqrt rsqr
))
2640 (sign1 (math-solve-get-sign 1))
2641 (de (list 'calcFunc-sqrt
2643 (math-sub (math-mul 3 asqr4
)
2645 (if (Math-zerop rsqr
)
2649 (list 'calcFunc-sqrt
2650 (math-sub (math-sqr y
)
2656 (math-mul 4 (math-mul a b
))
2662 (math-sub (math-add (math-mul sign1
(math-div r
2))
2663 (math-solve-get-sign (math-div de
2)))
2667 (defvar math-symbolic-solve nil
)
2668 (defvar math-int-coefs nil
)
2669 (defun math-poly-all-roots (var p
&optional math-factoring
)
2671 (let* ((math-symbolic-solve calc-symbolic-mode
)
2673 (deg (1- (length p
)))
2674 (orig-p (reverse p
))
2675 (math-int-coefs nil
)
2676 (math-int-scale nil
)
2677 (math-double-roots nil
)
2678 (math-int-factors nil
)
2679 (math-int-threshold nil
)
2681 ;; If rational coefficients, look for exact rational factors.
2682 (while (and pp
(Math-ratp (car pp
)))
2685 (if (or math-factoring math-symbolic-solve
)
2687 (let ((lead (car orig-p
))
2688 (calc-prefer-frac t
)
2689 (scale (apply 'math-lcm-denoms p
)))
2690 (setq math-int-scale
(math-abs (math-mul scale lead
))
2691 math-int-threshold
(math-div '(float 5 -
2) math-int-scale
)
2692 math-int-coefs
(cdr (math-div (cons 'vec orig-p
) lead
)))))
2694 (let ((calc-prefer-frac nil
)
2695 (calc-symbolic-mode nil
)
2697 (def-p (copy-sequence orig-p
)))
2699 (if (Math-numberp (car pp
))
2702 (while (> deg
(if math-symbolic-solve
2 4))
2703 (let* ((x (math-poly-any-root def-p
'(float 0 0) nil
))
2705 (if (and (eq (car-safe x
) 'cplx
)
2706 (math-nearly-zerop (nth 2 x
) (nth 1 x
)))
2707 (setq x
(calcFunc-re x
)))
2709 (setq roots
(cons x roots
)))
2710 (or (math-numberp x
)
2711 (setq x
(math-evaluate-expr x
)))
2714 (while (setq pp
(cdr pp
))
2717 (setq b
(math-add (math-mul x b
) c
)))
2718 (setq def-p
(cdr def-p
)
2720 (setq p
(reverse def-p
))))
2722 (let ((solve-var '(var DUMMY var-DUMMY
))
2723 (math-solve-sign nil
)
2724 (math-solve-ranges nil
)
2726 (if (= (length p
) (length math-int-coefs
))
2727 (setq p
(reverse math-int-coefs
)))
2728 (setq roots
(append (cdr (apply (cond ((= deg
2)
2729 'math-solve-quadratic
)
2733 'math-solve-quartic
))
2737 (setq roots
(cons (math-div (math-neg (car p
)) (nth 1 p
))
2742 (math-poly-integer-root (car roots
))
2743 (setq roots
(cdr roots
)))
2744 (list math-int-factors
(nreverse math-int-coefs
) math-int-scale
))
2745 (let ((vec nil
) res
)
2747 (let ((root (car roots
))
2748 (solve-full (and solve-full
'all
)))
2749 (if (math-floatp root
)
2750 (setq root
(math-poly-any-root orig-p root t
)))
2751 (setq vec
(append vec
2752 (cdr (or (math-try-solve-for var root nil t
)
2753 (throw 'ouch nil
))))))
2754 (setq roots
(cdr roots
)))
2755 (setq vec
(cons 'vec
(nreverse vec
)))
2756 (if math-symbolic-solve
2757 (setq vec
(math-normalize vec
)))
2758 (if (eq solve-full t
)
2759 (list 'calcFunc-subscr
2761 (math-solve-get-int 1 (1- (length orig-p
)) 1))
2764 (defun math-lcm-denoms (&rest fracs
)
2767 (if (eq (car-safe (car fracs
)) 'frac
)
2768 (setq den
(calcFunc-lcm den
(nth 2 (car fracs
)))))
2769 (setq fracs
(cdr fracs
)))
2772 (defun math-poly-any-root (p x polish
) ; p is a reverse poly coeff list
2773 (let* ((newt (if (math-zerop x
)
2774 (math-poly-newton-root
2775 p
'(cplx (float 123 -
6) (float 1 -
4)) 4)
2776 (math-poly-newton-root p x
4)))
2777 (res (if (math-zerop (cdr newt
))
2779 (if (and (math-lessp (cdr newt
) '(float 1 -
3)) (not polish
))
2780 (setq newt
(math-poly-newton-root p
(car newt
) 30)))
2781 (if (math-zerop (cdr newt
))
2783 (math-poly-laguerre-root p x polish
)))))
2784 (and math-symbolic-solve
(math-floatp res
)
2788 (defun math-poly-newton-root (p x iters
)
2789 (let* ((calc-prefer-frac nil
)
2790 (calc-symbolic-mode nil
)
2791 (try-integer math-int-coefs
)
2793 (while (and (> (setq iters
(1- iters
)) 0)
2795 (math-working "newton" x
)
2798 (while (setq pp
(cdr pp
))
2799 (setq d
(math-add (math-mul x d
) b
)
2800 b
(math-add (math-mul x b
) (car pp
))))
2801 (not (math-zerop d
)))
2803 (setq dx
(math-div b d
)
2806 (let ((adx (math-abs-approx dx
)))
2807 (and (math-lessp adx math-int-threshold
)
2808 (let ((iroot (math-poly-integer-root x
)))
2811 (setq try-integer nil
))))))
2812 (or (not (or (eq dx
0)
2813 (math-nearly-zerop dx
(math-abs-approx x
))))
2814 (progn (setq dx
0) nil
)))))
2815 (cons x
(if (math-zerop x
)
2816 1 (math-div (math-abs-approx dx
) (math-abs-approx x
))))))
2818 (defun math-poly-integer-root (x)
2819 (and (math-lessp (calcFunc-xpon (math-abs-approx x
)) calc-internal-prec
)
2821 (let* ((calc-prefer-frac t
)
2822 (xre (calcFunc-re x
))
2823 (xim (calcFunc-im x
))
2824 (xresq (math-sqr xre
))
2825 (ximsq (math-sqr xim
)))
2826 (if (math-lessp ximsq
(calcFunc-scf xresq -
1))
2827 ;; Look for linear factor
2828 (let* ((rnd (math-div (math-round (math-mul xre math-int-scale
))
2830 (icp math-int-coefs
)
2833 (while (setq icp
(cdr icp
))
2834 (setq newcoef
(cons rem newcoef
)
2835 rem
(math-add (car icp
)
2836 (math-mul rem rnd
))))
2837 (and (math-zerop rem
)
2839 (setq math-int-coefs
(nreverse newcoef
)
2840 math-int-factors
(cons (list (math-neg rnd
))
2843 ;; Look for irreducible quadratic factor
2844 (let* ((rnd1 (math-div (math-round
2845 (math-mul xre
(math-mul -
2 math-int-scale
)))
2847 (sqscale (math-sqr math-int-scale
))
2848 (rnd0 (math-div (math-round (math-mul (math-add xresq ximsq
)
2851 (rem1 (car math-int-coefs
))
2852 (icp (cdr math-int-coefs
))
2855 (found (assoc (list rnd0 rnd1
(math-posp xim
))
2859 (setq math-double-roots
(delq found math-double-roots
)
2861 (while (setq icp
(cdr icp
))
2863 newcoef
(cons rem1 newcoef
)
2864 rem1
(math-sub rem0
(math-mul this rnd1
))
2865 rem0
(math-sub (car icp
) (math-mul this rnd0
)))))
2866 (and (math-zerop rem0
)
2868 (let ((aa (math-div rnd1 -
2)))
2869 (or found
(setq math-int-coefs
(reverse newcoef
)
2870 math-double-roots
(cons (list
2875 math-int-factors
(cons (cons rnd0 rnd1
)
2878 (let ((calc-symbolic-mode math-symbolic-solve
))
2879 (math-mul (math-sqrt (math-sub (math-sqr aa
)
2881 (if (math-negp xim
) -
1 1)))))))))))
2883 ;;; The following routine is from Numerical Recipes, section 9.5.
2884 (defun math-poly-laguerre-root (p x polish
)
2885 (let* ((calc-prefer-frac nil
)
2886 (calc-symbolic-mode nil
)
2889 (try-newt (not polish
))
2893 (and (or (< (setq iters
(1+ iters
)) 50)
2894 (math-reject-arg x
"*Laguerre's method failed to converge"))
2895 (let ((err (math-abs-approx (car p
)))
2896 (abx (math-abs-approx x
))
2900 (while (setq pp
(cdr pp
))
2901 (setq f
(math-add (math-mul x f
) d
)
2902 d
(math-add (math-mul x d
) b
)
2903 b
(math-add (math-mul x b
) (car pp
))
2904 err
(math-add (math-abs-approx b
) (math-mul abx err
))))
2905 (math-lessp (calcFunc-scf err
(- -
2 calc-internal-prec
))
2906 (math-abs-approx b
)))
2907 (or (not (math-zerop d
))
2908 (not (math-zerop f
))
2910 (setq x
(math-pow (math-neg b
) (list 'frac
1 m
)))
2912 (let* ((g (math-div d b
))
2914 (h (math-sub g2
(math-mul 2 (math-div f b
))))
2916 (math-mul (1- m
) (math-sub (math-mul m h
) g2
))))
2917 (gp (math-add g sq
))
2918 (gm (math-sub g sq
)))
2919 (if (math-lessp (calcFunc-abssqr gp
) (calcFunc-abssqr gm
))
2921 (setq dx
(math-div m gp
)
2924 (math-lessp (math-abs-approx dx
)
2925 (calcFunc-scf (math-abs-approx x
) -
3)))
2926 (let ((newt (math-poly-newton-root p x1
7)))
2929 (if (math-zerop (cdr newt
))
2930 (setq x
(car newt
) x1 x
)
2931 (if (math-lessp (cdr newt
) '(float 1 -
6))
2932 (let ((newt2 (math-poly-newton-root
2934 (if (math-zerop (cdr newt2
))
2935 (setq x
(car newt2
) x1 x
)
2936 (setq x
(car newt
))))))))
2938 (math-nearly-equal x x1
))))
2939 (let ((cdx (math-abs-approx dx
)))
2944 (math-lessp cdx dxold
)
2947 (let ((digs (calcFunc-xpon
2948 (math-div (math-abs-approx x
) cdx
))))
2950 "*Could not attain full precision")
2952 (let ((calc-internal-prec (max 3 digs
)))
2953 (setq x
(math-normalize x
))))))
2957 (math-lessp (calcFunc-scf (math-abs-approx x
)
2958 (- calc-internal-prec
))
2960 (or (and (math-floatp x
)
2961 (math-poly-integer-root x
))
2964 (defun math-solve-above-dummy (x)
2965 (and (not (Math-primp x
))
2966 (if (and (equal (nth 1 x
) '(var SOLVEDUM SOLVEDUM
))
2970 (while (and (setq x
(cdr x
))
2971 (not (setq res
(math-solve-above-dummy (car x
))))))
2974 (defun math-solve-find-root-term (x neg
) ; sets "t2", "t3"
2975 (if (math-solve-find-root-in-prod x
)
2978 (and (memq (car-safe x
) '(+ -
))
2979 (or (math-solve-find-root-term (nth 1 x
) neg
)
2980 (math-solve-find-root-term (nth 2 x
)
2981 (if (eq (car x
) '-
) (not neg
) neg
))))))
2983 (defun math-solve-find-root-in-prod (x)
2985 (math-expr-contains x solve-var
)
2986 (or (and (eq (car x
) 'calcFunc-sqrt
)
2988 (and (eq (car x
) '^
)
2989 (or (and (memq (math-quarter-integer (nth 2 x
)) '(1 2 3))
2991 (and (eq (car-safe (nth 2 x
)) 'frac
)
2992 (eq (nth 2 (nth 2 x
)) 3)
2994 (and (memq (car x
) '(* /))
2995 (or (and (not (math-expr-contains (nth 1 x
) solve-var
))
2996 (math-solve-find-root-in-prod (nth 2 x
)))
2997 (and (not (math-expr-contains (nth 2 x
) solve-var
))
2998 (math-solve-find-root-in-prod (nth 1 x
))))))))
3001 (defun math-solve-system (exprs solve-vars solve-full
)
3002 (setq exprs
(mapcar 'list
(if (Math-vectorp exprs
)
3005 solve-vars
(if (Math-vectorp solve-vars
)
3008 (or (let ((math-solve-simplifying nil
))
3009 (math-solve-system-rec exprs solve-vars nil
))
3010 (let ((math-solve-simplifying t
))
3011 (math-solve-system-rec exprs solve-vars nil
))))
3013 ;;; The following backtracking solver works by choosing a variable
3014 ;;; and equation, and trying to solve the equation for the variable.
3015 ;;; If it succeeds it calls itself recursively with that variable and
3016 ;;; equation removed from their respective lists, and with the solution
3017 ;;; added to solns as well as being substituted into all existing
3018 ;;; equations. The algorithm terminates when any solution path
3019 ;;; manages to remove all the variables from var-list.
3021 ;;; To support calcFunc-roots, entries in eqn-list and solns are
3022 ;;; actually lists of equations.
3024 (defun math-solve-system-rec (eqn-list var-list solns
)
3029 ;; Try each variable in turn.
3035 (elim (eq (car-safe vv
) 'calcFunc-elim
)))
3037 (setq vv
(nth 1 vv
)))
3039 ;; Try each equation in turn.
3048 ;; Try to solve for vv the list of equations e2.
3050 (setq res2
(or (and (eq (car e2
) eprev
)
3052 (math-solve-for (car e2
) 0 vv
3054 (setq eprev
(car e2
)
3055 res
(cons (if (eq solve-full
'all
)
3063 ;; Found a solution. Now try other variables.
3064 (setq res
(nreverse res
)
3065 res
(math-solve-system-rec
3067 'math-solve-system-subst
3069 (copy-sequence eqn-list
)))
3070 (delq (car v
) (copy-sequence var-list
))
3071 (let ((math-solve-simplifying nil
)
3077 (math-solve-system-subst
3082 (cons (cons vv
(apply 'append res
))
3090 ;; Eliminated all variables, so now put solution into the proper format.
3091 (setq solns
(sort solns
3094 (not (memq (car x
) (memq (car y
) solve-vars
)))))))
3095 (if (eq solve-full
'all
)
3100 (mapcar (function (lambda (x) (cons 'vec
(cdr x
)))) solns
)
3101 (mapcar (function (lambda (x) (cons 'vec x
))) eqn-list
)))))
3105 (mapcar (function (lambda (x) (cons 'calcFunc-eq x
))) solns
)
3106 (mapcar 'car eqn-list
)))))))
3108 (defun math-solve-system-subst (x) ; uses "res" and "v"
3112 (setq accum
(nconc accum
3115 (if math-solve-simplifying
3117 (math-expr-subst (car x
) vv r
))
3118 (math-expr-subst (car x
) vv r
))))
3125 (defun math-get-from-counter (name)
3126 (let ((ctr (assq name calc-command-flags
)))
3128 (setcdr ctr
(1+ (cdr ctr
)))
3129 (setq ctr
(cons name
1)
3130 calc-command-flags
(cons ctr calc-command-flags
)))
3133 (defun math-solve-get-sign (val)
3134 (setq val
(math-simplify val
))
3135 (if (and (eq (car-safe val
) '*)
3136 (Math-numberp (nth 1 val
)))
3137 (list '* (nth 1 val
) (math-solve-get-sign (nth 2 val
)))
3138 (and (eq (car-safe val
) 'calcFunc-sqrt
)
3139 (eq (car-safe (nth 1 val
)) '^
)
3140 (setq val
(math-normalize (list '^
3142 (math-div (nth 2 (nth 1 val
)) 2)))))
3144 (if (and (calc-var-value 'var-GenCount
)
3145 (Math-natnump var-GenCount
)
3146 (not (eq solve-full
'all
)))
3148 (math-mul (list 'calcFunc-as var-GenCount
) val
)
3149 (setq var-GenCount
(math-add var-GenCount
1))
3150 (calc-refresh-evaltos 'var-GenCount
))
3151 (let* ((var (concat "s" (int-to-string (math-get-from-counter 'solve-sign
))))
3152 (var2 (list 'var
(intern var
) (intern (concat "var-" var
)))))
3153 (if (eq solve-full
'all
)
3154 (setq math-solve-ranges
(cons (list var2
1 -
1)
3155 math-solve-ranges
)))
3156 (math-mul var2 val
)))
3157 (calc-record-why "*Choosing positive solution")
3160 (defun math-solve-get-int (val &optional range first
)
3162 (if (and (calc-var-value 'var-GenCount
)
3163 (Math-natnump var-GenCount
)
3164 (not (eq solve-full
'all
)))
3166 (math-mul val
(list 'calcFunc-an var-GenCount
))
3167 (setq var-GenCount
(math-add var-GenCount
1))
3168 (calc-refresh-evaltos 'var-GenCount
))
3169 (let* ((var (concat "n" (int-to-string
3170 (math-get-from-counter 'solve-int
))))
3171 (var2 (list 'var
(intern var
) (intern (concat "var-" var
)))))
3172 (if (and range
(eq solve-full
'all
))
3173 (setq math-solve-ranges
(cons (cons var2
3174 (cdr (calcFunc-index
3175 range
(or first
0))))
3176 math-solve-ranges
)))
3177 (math-mul val var2
)))
3178 (calc-record-why "*Choosing 0 for arbitrary integer in solution")
3181 (defun math-solve-sign (sign expr
)
3183 (let ((s1 (math-possible-signs expr
)))
3184 (cond ((memq s1
'(4 6))
3189 (defun math-looks-evenp (expr)
3190 (if (Math-integerp expr
)
3192 (if (memq (car expr
) '(* /))
3193 (math-looks-evenp (nth 1 expr
)))))
3195 (defun math-solve-for (lhs rhs solve-var solve-full
&optional sign
)
3196 (if (math-expr-contains rhs solve-var
)
3197 (math-solve-for (math-sub lhs rhs
) 0 solve-var solve-full
)
3198 (and (math-expr-contains lhs solve-var
)
3199 (math-with-extra-prec 1
3200 (let* ((math-poly-base-variable solve-var
)
3201 (res (math-try-solve-for lhs rhs sign
)))
3202 (if (and (eq solve-full
'all
)
3203 (math-known-realp solve-var
))
3204 (let ((old-len (length res
))
3209 (and (not (memq (car-safe x
)
3213 new-len
(length res
))
3214 (if (< new-len old-len
)
3215 (calc-record-why (if (= new-len
1)
3216 "*All solutions were complex"
3218 "*Omitted %d complex solutions"
3219 (- old-len new-len
)))))))
3222 (defun math-solve-eqn (expr var full
)
3223 (if (memq (car-safe expr
) '(calcFunc-neq calcFunc-lt calcFunc-gt
3224 calcFunc-leq calcFunc-geq
))
3225 (let ((res (math-solve-for (cons '-
(cdr expr
))
3227 (if (eq (car expr
) 'calcFunc-neq
) nil
1))))
3229 (if (eq math-solve-sign
1)
3230 (list (car expr
) var res
)
3231 (if (eq math-solve-sign -
1)
3232 (list (car expr
) res var
)
3233 (or (eq (car expr
) 'calcFunc-neq
)
3235 "*Can't determine direction of inequality"))
3236 (and (memq (car expr
) '(calcFunc-neq calcFunc-lt calcFunc-gt
))
3237 (list 'calcFunc-neq var res
))))))
3238 (let ((res (math-solve-for expr
0 var full
)))
3240 (list 'calcFunc-eq var res
)))))
3242 (defun math-reject-solution (expr var func
)
3243 (if (math-expr-contains expr var
)
3244 (or (equal (car calc-next-why
) '(* "Unable to find a symbolic solution"))
3245 (calc-record-why "*Unable to find a solution")))
3246 (list func expr var
))
3248 (defun calcFunc-solve (expr var
)
3249 (or (if (or (Math-vectorp expr
) (Math-vectorp var
))
3250 (math-solve-system expr var nil
)
3251 (math-solve-eqn expr var nil
))
3252 (math-reject-solution expr var
'calcFunc-solve
)))
3254 (defun calcFunc-fsolve (expr var
)
3255 (or (if (or (Math-vectorp expr
) (Math-vectorp var
))
3256 (math-solve-system expr var t
)
3257 (math-solve-eqn expr var t
))
3258 (math-reject-solution expr var
'calcFunc-fsolve
)))
3260 (defun calcFunc-roots (expr var
)
3261 (let ((math-solve-ranges nil
))
3262 (or (if (or (Math-vectorp expr
) (Math-vectorp var
))
3263 (math-solve-system expr var
'all
)
3264 (math-solve-for expr
0 var
'all
))
3265 (math-reject-solution expr var
'calcFunc-roots
))))
3267 (defun calcFunc-finv (expr var
)
3268 (let ((res (math-solve-for expr math-integ-var var nil
)))
3270 (math-normalize (math-expr-subst res math-integ-var var
))
3271 (math-reject-solution expr var
'calcFunc-finv
))))
3273 (defun calcFunc-ffinv (expr var
)
3274 (let ((res (math-solve-for expr math-integ-var var t
)))
3276 (math-normalize (math-expr-subst res math-integ-var var
))
3277 (math-reject-solution expr var
'calcFunc-finv
))))
3280 (put 'calcFunc-inv
'math-inverse
3281 (function (lambda (x) (math-div 1 x
))))
3282 (put 'calcFunc-inv
'math-inverse-sign -
1)
3284 (put 'calcFunc-sqrt
'math-inverse
3285 (function (lambda (x) (math-sqr x
))))
3287 (put 'calcFunc-conj
'math-inverse
3288 (function (lambda (x) (list 'calcFunc-conj x
))))
3290 (put 'calcFunc-abs
'math-inverse
3291 (function (lambda (x) (math-solve-get-sign x
))))
3293 (put 'calcFunc-deg
'math-inverse
3294 (function (lambda (x) (list 'calcFunc-rad x
))))
3295 (put 'calcFunc-deg
'math-inverse-sign
1)
3297 (put 'calcFunc-rad
'math-inverse
3298 (function (lambda (x) (list 'calcFunc-deg x
))))
3299 (put 'calcFunc-rad
'math-inverse-sign
1)
3301 (put 'calcFunc-ln
'math-inverse
3302 (function (lambda (x) (list 'calcFunc-exp x
))))
3303 (put 'calcFunc-ln
'math-inverse-sign
1)
3305 (put 'calcFunc-log10
'math-inverse
3306 (function (lambda (x) (list 'calcFunc-exp10 x
))))
3307 (put 'calcFunc-log10
'math-inverse-sign
1)
3309 (put 'calcFunc-lnp1
'math-inverse
3310 (function (lambda (x) (list 'calcFunc-expm1 x
))))
3311 (put 'calcFunc-lnp1
'math-inverse-sign
1)
3313 (put 'calcFunc-exp
'math-inverse
3314 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-ln x
))
3316 (math-mul '(var pi var-pi
)
3318 '(var i var-i
))))))))
3319 (put 'calcFunc-exp
'math-inverse-sign
1)
3321 (put 'calcFunc-expm1
'math-inverse
3322 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-lnp1 x
))
3324 (math-mul '(var pi var-pi
)
3326 '(var i var-i
))))))))
3327 (put 'calcFunc-expm1
'math-inverse-sign
1)
3329 (put 'calcFunc-sin
'math-inverse
3330 (function (lambda (x) (let ((n (math-solve-get-int 1)))
3331 (math-add (math-mul (math-normalize
3332 (list 'calcFunc-arcsin x
))
3334 (math-mul (math-half-circle t
)
3337 (put 'calcFunc-cos
'math-inverse
3338 (function (lambda (x) (math-add (math-solve-get-sign
3340 (list 'calcFunc-arccos x
)))
3342 (math-full-circle t
))))))
3344 (put 'calcFunc-tan
'math-inverse
3345 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-arctan x
))
3347 (math-half-circle t
))))))
3349 (put 'calcFunc-arcsin
'math-inverse
3350 (function (lambda (x) (math-normalize (list 'calcFunc-sin x
)))))
3352 (put 'calcFunc-arccos
'math-inverse
3353 (function (lambda (x) (math-normalize (list 'calcFunc-cos x
)))))
3355 (put 'calcFunc-arctan
'math-inverse
3356 (function (lambda (x) (math-normalize (list 'calcFunc-tan x
)))))
3358 (put 'calcFunc-sinh
'math-inverse
3359 (function (lambda (x) (let ((n (math-solve-get-int 1)))
3360 (math-add (math-mul (math-normalize
3361 (list 'calcFunc-arcsinh x
))
3363 (math-mul (math-half-circle t
)
3367 (put 'calcFunc-sinh
'math-inverse-sign
1)
3369 (put 'calcFunc-cosh
'math-inverse
3370 (function (lambda (x) (math-add (math-solve-get-sign
3372 (list 'calcFunc-arccosh x
)))
3373 (math-mul (math-full-circle t
)
3375 '(var i var-i
)))))))
3377 (put 'calcFunc-tanh
'math-inverse
3378 (function (lambda (x) (math-add (math-normalize
3379 (list 'calcFunc-arctanh x
))
3380 (math-mul (math-half-circle t
)
3382 '(var i var-i
)))))))
3383 (put 'calcFunc-tanh
'math-inverse-sign
1)
3385 (put 'calcFunc-arcsinh
'math-inverse
3386 (function (lambda (x) (math-normalize (list 'calcFunc-sinh x
)))))
3387 (put 'calcFunc-arcsinh
'math-inverse-sign
1)
3389 (put 'calcFunc-arccosh
'math-inverse
3390 (function (lambda (x) (math-normalize (list 'calcFunc-cosh x
)))))
3392 (put 'calcFunc-arctanh
'math-inverse
3393 (function (lambda (x) (math-normalize (list 'calcFunc-tanh x
)))))
3394 (put 'calcFunc-arctanh
'math-inverse-sign
1)
3398 (defun calcFunc-taylor (expr var num
)
3399 (let ((x0 0) (v var
))
3400 (if (memq (car-safe var
) '(+ - calcFunc-eq
))
3401 (setq x0
(if (eq (car var
) '+) (math-neg (nth 2 var
)) (nth 2 var
))
3403 (or (and (eq (car-safe v
) 'var
)
3404 (math-expr-contains expr v
)
3406 (let ((accum (math-expr-subst expr v x0
))
3407 (var2 (if (eq (car var
) 'calcFunc-eq
)
3413 (while (and (<= (setq n
(1+ n
)) num
)
3414 (setq fprime
(calcFunc-deriv fprime v nil t
)))
3415 (setq fprime
(math-simplify fprime
)
3416 nfac
(math-mul nfac n
)
3417 accum
(math-add accum
3418 (math-div (math-mul (math-pow var2 n
)
3423 (math-normalize accum
))))
3424 (list 'calcFunc-taylor expr var num
))))
3426 ;;; arch-tag: f2932ec8-dd63-418b-a542-11a644b9d4c4
3427 ;;; calcalg2.el ends here
