1 ;;; calcalg2.el --- more algebraic functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
4 ;; 2005, 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
6 ;; Author: David Gillespie <daveg@synaptics.com>
7 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
9 ;; This file is part of GNU Emacs.
11 ;; GNU Emacs is free software: you can redistribute it and/or modify
12 ;; it under the terms of the GNU General Public License as published by
13 ;; the Free Software Foundation, either version 3 of the License, or
14 ;; (at your option) any later version.
16 ;; GNU Emacs is distributed in the hope that it will be useful,
17 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
18 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 ;; GNU General Public License for more details.
21 ;; You should have received a copy of the GNU General Public License
22 ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
28 ;; This file is autoloaded from calc-ext.el.
33 (defun calc-derivative (var num
)
34 (interactive "sDifferentiate with respect to: \np")
37 (error "Order of derivative must be positive"))
38 (let ((func (if (calc-is-hyperbolic) 'calcFunc-tderiv
'calcFunc-deriv
))
40 (if (or (equal var
"") (equal var
"$"))
44 (setq var
(math-read-expr var
))
45 (when (eq (car-safe var
) 'error
)
46 (error "Bad format in expression: %s" (nth 1 var
)))
49 (while (>= (setq num
(1- num
)) 0)
50 (setq expr
(list func expr var
)))
51 (calc-enter-result n
"derv" expr
))))
53 (defun calc-integral (var &optional arg
)
54 (interactive "sIntegration variable: \nP")
56 (calc-tabular-command 'calcFunc-integ
"Integration" "intg" nil var nil nil
)
58 (if (or (equal var
"") (equal var
"$"))
59 (calc-enter-result 2 "intg" (list 'calcFunc-integ
62 (let ((var (math-read-expr var
)))
63 (if (eq (car-safe var
) 'error
)
64 (error "Bad format in expression: %s" (nth 1 var
)))
65 (calc-enter-result 1 "intg" (list 'calcFunc-integ
69 (defun calc-num-integral (&optional varname lowname highname
)
70 (interactive "sIntegration variable: ")
71 (calc-tabular-command 'calcFunc-ninteg
"Integration" "nint"
72 nil varname lowname highname
))
74 (defun calc-summation (arg &optional varname lowname highname
)
75 (interactive "P\nsSummation variable: ")
76 (calc-tabular-command 'calcFunc-sum
"Summation" "sum"
77 arg varname lowname highname
))
79 (defun calc-alt-summation (arg &optional varname lowname highname
)
80 (interactive "P\nsSummation variable: ")
81 (calc-tabular-command 'calcFunc-asum
"Summation" "asum"
82 arg varname lowname highname
))
84 (defun calc-product (arg &optional varname lowname highname
)
85 (interactive "P\nsIndex variable: ")
86 (calc-tabular-command 'calcFunc-prod
"Index" "prod"
87 arg varname lowname highname
))
89 (defun calc-tabulate (arg &optional varname lowname highname
)
90 (interactive "P\nsIndex variable: ")
91 (calc-tabular-command 'calcFunc-table
"Index" "tabl"
92 arg varname lowname highname
))
94 (defun calc-tabular-command (func prompt prefix arg varname lowname highname
)
96 (let (var (low nil
) (high nil
) (step nil
) stepname stepnum
(num 1) expr
)
100 (if (or (equal varname
"") (equal varname
"$") (null varname
))
101 (setq high
(calc-top-n (+ stepnum
1))
102 low
(calc-top-n (+ stepnum
2))
103 var
(calc-top-n (+ stepnum
3))
105 (setq var
(if (stringp varname
) (math-read-expr varname
) varname
))
106 (if (eq (car-safe var
) 'error
)
107 (error "Bad format in expression: %s" (nth 1 var
)))
109 (setq lowname
(read-string (concat prompt
" variable: " varname
111 (if (or (equal lowname
"") (equal lowname
"$"))
112 (setq high
(calc-top-n (+ stepnum
1))
113 low
(calc-top-n (+ stepnum
2))
115 (setq low
(if (stringp lowname
) (math-read-expr lowname
) lowname
))
116 (if (eq (car-safe low
) 'error
)
117 (error "Bad format in expression: %s" (nth 1 low
)))
119 (setq highname
(read-string (concat prompt
" variable: " varname
122 (if (or (equal highname
"") (equal highname
"$"))
123 (setq high
(calc-top-n (+ stepnum
1))
125 (setq high
(if (stringp highname
) (math-read-expr highname
)
127 (if (eq (car-safe high
) 'error
)
128 (error "Bad format in expression: %s" (nth 1 high
)))
131 (setq stepname
(read-string (concat prompt
" variable: "
136 (if (or (equal stepname
"") (equal stepname
"$"))
137 (setq step
(calc-top-n 1)
139 (setq step
(math-read-expr stepname
))
140 (if (eq (car-safe step
) 'error
)
141 (error "Bad format in expression: %s"
145 (setq step
(calc-top-n 1))
147 (setq step
(prefix-numeric-value arg
)))))
148 (setq expr
(calc-top-n num
))
149 (calc-enter-result num prefix
(append (list func expr var low high
)
150 (and step
(list step
)))))))
152 (defun calc-solve-for (var)
153 (interactive "sVariable(s) to solve for: ")
155 (let ((func (if (calc-is-inverse)
156 (if (calc-is-hyperbolic) 'calcFunc-ffinv
'calcFunc-finv
)
157 (if (calc-is-hyperbolic) 'calcFunc-fsolve
'calcFunc-solve
))))
158 (if (or (equal var
"") (equal var
"$"))
159 (calc-enter-result 2 "solv" (list func
162 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var
)
163 (not (string-match "\\[" var
)))
164 (math-read-expr (concat "[" var
"]"))
165 (math-read-expr var
))))
166 (if (eq (car-safe var
) 'error
)
167 (error "Bad format in expression: %s" (nth 1 var
)))
168 (calc-enter-result 1 "solv" (list func
172 (defun calc-poly-roots (var)
173 (interactive "sVariable to solve for: ")
175 (if (or (equal var
"") (equal var
"$"))
176 (calc-enter-result 2 "prts" (list 'calcFunc-roots
179 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var
)
180 (not (string-match "\\[" var
)))
181 (math-read-expr (concat "[" var
"]"))
182 (math-read-expr var
))))
183 (if (eq (car-safe var
) 'error
)
184 (error "Bad format in expression: %s" (nth 1 var
)))
185 (calc-enter-result 1 "prts" (list 'calcFunc-roots
189 (defun calc-taylor (var nterms
)
190 (interactive "sTaylor expansion variable: \nNNumber of terms: ")
192 (let ((var (math-read-expr var
)))
193 (if (eq (car-safe var
) 'error
)
194 (error "Bad format in expression: %s" (nth 1 var
)))
195 (calc-enter-result 1 "tylr" (list 'calcFunc-taylor
198 (prefix-numeric-value nterms
))))))
201 ;; The following are global variables used by math-derivative and some
203 (defvar math-deriv-var
)
204 (defvar math-deriv-total
)
205 (defvar math-deriv-symb
)
206 (defvar math-decls-cache
)
207 (defvar math-decls-all
)
209 (defun math-derivative (expr)
210 (cond ((equal expr math-deriv-var
)
212 ((or (Math-scalarp expr
)
213 (eq (car expr
) 'sdev
)
214 (and (eq (car expr
) 'var
)
215 (or (not math-deriv-total
)
216 (math-const-var expr
)
218 (math-setup-declarations)
219 (memq 'const
(nth 1 (or (assq (nth 2 expr
)
221 math-decls-all
)))))))
224 (math-add (math-derivative (nth 1 expr
))
225 (math-derivative (nth 2 expr
))))
227 (math-sub (math-derivative (nth 1 expr
))
228 (math-derivative (nth 2 expr
))))
229 ((memq (car expr
) '(calcFunc-eq calcFunc-neq calcFunc-lt
230 calcFunc-gt calcFunc-leq calcFunc-geq
))
232 (math-derivative (nth 1 expr
))
233 (math-derivative (nth 2 expr
))))
234 ((eq (car expr
) 'neg
)
235 (math-neg (math-derivative (nth 1 expr
))))
237 (math-add (math-mul (nth 2 expr
)
238 (math-derivative (nth 1 expr
)))
239 (math-mul (nth 1 expr
)
240 (math-derivative (nth 2 expr
)))))
242 (math-sub (math-div (math-derivative (nth 1 expr
))
244 (math-div (math-mul (nth 1 expr
)
245 (math-derivative (nth 2 expr
)))
246 (math-sqr (nth 2 expr
)))))
248 (let ((du (math-derivative (nth 1 expr
)))
249 (dv (math-derivative (nth 2 expr
))))
251 (setq du
(math-mul (nth 2 expr
)
252 (math-mul (math-normalize
255 (math-add (nth 2 expr
) -
1)))
258 (setq dv
(math-mul (math-normalize
259 (list 'calcFunc-ln
(nth 1 expr
)))
260 (math-mul expr dv
))))
263 (math-derivative (nth 1 expr
))) ; a reasonable definition
264 ((eq (car expr
) 'vec
)
265 (math-map-vec 'math-derivative expr
))
266 ((and (memq (car expr
) '(calcFunc-conj calcFunc-re calcFunc-im
))
268 (list (car expr
) (math-derivative (nth 1 expr
))))
269 ((and (memq (car expr
) '(calcFunc-subscr calcFunc-mrow calcFunc-mcol
))
271 (let ((d (math-derivative (nth 1 expr
))))
273 0 ; assume x and x_1 are independent vars
274 (list (car expr
) d
(nth 2 expr
)))))
275 (t (or (and (symbolp (car expr
))
276 (if (= (length expr
) 2)
277 (let ((handler (get (car expr
) 'math-derivative
)))
279 (let ((deriv (math-derivative (nth 1 expr
))))
280 (if (Math-zerop deriv
)
282 (math-mul (funcall handler
(nth 1 expr
))
284 (let ((handler (get (car expr
) 'math-derivative-n
)))
286 (funcall handler expr
)))))
287 (and (not (eq math-deriv-symb
'pre-expand
))
288 (let ((exp (math-expand-formula expr
)))
290 (or (let ((math-deriv-symb 'pre-expand
))
291 (catch 'math-deriv
(math-derivative expr
)))
292 (math-derivative exp
)))))
293 (if (or (Math-objvecp expr
)
295 (not (symbolp (car expr
))))
297 (throw 'math-deriv nil
)
298 (list (if math-deriv-total
'calcFunc-tderiv
'calcFunc-deriv
)
305 (while (setq arg
(cdr arg
))
306 (or (Math-zerop (setq derv
(math-derivative (car arg
))))
307 (let ((func (intern (concat (symbol-name (car expr
))
312 (prop (cond ((= (length expr
) 2)
321 'math-derivative-5
))))
327 (let ((handler (get func prop
)))
328 (or (and prop handler
329 (apply handler
(cdr expr
)))
330 (if (and math-deriv-symb
333 (throw 'math-deriv nil
)
334 (cons func
(cdr expr
))))))))))
338 (defun calcFunc-deriv (expr math-deriv-var
&optional deriv-value math-deriv-symb
)
339 (let* ((math-deriv-total nil
)
340 (res (catch 'math-deriv
(math-derivative expr
))))
341 (or (eq (car-safe res
) 'calcFunc-deriv
)
343 (setq res
(math-normalize res
)))
346 (math-expr-subst res math-deriv-var deriv-value
)
349 (defun calcFunc-tderiv (expr math-deriv-var
&optional deriv-value math-deriv-symb
)
350 (math-setup-declarations)
351 (let* ((math-deriv-total t
)
352 (res (catch 'math-deriv
(math-derivative expr
))))
353 (or (eq (car-safe res
) 'calcFunc-tderiv
)
355 (setq res
(math-normalize res
)))
358 (math-expr-subst res math-deriv-var deriv-value
)
361 (put 'calcFunc-inv
\' 'math-derivative-1
362 (function (lambda (u) (math-neg (math-div 1 (math-sqr u
))))))
364 (put 'calcFunc-sqrt
\' 'math-derivative-1
365 (function (lambda (u) (math-div 1 (math-mul 2 (list 'calcFunc-sqrt u
))))))
367 (put 'calcFunc-deg
\' 'math-derivative-1
368 (function (lambda (u) (math-div-float '(float 18 1) (math-pi)))))
370 (put 'calcFunc-rad
\' 'math-derivative-1
371 (function (lambda (u) (math-pi-over-180))))
373 (put 'calcFunc-ln
\' 'math-derivative-1
374 (function (lambda (u) (math-div 1 u
))))
376 (put 'calcFunc-log10
\' 'math-derivative-1
377 (function (lambda (u)
378 (math-div (math-div 1 (math-normalize '(calcFunc-ln 10)))
381 (put 'calcFunc-lnp1
\' 'math-derivative-1
382 (function (lambda (u) (math-div 1 (math-add u
1)))))
384 (put 'calcFunc-log
\' 'math-derivative-2
385 (function (lambda (x b
)
386 (and (not (Math-zerop b
))
387 (let ((lnv (math-normalize
388 (list 'calcFunc-ln b
))))
389 (math-div 1 (math-mul lnv x
)))))))
391 (put 'calcFunc-log
\'2 'math-derivative-2
392 (function (lambda (x b
)
393 (let ((lnv (list 'calcFunc-ln b
)))
394 (math-neg (math-div (list 'calcFunc-log x b
)
395 (math-mul lnv b
)))))))
397 (put 'calcFunc-exp
\' 'math-derivative-1
398 (function (lambda (u) (math-normalize (list 'calcFunc-exp u
)))))
400 (put 'calcFunc-expm1
\' 'math-derivative-1
401 (function (lambda (u) (math-normalize (list 'calcFunc-expm1 u
)))))
403 (put 'calcFunc-sin
\' 'math-derivative-1
404 (function (lambda (u) (math-to-radians-2 (math-normalize
405 (list 'calcFunc-cos u
))))))
407 (put 'calcFunc-cos
\' 'math-derivative-1
408 (function (lambda (u) (math-neg (math-to-radians-2
410 (list 'calcFunc-sin u
)))))))
412 (put 'calcFunc-tan
\' 'math-derivative-1
413 (function (lambda (u) (math-to-radians-2
416 (list 'calcFunc-sec u
)))))))
418 (put 'calcFunc-sec
\' 'math-derivative-1
419 (function (lambda (u) (math-to-radians-2
422 (list 'calcFunc-sec u
))
424 (list 'calcFunc-tan u
)))))))
426 (put 'calcFunc-csc
\' 'math-derivative-1
427 (function (lambda (u) (math-neg
431 (list 'calcFunc-csc u
))
433 (list 'calcFunc-cot u
))))))))
435 (put 'calcFunc-cot
\' 'math-derivative-1
436 (function (lambda (u) (math-neg
440 (list 'calcFunc-csc u
))))))))
442 (put 'calcFunc-arcsin
\' 'math-derivative-1
443 (function (lambda (u)
445 (math-div 1 (math-normalize
447 (math-sub 1 (math-sqr u
)))))))))
449 (put 'calcFunc-arccos
\' 'math-derivative-1
450 (function (lambda (u)
452 (math-div -
1 (math-normalize
454 (math-sub 1 (math-sqr u
)))))))))
456 (put 'calcFunc-arctan
\' 'math-derivative-1
457 (function (lambda (u) (math-from-radians-2
458 (math-div 1 (math-add 1 (math-sqr u
)))))))
460 (put 'calcFunc-sinh
\' 'math-derivative-1
461 (function (lambda (u) (math-normalize (list 'calcFunc-cosh u
)))))
463 (put 'calcFunc-cosh
\' 'math-derivative-1
464 (function (lambda (u) (math-normalize (list 'calcFunc-sinh u
)))))
466 (put 'calcFunc-tanh
\' 'math-derivative-1
467 (function (lambda (u) (math-sqr
469 (list 'calcFunc-sech u
))))))
471 (put 'calcFunc-sech
\' 'math-derivative-1
472 (function (lambda (u) (math-neg
474 (math-normalize (list 'calcFunc-sech u
))
475 (math-normalize (list 'calcFunc-tanh u
)))))))
477 (put 'calcFunc-csch
\' 'math-derivative-1
478 (function (lambda (u) (math-neg
480 (math-normalize (list 'calcFunc-csch u
))
481 (math-normalize (list 'calcFunc-coth u
)))))))
483 (put 'calcFunc-coth
\' 'math-derivative-1
484 (function (lambda (u) (math-neg
487 (list 'calcFunc-csch u
)))))))
489 (put 'calcFunc-arcsinh
\' 'math-derivative-1
490 (function (lambda (u)
491 (math-div 1 (math-normalize
493 (math-add (math-sqr u
) 1)))))))
495 (put 'calcFunc-arccosh
\' 'math-derivative-1
496 (function (lambda (u)
497 (math-div 1 (math-normalize
499 (math-add (math-sqr u
) -
1)))))))
501 (put 'calcFunc-arctanh
\' 'math-derivative-1
502 (function (lambda (u) (math-div 1 (math-sub 1 (math-sqr u
))))))
504 (put 'calcFunc-bern
\'2 'math-derivative-2
505 (function (lambda (n x
)
506 (math-mul n
(list 'calcFunc-bern
(math-add n -
1) x
)))))
508 (put 'calcFunc-euler
\'2 'math-derivative-2
509 (function (lambda (n x
)
510 (math-mul n
(list 'calcFunc-euler
(math-add n -
1) x
)))))
512 (put 'calcFunc-gammag
\'2 'math-derivative-2
513 (function (lambda (a x
) (math-deriv-gamma a x
1))))
515 (put 'calcFunc-gammaG
\'2 'math-derivative-2
516 (function (lambda (a x
) (math-deriv-gamma a x -
1))))
518 (put 'calcFunc-gammaP
\'2 'math-derivative-2
519 (function (lambda (a x
) (math-deriv-gamma a x
522 (list 'calcFunc-gamma
525 (put 'calcFunc-gammaQ
\'2 'math-derivative-2
526 (function (lambda (a x
) (math-deriv-gamma a x
529 (list 'calcFunc-gamma
532 (defun math-deriv-gamma (a x scale
)
534 (math-mul (math-pow x
(math-add a -
1))
535 (list 'calcFunc-exp
(math-neg x
)))))
537 (put 'calcFunc-betaB
\' 'math-derivative-3
538 (function (lambda (x a b
) (math-deriv-beta x a b
1))))
540 (put 'calcFunc-betaI
\' 'math-derivative-3
541 (function (lambda (x a b
) (math-deriv-beta x a b
543 1 (list 'calcFunc-beta
546 (defun math-deriv-beta (x a b scale
)
547 (math-mul (math-mul (math-pow x
(math-add a -
1))
548 (math-pow (math-sub 1 x
) (math-add b -
1)))
551 (put 'calcFunc-erf
\' 'math-derivative-1
552 (function (lambda (x) (math-div 2
553 (math-mul (list 'calcFunc-exp
555 (if calc-symbolic-mode
560 (put 'calcFunc-erfc
\' 'math-derivative-1
561 (function (lambda (x) (math-div -
2
562 (math-mul (list 'calcFunc-exp
564 (if calc-symbolic-mode
569 (put 'calcFunc-besJ
\'2 'math-derivative-2
570 (function (lambda (v z
) (math-div (math-sub (list 'calcFunc-besJ
578 (put 'calcFunc-besY
\'2 'math-derivative-2
579 (function (lambda (v z
) (math-div (math-sub (list 'calcFunc-besY
587 (put 'calcFunc-sum
'math-derivative-n
590 (if (math-expr-contains (cons 'vec
(cdr (cdr expr
))) math-deriv-var
)
591 (throw 'math-deriv nil
)
593 (cons (math-derivative (nth 1 expr
))
594 (cdr (cdr expr
))))))))
596 (put 'calcFunc-prod
'math-derivative-n
599 (if (math-expr-contains (cons 'vec
(cdr (cdr expr
))) math-deriv-var
)
600 (throw 'math-deriv nil
)
603 (cons (math-div (math-derivative (nth 1 expr
))
605 (cdr (cdr expr
)))))))))
607 (put 'calcFunc-integ
'math-derivative-n
610 (if (= (length expr
) 3)
611 (if (equal (nth 2 expr
) math-deriv-var
)
614 (list 'calcFunc-integ
615 (math-derivative (nth 1 expr
))
617 (if (= (length expr
) 5)
618 (let ((lower (math-expr-subst (nth 1 expr
) (nth 2 expr
)
620 (upper (math-expr-subst (nth 1 expr
) (nth 2 expr
)
622 (math-add (math-sub (math-mul upper
623 (math-derivative (nth 4 expr
)))
625 (math-derivative (nth 3 expr
))))
626 (if (equal (nth 2 expr
) math-deriv-var
)
629 (list 'calcFunc-integ
630 (math-derivative (nth 1 expr
)) (nth 2 expr
)
631 (nth 3 expr
) (nth 4 expr
)))))))))))
633 (put 'calcFunc-if
'math-derivative-n
636 (and (= (length expr
) 4)
637 (list 'calcFunc-if
(nth 1 expr
)
638 (math-derivative (nth 2 expr
))
639 (math-derivative (nth 3 expr
)))))))
641 (put 'calcFunc-subscr
'math-derivative-n
644 (and (= (length expr
) 3)
645 (list 'calcFunc-subscr
(nth 1 expr
)
646 (math-derivative (nth 2 expr
)))))))
649 (defvar math-integ-var
'(var X ---
))
650 (defvar math-integ-var-2
'(var Y ---
))
651 (defvar math-integ-vars
(list 'f math-integ-var math-integ-var-2
))
652 (defvar math-integ-var-list
(list math-integ-var
))
653 (defvar math-integ-var-list-list
(list math-integ-var-list
))
655 ;; math-integ-depth is a local variable for math-try-integral, but is used
656 ;; by math-integral and math-tracing-integral
657 ;; which are called (directly or indirectly) by math-try-integral.
658 (defvar math-integ-depth
)
659 ;; math-integ-level is a local variable for math-try-integral, but is used
660 ;; by math-integral, math-do-integral, math-tracing-integral,
661 ;; math-sub-integration, math-integrate-by-parts and
662 ;; math-integrate-by-substitution, which are called (directly or
663 ;; indirectly) by math-try-integral.
664 (defvar math-integ-level
)
665 ;; math-integral-limit is a local variable for calcFunc-integ, but is
666 ;; used by math-tracing-integral, math-sub-integration and
667 ;; math-try-integration.
668 (defvar math-integral-limit
)
670 (defmacro math-tracing-integral
(&rest parts
)
673 (list 'with-current-buffer
675 '(goto-char (point-max))
678 '(insert (make-string (- math-integral-limit
679 math-integ-level
) 32)
680 (format "%2d " math-integ-depth
)
681 (make-string math-integ-level
32)))
682 ;;(list 'condition-case 'err
684 ;; '(error (insert (prin1-to-string err))))
687 ;;; The following wrapper caches results and avoids infinite recursion.
688 ;;; Each cache entry is: ( A B ) Integral of A is B;
689 ;;; ( A N ) Integral of A failed at level N;
690 ;;; ( A busy ) Currently working on integral of A;
691 ;;; ( A parts ) Currently working, integ-by-parts;
692 ;;; ( A parts2 ) Currently working, integ-by-parts;
693 ;;; ( A cancelled ) Ignore this cache entry;
694 ;;; ( A [B] ) Same result as for math-cur-record = B.
696 ;; math-cur-record is a local variable for math-try-integral, but is used
697 ;; by math-integral, math-replace-integral-parts and math-integrate-by-parts
698 ;; which are called (directly or indirectly) by math-try-integral, as well as
699 ;; by calc-dump-integral-cache
700 (defvar math-cur-record
)
701 ;; math-enable-subst and math-any-substs are local variables for
702 ;; calcFunc-integ, but are used by math-integral and math-try-integral.
703 (defvar math-enable-subst
)
704 (defvar math-any-substs
)
706 ;; math-integ-msg is a local variable for math-try-integral, but is
707 ;; used (both locally and non-locally) by math-integral.
708 (defvar math-integ-msg
)
710 (defvar math-integral-cache nil
)
711 (defvar math-integral-cache-state nil
)
713 (defun math-integral (expr &optional simplify same-as-above
)
714 (let* ((simp math-cur-record
)
715 (math-cur-record (assoc expr math-integral-cache
))
716 (math-integ-depth (1+ math-integ-depth
))
718 (math-tracing-integral "Integrating "
719 (math-format-value expr
1000)
723 (math-tracing-integral "Found "
724 (math-format-value (nth 1 math-cur-record
) 1000))
725 (and (consp (nth 1 math-cur-record
))
726 (math-replace-integral-parts math-cur-record
))
727 (math-tracing-integral " => "
728 (math-format-value (nth 1 math-cur-record
) 1000)
730 (or (and math-cur-record
731 (not (eq (nth 1 math-cur-record
) 'cancelled
))
732 (or (not (integerp (nth 1 math-cur-record
)))
733 (>= (nth 1 math-cur-record
) math-integ-level
)))
734 (and (math-integral-contains-parts expr
)
740 (let (math-integ-msg)
741 (if (eq calc-display-working-message
'lots
)
743 (calc-set-command-flag 'clear-message
)
744 (setq math-integ-msg
(format
745 "Working... Integrating %s"
746 (math-format-flat-expr expr
0)))
747 (message "%s" math-integ-msg
)))
749 (setcar (cdr math-cur-record
)
750 (if same-as-above
(vector simp
) 'busy
))
751 (setq math-cur-record
752 (list expr
(if same-as-above
(vector simp
) 'busy
))
753 math-integral-cache
(cons math-cur-record
754 math-integral-cache
)))
755 (if (eq simplify
'yes
)
757 (math-tracing-integral "Simplifying...")
758 (setq simp
(math-simplify expr
))
759 (setq val
(if (equal simp expr
)
761 (math-tracing-integral " no change\n")
762 (math-do-integral expr
))
763 (math-tracing-integral " simplified\n")
764 (math-integral simp
'no t
))))
765 (or (setq val
(math-do-integral expr
))
767 (let ((simp (math-simplify expr
)))
768 (or (equal simp expr
)
770 (math-tracing-integral "Trying again after "
771 "simplification...\n")
772 (setq val
(math-integral simp
'no t
))))))))
773 (if (eq calc-display-working-message
'lots
)
774 (message "%s" math-integ-msg
)))
775 (setcar (cdr math-cur-record
) (or val
776 (if (or math-enable-subst
777 (not math-any-substs
))
780 (setq val math-cur-record
)
781 (while (vectorp (nth 1 val
))
782 (setq val
(aref (nth 1 val
) 0)))
783 (setq val
(if (memq (nth 1 val
) '(parts parts2
))
785 (setcar (cdr val
) 'parts2
)
786 (list 'var
'PARTS val
))
787 (and (consp (nth 1 val
))
789 (math-tracing-integral "Integral of "
790 (math-format-value expr
1000)
792 (math-format-value val
1000)
796 (defun math-integral-contains-parts (expr)
797 (if (Math-primp expr
)
798 (and (eq (car-safe expr
) 'var
)
799 (eq (nth 1 expr
) 'PARTS
)
800 (listp (nth 2 expr
)))
801 (while (and (setq expr
(cdr expr
))
802 (not (math-integral-contains-parts (car expr
)))))
805 (defun math-replace-integral-parts (expr)
806 (or (Math-primp expr
)
807 (while (setq expr
(cdr expr
))
808 (and (consp (car expr
))
809 (if (eq (car (car expr
)) 'var
)
810 (and (eq (nth 1 (car expr
)) 'PARTS
)
811 (consp (nth 2 (car expr
)))
812 (if (listp (nth 1 (nth 2 (car expr
))))
814 (setcar expr
(nth 1 (nth 2 (car expr
))))
815 (math-replace-integral-parts (cons 'foo expr
)))
816 (setcar (cdr math-cur-record
) 'cancelled
)))
817 (math-replace-integral-parts (car expr
)))))))
819 (defvar math-linear-subst-tried t
820 "Non-nil means that a linear substitution has been tried.")
822 ;; The variable math-has-rules is a local variable for math-try-integral,
823 ;; but is used by math-do-integral, which is called (non-directly) by
824 ;; math-try-integral.
825 (defvar math-has-rules
)
827 ;; math-old-integ is a local variable for math-do-integral, but is
828 ;; used by math-sub-integration.
829 (defvar math-old-integ
)
831 ;; The variables math-t1, math-t2 and math-t3 are local to
832 ;; math-do-integral, math-try-solve-for and math-decompose-poly, but
833 ;; are used by functions they call (directly or indirectly);
834 ;; math-do-integral calls math-do-integral-methods;
835 ;; math-try-solve-for calls math-try-solve-prod,
836 ;; math-solve-find-root-term and math-solve-find-root-in-prod;
837 ;; math-decompose-poly calls math-solve-poly-funny-powers and
838 ;; math-solve-crunch-poly.
843 (defun math-do-integral (expr)
844 (let ((math-linear-subst-tried nil
)
846 (or (cond ((not (math-expr-contains expr math-integ-var
))
847 (math-mul expr math-integ-var
))
848 ((equal expr math-integ-var
)
849 (math-div (math-sqr expr
) 2))
851 (and (setq math-t1
(math-integral (nth 1 expr
)))
852 (setq math-t2
(math-integral (nth 2 expr
)))
853 (math-add math-t1 math-t2
)))
855 (and (setq math-t1
(math-integral (nth 1 expr
)))
856 (setq math-t2
(math-integral (nth 2 expr
)))
857 (math-sub math-t1 math-t2
)))
858 ((eq (car expr
) 'neg
)
859 (and (setq math-t1
(math-integral (nth 1 expr
)))
862 (cond ((not (math-expr-contains (nth 1 expr
) math-integ-var
))
863 (and (setq math-t1
(math-integral (nth 2 expr
)))
864 (math-mul (nth 1 expr
) math-t1
)))
865 ((not (math-expr-contains (nth 2 expr
) math-integ-var
))
866 (and (setq math-t1
(math-integral (nth 1 expr
)))
867 (math-mul math-t1
(nth 2 expr
))))
868 ((memq (car-safe (nth 1 expr
)) '(+ -
))
869 (math-integral (list (car (nth 1 expr
))
870 (math-mul (nth 1 (nth 1 expr
))
872 (math-mul (nth 2 (nth 1 expr
))
875 ((memq (car-safe (nth 2 expr
)) '(+ -
))
876 (math-integral (list (car (nth 2 expr
))
877 (math-mul (nth 1 (nth 2 expr
))
879 (math-mul (nth 2 (nth 2 expr
))
883 (cond ((and (not (math-expr-contains (nth 1 expr
)
885 (not (math-equal-int (nth 1 expr
) 1)))
886 (and (setq math-t1
(math-integral (math-div 1 (nth 2 expr
))))
887 (math-mul (nth 1 expr
) math-t1
)))
888 ((not (math-expr-contains (nth 2 expr
) math-integ-var
))
889 (and (setq math-t1
(math-integral (nth 1 expr
)))
890 (math-div math-t1
(nth 2 expr
))))
891 ((and (eq (car-safe (nth 1 expr
)) '*)
892 (not (math-expr-contains (nth 1 (nth 1 expr
))
894 (and (setq math-t1
(math-integral
895 (math-div (nth 2 (nth 1 expr
))
897 (math-mul math-t1
(nth 1 (nth 1 expr
)))))
898 ((and (eq (car-safe (nth 1 expr
)) '*)
899 (not (math-expr-contains (nth 2 (nth 1 expr
))
901 (and (setq math-t1
(math-integral
902 (math-div (nth 1 (nth 1 expr
))
904 (math-mul math-t1
(nth 2 (nth 1 expr
)))))
905 ((and (eq (car-safe (nth 2 expr
)) '*)
906 (not (math-expr-contains (nth 1 (nth 2 expr
))
908 (and (setq math-t1
(math-integral
909 (math-div (nth 1 expr
)
910 (nth 2 (nth 2 expr
)))))
911 (math-div math-t1
(nth 1 (nth 2 expr
)))))
912 ((and (eq (car-safe (nth 2 expr
)) '*)
913 (not (math-expr-contains (nth 2 (nth 2 expr
))
915 (and (setq math-t1
(math-integral
916 (math-div (nth 1 expr
)
917 (nth 1 (nth 2 expr
)))))
918 (math-div math-t1
(nth 2 (nth 2 expr
)))))
919 ((eq (car-safe (nth 2 expr
)) 'calcFunc-exp
)
921 (math-mul (nth 1 expr
)
923 (math-neg (nth 1 (nth 2 expr
)))))))))
925 (cond ((not (math-expr-contains (nth 1 expr
) math-integ-var
))
926 (or (and (setq math-t1
(math-is-polynomial (nth 2 expr
)
929 (math-mul (nth 1 math-t1
)
935 (math-mul (nth 2 expr
)
940 ((not (math-expr-contains (nth 2 expr
) math-integ-var
))
941 (if (and (integerp (nth 2 expr
)) (< (nth 2 expr
) 0))
943 (list '/ 1 (math-pow (nth 1 expr
) (- (nth 2 expr
))))
945 (or (and (setq math-t1
(math-is-polynomial (nth 1 expr
)
948 (setq math-t2
(math-add (nth 2 expr
) 1))
949 (math-div (math-pow (nth 1 expr
) math-t2
)
950 (math-mul math-t2
(nth 1 math-t1
))))
951 (and (Math-negp (nth 2 expr
))
954 (math-pow (nth 1 expr
)
960 ;; Integral of a polynomial.
961 (and (setq math-t1
(math-is-polynomial expr math-integ-var
20))
965 (if (setq accum
(math-add accum
966 (math-div (math-mul (car math-t1
)
971 math-t1
(cdr math-t1
))
975 ;; Try looking it up!
976 (cond ((= (length expr
) 2)
977 (and (symbolp (car expr
))
978 (setq math-t1
(get (car expr
) 'math-integral
))
981 (not (setq math-t2
(funcall (car math-t1
)
983 (setq math-t1
(cdr math-t1
)))
984 (and math-t2
(math-normalize math-t2
)))))
986 (and (symbolp (car expr
))
987 (setq math-t1
(get (car expr
) 'math-integral-2
))
990 (not (setq math-t2
(funcall (car math-t1
)
993 (setq math-t1
(cdr math-t1
)))
994 (and math-t2
(math-normalize math-t2
))))))
996 ;; Integral of a rational function.
997 (and (math-ratpoly-p expr math-integ-var
)
998 (setq math-t1
(calcFunc-apart expr math-integ-var
))
999 (not (equal math-t1 expr
))
1000 (math-integral math-t1
))
1002 ;; Try user-defined integration rules.
1004 (let ((math-old-integ (symbol-function 'calcFunc-integ
))
1005 (input (list 'calcFunc-integtry expr math-integ-var
))
1009 (fset 'calcFunc-integ
'math-sub-integration
)
1010 (setq res
(math-rewrite input
1011 '(var IntegRules var-IntegRules
)
1013 (fset 'calcFunc-integ math-old-integ
)
1014 (and (not (equal res input
))
1015 (if (setq part
(math-expr-calls
1016 res
'(calcFunc-integsubst)))
1017 (and (memq (length part
) '(3 4 5))
1018 (let ((parts (mapcar
1025 (math-integrate-by-substitution
1028 (list 'calcFunc-integfailed
1031 (if (not (math-expr-calls res
1033 calcFunc-integfailed
)))
1035 (fset 'calcFunc-integ math-old-integ
))))
1037 ;; See if the function is a symbolic derivative.
1038 (and (string-match "'" (symbol-name (car expr
)))
1039 (let ((name (symbol-name (car expr
)))
1040 (p expr
) (n 0) (which nil
) (bad nil
))
1041 (while (setq n
(1+ n
) p
(cdr p
))
1042 (if (equal (car p
) math-integ-var
)
1043 (if which
(setq bad t
) (setq which n
))
1044 (if (math-expr-contains (car p
) math-integ-var
)
1046 (and which
(not bad
)
1047 (let ((prime (if (= which
1) "'" (format "'%d" which
))))
1048 (and (string-match (concat prime
"\\('['0-9]*\\|$\\)")
1052 (substring name
0 (match-beginning 0))
1053 (substring name
(+ (match-beginning 0)
1057 ;; Try transformation methods (parts, substitutions).
1058 (and (> math-integ-level
0)
1059 (math-do-integral-methods expr
))
1061 ;; Try expanding the function's definition.
1062 (let ((res (math-expand-formula expr
)))
1064 (math-integral res
))))))
1066 (defun math-sub-integration (expr &rest rest
)
1067 (or (if (or (not rest
)
1068 (and (< math-integ-level math-integral-limit
)
1069 (eq (car rest
) math-integ-var
)))
1070 (math-integral expr
)
1071 (let ((res (apply math-old-integ expr rest
)))
1072 (and (or (= math-integ-level math-integral-limit
)
1073 (not (math-expr-calls res
'calcFunc-integ
)))
1075 (list 'calcFunc-integfailed expr
)))
1077 ;; math-so-far is a local variable for math-do-integral-methods, but
1078 ;; is used by math-integ-try-linear-substitutions and
1079 ;; math-integ-try-substitutions.
1080 (defvar math-so-far
)
1082 ;; math-integ-expr is a local variable for math-do-integral-methods,
1083 ;; but is used by math-integ-try-linear-substitutions and
1084 ;; math-integ-try-substitutions.
1085 (defvar math-integ-expr
)
1087 (defun math-do-integral-methods (math-integ-expr)
1088 (let ((math-so-far math-integ-var-list-list
)
1091 ;; Integration by substitution, for various likely sub-expressions.
1092 ;; (In first pass, we look only for sub-exprs that are linear in X.)
1093 (or (math-integ-try-linear-substitutions math-integ-expr
)
1094 (math-integ-try-substitutions math-integ-expr
)
1096 ;; If function has sines and cosines, try tan(x/2) substitution.
1097 (and (let ((p (setq rat-in
(math-expr-rational-in math-integ-expr
))))
1099 (memq (car (car p
)) '(calcFunc-sin
1105 (equal (nth 1 (car p
)) math-integ-var
))
1108 (or (and (math-integ-parts-easy math-integ-expr
)
1109 (math-integ-try-parts math-integ-expr t
))
1110 (math-integrate-by-good-substitution
1111 math-integ-expr
(list 'calcFunc-tan
(math-div math-integ-var
2)))))
1113 ;; If function has sinh and cosh, try tanh(x/2) substitution.
1114 (and (let ((p rat-in
))
1116 (memq (car (car p
)) '(calcFunc-sinh
1123 (equal (nth 1 (car p
)) math-integ-var
))
1126 (or (and (math-integ-parts-easy math-integ-expr
)
1127 (math-integ-try-parts math-integ-expr t
))
1128 (math-integrate-by-good-substitution
1129 math-integ-expr
(list 'calcFunc-tanh
(math-div math-integ-var
2)))))
1131 ;; If function has square roots, try sin, tan, or sec substitution.
1132 (and (let ((p rat-in
))
1135 (or (equal (car p
) math-integ-var
)
1136 (and (eq (car (car p
)) 'calcFunc-sqrt
)
1137 (setq math-t1
(math-is-polynomial
1138 (nth 1 (setq math-t2
(car p
)))
1139 math-integ-var
2)))))
1141 (and (null p
) math-t1
))
1142 (if (cdr (cdr math-t1
))
1143 (if (math-guess-if-neg (nth 2 math-t1
))
1144 (let* ((c (math-sqrt (math-neg (nth 2 math-t1
))))
1145 (d (math-div (nth 1 math-t1
) (math-mul -
2 c
)))
1146 (a (math-sqrt (math-add (car math-t1
) (math-sqr d
)))))
1147 (math-integrate-by-good-substitution
1148 math-integ-expr
(list 'calcFunc-arcsin
1150 (math-add (math-mul c math-integ-var
) d
)
1152 (let* ((c (math-sqrt (nth 2 math-t1
)))
1153 (d (math-div (nth 1 math-t1
) (math-mul 2 c
)))
1154 (aa (math-sub (car math-t1
) (math-sqr d
))))
1155 (if (and nil
(not (and (eq d
0) (eq c
1))))
1156 (math-integrate-by-good-substitution
1157 math-integ-expr
(math-add (math-mul c math-integ-var
) d
))
1158 (if (math-guess-if-neg aa
)
1159 (math-integrate-by-good-substitution
1160 math-integ-expr
(list 'calcFunc-arccosh
1162 (math-add (math-mul c math-integ-var
)
1164 (math-sqrt (math-neg aa
)))))
1165 (math-integrate-by-good-substitution
1166 math-integ-expr
(list 'calcFunc-arcsinh
1168 (math-add (math-mul c math-integ-var
)
1170 (math-sqrt aa
))))))))
1171 (math-integrate-by-good-substitution math-integ-expr math-t2
)) )
1173 ;; Try integration by parts.
1174 (math-integ-try-parts math-integ-expr
)
1179 (defun math-integ-parts-easy (expr)
1180 (cond ((Math-primp expr
) t
)
1181 ((memq (car expr
) '(+ -
*))
1182 (and (math-integ-parts-easy (nth 1 expr
))
1183 (math-integ-parts-easy (nth 2 expr
))))
1185 (and (math-integ-parts-easy (nth 1 expr
))
1186 (math-atomic-factorp (nth 2 expr
))))
1188 (and (natnump (nth 2 expr
))
1189 (math-integ-parts-easy (nth 1 expr
))))
1190 ((eq (car expr
) 'neg
)
1191 (math-integ-parts-easy (nth 1 expr
)))
1194 ;; math-prev-parts-v is local to calcFunc-integ (as well as
1195 ;; math-integrate-by-parts), but is used by math-integ-try-parts.
1196 (defvar math-prev-parts-v
)
1198 ;; math-good-parts is local to calcFunc-integ (as well as
1199 ;; math-integ-try-parts), but is used by math-integrate-by-parts.
1200 (defvar math-good-parts
)
1203 (defun math-integ-try-parts (expr &optional math-good-parts
)
1204 ;; Integration by parts:
1205 ;; integ(f(x) g(x),x) = f(x) h(x) - integ(h(x) f'(x),x)
1206 ;; where h(x) = integ(g(x),x).
1207 (or (let ((exp (calcFunc-expand expr
)))
1208 (and (not (equal exp expr
))
1209 (math-integral exp
)))
1210 (and (eq (car expr
) '*)
1211 (let ((first-bad (or (math-polynomial-p (nth 1 expr
)
1213 (equal (nth 2 expr
) math-prev-parts-v
))))
1214 (or (and first-bad
; so try this one first
1215 (math-integrate-by-parts (nth 1 expr
) (nth 2 expr
)))
1216 (math-integrate-by-parts (nth 2 expr
) (nth 1 expr
))
1217 (and (not first-bad
)
1218 (math-integrate-by-parts (nth 1 expr
) (nth 2 expr
))))))
1219 (and (eq (car expr
) '/)
1220 (math-expr-contains (nth 1 expr
) math-integ-var
)
1221 (let ((recip (math-div 1 (nth 2 expr
))))
1222 (or (math-integrate-by-parts (nth 1 expr
) recip
)
1223 (math-integrate-by-parts recip
(nth 1 expr
)))))
1224 (and (eq (car expr
) '^
)
1225 (math-integrate-by-parts (math-pow (nth 1 expr
)
1226 (math-sub (nth 2 expr
) 1))
1229 (defun math-integrate-by-parts (u vprime
)
1230 (let ((math-integ-level (if (or math-good-parts
1231 (math-polynomial-p u math-integ-var
))
1233 (1- math-integ-level
)))
1234 (math-doing-parts t
)
1236 (and (>= math-integ-level
0)
1239 (setcar (cdr math-cur-record
) 'parts
)
1240 (math-tracing-integral "Integrating by parts, u = "
1241 (math-format-value u
1000)
1243 (math-format-value vprime
1000)
1245 (and (setq v
(math-integral vprime
))
1246 (setq temp
(calcFunc-deriv u math-integ-var nil t
))
1247 (setq temp
(let ((math-prev-parts-v v
))
1248 (math-integral (math-mul v temp
) 'yes
)))
1249 (setq temp
(math-sub (math-mul u v
) temp
))
1250 (if (eq (nth 1 math-cur-record
) 'parts
)
1251 (calcFunc-expand temp
)
1252 (setq v
(list 'var
'PARTS math-cur-record
)
1253 temp
(let (calc-next-why)
1254 (math-simplify-extended
1255 (math-solve-for (math-sub v temp
) 0 v nil
)))
1256 temp
(if (and (eq (car-safe temp
) '/)
1257 (math-zerop (nth 2 temp
)))
1259 (setcar (cdr math-cur-record
) 'busy
)))))
1261 ;;; This tries two different formulations, hoping the algebraic simplifier
1262 ;;; will be strong enough to handle at least one.
1263 (defun math-integrate-by-substitution (expr u
&optional user uinv uinvprime
)
1264 (and (> math-integ-level
0)
1265 (let ((math-integ-level (max (- math-integ-level
2) 0)))
1266 (math-integrate-by-good-substitution expr u user uinv uinvprime
))))
1268 (defun math-integrate-by-good-substitution (expr u
&optional user
1270 (let ((math-living-dangerously t
)
1272 (and (setq uinv
(if uinv
1273 (math-expr-subst uinv math-integ-var
1275 (let (calc-next-why)
1278 math-integ-var nil
))))
1280 (math-tracing-integral "Integrating by substitution, u = "
1281 (math-format-value u
1000)
1283 (or (and (setq deriv
(calcFunc-deriv u
1286 (setq temp
(math-integral (math-expr-subst
1289 (math-div expr deriv
)
1297 (and (setq deriv
(or uinvprime
1298 (calcFunc-deriv uinv
1302 (setq temp
(math-integral (math-mul
1315 (math-simplify-extended
1316 (math-expr-subst temp math-integ-var u
)))))
1318 ;;; Look for substitutions of the form u = a x + b.
1319 (defun math-integ-try-linear-substitutions (sub-expr)
1320 (setq math-linear-subst-tried t
)
1321 (and (not (Math-primp sub-expr
))
1322 (or (and (not (memq (car sub-expr
) '(+ -
* / neg
)))
1323 (not (and (eq (car sub-expr
) '^
)
1324 (integerp (nth 2 sub-expr
))))
1325 (math-expr-contains sub-expr math-integ-var
)
1327 (while (and (setq sub-expr
(cdr sub-expr
))
1328 (or (not (math-linear-in (car sub-expr
)
1330 (assoc (car sub-expr
) math-so-far
)
1332 (setq math-so-far
(cons (list (car sub-expr
))
1335 (math-integrate-by-substitution
1336 math-integ-expr
(car sub-expr
))))))))
1339 (while (and (setq sub-expr
(cdr sub-expr
))
1340 (not (setq res
(math-integ-try-linear-substitutions
1344 ;;; Recursively try different substitutions based on various sub-expressions.
1345 (defun math-integ-try-substitutions (sub-expr &optional allow-rat
)
1346 (and (not (Math-primp sub-expr
))
1347 (not (assoc sub-expr math-so-far
))
1348 (math-expr-contains sub-expr math-integ-var
)
1349 (or (and (if (and (not (memq (car sub-expr
) '(+ -
* / neg
)))
1350 (not (and (eq (car sub-expr
) '^
)
1351 (integerp (nth 2 sub-expr
)))))
1353 (prog1 allow-rat
(setq allow-rat nil
)))
1354 (not (eq sub-expr math-integ-expr
))
1355 (or (math-integrate-by-substitution math-integ-expr sub-expr
)
1356 (and (eq (car sub-expr
) '^
)
1357 (integerp (nth 2 sub-expr
))
1358 (< (nth 2 sub-expr
) 0)
1359 (math-integ-try-substitutions
1360 (math-pow (nth 1 sub-expr
) (- (nth 2 sub-expr
)))
1363 (setq math-so-far
(cons (list sub-expr
) math-so-far
))
1364 (while (and (setq sub-expr
(cdr sub-expr
))
1365 (not (setq res
(math-integ-try-substitutions
1366 (car sub-expr
) allow-rat
)))))
1369 ;; The variable math-expr-parts is local to math-expr-rational-in,
1370 ;; but is used by math-expr-rational-in-rec
1371 (defvar math-expr-parts
)
1373 (defun math-expr-rational-in (expr)
1374 (let ((math-expr-parts nil
))
1375 (math-expr-rational-in-rec expr
)
1376 (mapcar 'car math-expr-parts
)))
1378 (defun math-expr-rational-in-rec (expr)
1379 (cond ((Math-primp expr
)
1380 (and (equal expr math-integ-var
)
1381 (not (assoc expr math-expr-parts
))
1382 (setq math-expr-parts
(cons (list expr
) math-expr-parts
))))
1383 ((or (memq (car expr
) '(+ -
* / neg
))
1384 (and (eq (car expr
) '^
) (integerp (nth 2 expr
))))
1385 (math-expr-rational-in-rec (nth 1 expr
))
1386 (and (nth 2 expr
) (math-expr-rational-in-rec (nth 2 expr
))))
1387 ((and (eq (car expr
) '^
)
1388 (eq (math-quarter-integer (nth 2 expr
)) 2))
1389 (math-expr-rational-in-rec (list 'calcFunc-sqrt
(nth 1 expr
))))
1391 (and (not (assoc expr math-expr-parts
))
1392 (math-expr-contains expr math-integ-var
)
1393 (setq math-expr-parts
(cons (list expr
) math-expr-parts
))))))
1395 (defun math-expr-calls (expr funcs
&optional arg-contains
)
1397 (if (or (memq (car expr
) funcs
)
1398 (and (eq (car expr
) '^
) (eq (car funcs
) 'calcFunc-sqrt
)
1399 (eq (math-quarter-integer (nth 2 expr
)) 2)))
1400 (and (or (not arg-contains
)
1401 (math-expr-contains expr arg-contains
))
1403 (and (not (Math-primp expr
))
1405 (while (and (setq expr
(cdr expr
))
1406 (not (setq res
(math-expr-calls
1407 (car expr
) funcs arg-contains
)))))
1410 (defun math-fix-const-terms (expr except-vars
)
1411 (cond ((not (math-expr-depends expr except-vars
)) 0)
1412 ((Math-primp expr
) expr
)
1414 (math-add (math-fix-const-terms (nth 1 expr
) except-vars
)
1415 (math-fix-const-terms (nth 2 expr
) except-vars
)))
1417 (math-sub (math-fix-const-terms (nth 1 expr
) except-vars
)
1418 (math-fix-const-terms (nth 2 expr
) except-vars
)))
1421 ;; Command for debugging the Calculator's symbolic integrator.
1422 (defun calc-dump-integral-cache (&optional arg
)
1424 (let ((buf (current-buffer)))
1426 (let ((p math-integral-cache
)
1428 (display-buffer (get-buffer-create "*Integral Cache*"))
1429 (set-buffer (get-buffer "*Integral Cache*"))
1432 (setq math-cur-record
(car p
))
1433 (or arg
(math-replace-integral-parts math-cur-record
))
1434 (insert (math-format-flat-expr (car math-cur-record
) 0)
1436 (if (symbolp (nth 1 math-cur-record
))
1437 (concat "(" (symbol-name (nth 1 math-cur-record
)) ")")
1438 (math-format-flat-expr (nth 1 math-cur-record
) 0))
1441 (goto-char (point-min)))
1444 ;; The variable math-max-integral-limit is local to calcFunc-integ,
1445 ;; but is used by math-try-integral.
1446 (defvar math-max-integral-limit
)
1448 (defun math-try-integral (expr)
1449 (let ((math-integ-level math-integral-limit
)
1450 (math-integ-depth 0)
1451 (math-integ-msg "Working...done")
1452 (math-cur-record nil
) ; a technicality
1453 (math-integrating t
)
1454 (calc-prefer-frac t
)
1455 (calc-symbolic-mode t
)
1456 (math-has-rules (calc-has-rules 'var-IntegRules
)))
1457 (or (math-integral expr
'yes
)
1458 (and math-any-substs
1459 (setq math-enable-subst t
)
1460 (math-integral expr
'yes
))
1461 (and (> math-max-integral-limit math-integral-limit
)
1462 (setq math-integral-limit math-max-integral-limit
1463 math-integ-level math-integral-limit
)
1464 (math-integral expr
'yes
)))))
1466 (defvar var-IntegLimit nil
)
1468 (defun calcFunc-integ (expr var
&optional low high
)
1470 ;; Do these even if the parts turn out not to be integrable.
1471 ((eq (car-safe expr
) '+)
1472 (math-add (calcFunc-integ (nth 1 expr
) var low high
)
1473 (calcFunc-integ (nth 2 expr
) var low high
)))
1474 ((eq (car-safe expr
) '-
)
1475 (math-sub (calcFunc-integ (nth 1 expr
) var low high
)
1476 (calcFunc-integ (nth 2 expr
) var low high
)))
1477 ((eq (car-safe expr
) 'neg
)
1478 (math-neg (calcFunc-integ (nth 1 expr
) var low high
)))
1479 ((and (eq (car-safe expr
) '*)
1480 (not (math-expr-contains (nth 1 expr
) var
)))
1481 (math-mul (nth 1 expr
) (calcFunc-integ (nth 2 expr
) var low high
)))
1482 ((and (eq (car-safe expr
) '*)
1483 (not (math-expr-contains (nth 2 expr
) var
)))
1484 (math-mul (calcFunc-integ (nth 1 expr
) var low high
) (nth 2 expr
)))
1485 ((and (eq (car-safe expr
) '/)
1486 (not (math-expr-contains (nth 1 expr
) var
))
1487 (not (math-equal-int (nth 1 expr
) 1)))
1488 (math-mul (nth 1 expr
)
1489 (calcFunc-integ (math-div 1 (nth 2 expr
)) var low high
)))
1490 ((and (eq (car-safe expr
) '/)
1491 (not (math-expr-contains (nth 2 expr
) var
)))
1492 (math-div (calcFunc-integ (nth 1 expr
) var low high
) (nth 2 expr
)))
1493 ((and (eq (car-safe expr
) '/)
1494 (eq (car-safe (nth 1 expr
)) '*)
1495 (not (math-expr-contains (nth 1 (nth 1 expr
)) var
)))
1496 (math-mul (nth 1 (nth 1 expr
))
1497 (calcFunc-integ (math-div (nth 2 (nth 1 expr
)) (nth 2 expr
))
1499 ((and (eq (car-safe expr
) '/)
1500 (eq (car-safe (nth 1 expr
)) '*)
1501 (not (math-expr-contains (nth 2 (nth 1 expr
)) var
)))
1502 (math-mul (nth 2 (nth 1 expr
))
1503 (calcFunc-integ (math-div (nth 1 (nth 1 expr
)) (nth 2 expr
))
1505 ((and (eq (car-safe expr
) '/)
1506 (eq (car-safe (nth 2 expr
)) '*)
1507 (not (math-expr-contains (nth 1 (nth 2 expr
)) var
)))
1508 (math-div (calcFunc-integ (math-div (nth 1 expr
) (nth 2 (nth 2 expr
)))
1510 (nth 1 (nth 2 expr
))))
1511 ((and (eq (car-safe expr
) '/)
1512 (eq (car-safe (nth 2 expr
)) '*)
1513 (not (math-expr-contains (nth 2 (nth 2 expr
)) var
)))
1514 (math-div (calcFunc-integ (math-div (nth 1 expr
) (nth 1 (nth 2 expr
)))
1516 (nth 2 (nth 2 expr
))))
1517 ((eq (car-safe expr
) 'vec
)
1518 (cons 'vec
(mapcar (function (lambda (x) (calcFunc-integ x var low high
)))
1521 (let ((state (list calc-angle-mode
1522 ;;calc-symbolic-mode
1525 (calc-var-value 'var-IntegRules
)
1526 (calc-var-value 'var-IntegSimpRules
))))
1527 (or (equal state math-integral-cache-state
)
1528 (setq math-integral-cache-state state
1529 math-integral-cache nil
)))
1530 (let* ((math-max-integral-limit (or (and (natnump var-IntegLimit
)
1533 (math-integral-limit 1)
1534 (sexpr (math-expr-subst expr var math-integ-var
))
1535 (trace-buffer (get-buffer "*Trace*"))
1536 (calc-language (if (eq calc-language
'big
) nil calc-language
))
1538 (math-enable-subst nil
)
1539 (math-prev-parts-v nil
)
1540 (math-doing-parts nil
)
1541 (math-good-parts nil
)
1544 (let ((calcbuf (current-buffer))
1545 (calcwin (selected-window)))
1548 (if (get-buffer-window trace-buffer
)
1549 (select-window (get-buffer-window trace-buffer
)))
1550 (set-buffer trace-buffer
)
1551 (goto-char (point-max))
1552 (or (assq 'scroll-stop
(buffer-local-variables))
1554 (make-local-variable 'scroll-step
)
1555 (setq scroll-step
3)))
1557 (set-buffer calcbuf
)
1558 (math-try-integral sexpr
))
1559 (select-window calcwin
)
1560 (set-buffer calcbuf
)))
1561 (math-try-integral sexpr
))))
1564 (if (calc-has-rules 'var-IntegAfterRules
)
1565 (setq res
(math-rewrite res
'(var IntegAfterRules
1566 var-IntegAfterRules
))))
1569 (math-sub (math-expr-subst res math-integ-var high
)
1570 (math-expr-subst res math-integ-var low
))
1571 (setq res
(math-fix-const-terms res math-integ-vars
))
1573 (math-expr-subst res math-integ-var low
)
1574 (math-expr-subst res math-integ-var var
)))))
1575 (append (list 'calcFunc-integ expr var
)
1576 (and low
(list low
))
1577 (and high
(list high
))))))))
1580 (math-defintegral calcFunc-inv
1581 (math-integral (math-div 1 u
)))
1583 (math-defintegral calcFunc-conj
1584 (let ((int (math-integral u
)))
1586 (list 'calcFunc-conj int
))))
1588 (math-defintegral calcFunc-deg
1589 (let ((int (math-integral u
)))
1591 (list 'calcFunc-deg int
))))
1593 (math-defintegral calcFunc-rad
1594 (let ((int (math-integral u
)))
1596 (list 'calcFunc-rad int
))))
1598 (math-defintegral calcFunc-re
1599 (let ((int (math-integral u
)))
1601 (list 'calcFunc-re int
))))
1603 (math-defintegral calcFunc-im
1604 (let ((int (math-integral u
)))
1606 (list 'calcFunc-im int
))))
1608 (math-defintegral calcFunc-sqrt
1609 (and (equal u math-integ-var
)
1610 (math-mul '(frac 2 3)
1611 (list 'calcFunc-sqrt
(math-pow u
3)))))
1613 (math-defintegral calcFunc-exp
1614 (or (and (equal u math-integ-var
)
1615 (list 'calcFunc-exp u
))
1616 (let ((p (math-is-polynomial u math-integ-var
2)))
1618 (let ((sqa (math-sqrt (math-neg (nth 2 p
)))))
1621 (math-mul (math-div (list 'calcFunc-sqrt
'(var pi var-pi
))
1625 (math-div (math-sub (math-mul (car p
)
1628 (math-sqr (nth 1 p
))
1632 (math-sub (math-mul sqa math-integ-var
)
1633 (math-div (nth 1 p
) (math-mul 2 sqa
)))))
1636 (math-defintegral calcFunc-ln
1637 (or (and (equal u math-integ-var
)
1638 (math-sub (math-mul u
(list 'calcFunc-ln u
)) u
))
1639 (and (eq (car u
) '*)
1640 (math-integral (math-add (list 'calcFunc-ln
(nth 1 u
))
1641 (list 'calcFunc-ln
(nth 2 u
)))))
1642 (and (eq (car u
) '/)
1643 (math-integral (math-sub (list 'calcFunc-ln
(nth 1 u
))
1644 (list 'calcFunc-ln
(nth 2 u
)))))
1645 (and (eq (car u
) '^
)
1646 (math-integral (math-mul (nth 2 u
)
1647 (list 'calcFunc-ln
(nth 1 u
)))))))
1649 (math-defintegral calcFunc-log10
1650 (and (equal u math-integ-var
)
1651 (math-sub (math-mul u
(list 'calcFunc-ln u
))
1652 (math-div u
(list 'calcFunc-ln
10)))))
1654 (math-defintegral-2 calcFunc-log
1655 (math-integral (math-div (list 'calcFunc-ln u
)
1656 (list 'calcFunc-ln v
))))
1658 (math-defintegral calcFunc-sin
1659 (or (and (equal u math-integ-var
)
1660 (math-neg (math-from-radians-2 (list 'calcFunc-cos u
))))
1661 (and (nth 2 (math-is-polynomial u math-integ-var
2))
1662 (math-integral (math-to-exponentials (list 'calcFunc-sin u
))))))
1664 (math-defintegral calcFunc-cos
1665 (or (and (equal u math-integ-var
)
1666 (math-from-radians-2 (list 'calcFunc-sin u
)))
1667 (and (nth 2 (math-is-polynomial u math-integ-var
2))
1668 (math-integral (math-to-exponentials (list 'calcFunc-cos u
))))))
1670 (math-defintegral calcFunc-tan
1671 (and (equal u math-integ-var
)
1672 (math-from-radians-2
1673 (list 'calcFunc-ln
(list 'calcFunc-sec u
)))))
1675 (math-defintegral calcFunc-sec
1676 (and (equal u math-integ-var
)
1677 (math-from-radians-2
1680 (list 'calcFunc-sec u
)
1681 (list 'calcFunc-tan u
))))))
1683 (math-defintegral calcFunc-csc
1684 (and (equal u math-integ-var
)
1685 (math-from-radians-2
1688 (list 'calcFunc-csc u
)
1689 (list 'calcFunc-cot u
))))))
1691 (math-defintegral calcFunc-cot
1692 (and (equal u math-integ-var
)
1693 (math-from-radians-2
1694 (list 'calcFunc-ln
(list 'calcFunc-sin u
)))))
1696 (math-defintegral calcFunc-arcsin
1697 (and (equal u math-integ-var
)
1698 (math-add (math-mul u
(list 'calcFunc-arcsin u
))
1699 (math-from-radians-2
1700 (list 'calcFunc-sqrt
(math-sub 1 (math-sqr u
)))))))
1702 (math-defintegral calcFunc-arccos
1703 (and (equal u math-integ-var
)
1704 (math-sub (math-mul u
(list 'calcFunc-arccos u
))
1705 (math-from-radians-2
1706 (list 'calcFunc-sqrt
(math-sub 1 (math-sqr u
)))))))
1708 (math-defintegral calcFunc-arctan
1709 (and (equal u math-integ-var
)
1710 (math-sub (math-mul u
(list 'calcFunc-arctan u
))
1711 (math-from-radians-2
1712 (math-div (list 'calcFunc-ln
(math-add 1 (math-sqr u
)))
1715 (math-defintegral calcFunc-sinh
1716 (and (equal u math-integ-var
)
1717 (list 'calcFunc-cosh u
)))
1719 (math-defintegral calcFunc-cosh
1720 (and (equal u math-integ-var
)
1721 (list 'calcFunc-sinh u
)))
1723 (math-defintegral calcFunc-tanh
1724 (and (equal u math-integ-var
)
1725 (list 'calcFunc-ln
(list 'calcFunc-cosh u
))))
1727 (math-defintegral calcFunc-sech
1728 (and (equal u math-integ-var
)
1729 (list 'calcFunc-arctan
(list 'calcFunc-sinh u
))))
1731 (math-defintegral calcFunc-csch
1732 (and (equal u math-integ-var
)
1733 (list 'calcFunc-ln
(list 'calcFunc-tanh
(math-div u
2)))))
1735 (math-defintegral calcFunc-coth
1736 (and (equal u math-integ-var
)
1737 (list 'calcFunc-ln
(list 'calcFunc-sinh u
))))
1739 (math-defintegral calcFunc-arcsinh
1740 (and (equal u math-integ-var
)
1741 (math-sub (math-mul u
(list 'calcFunc-arcsinh u
))
1742 (list 'calcFunc-sqrt
(math-add (math-sqr u
) 1)))))
1744 (math-defintegral calcFunc-arccosh
1745 (and (equal u math-integ-var
)
1746 (math-sub (math-mul u
(list 'calcFunc-arccosh u
))
1747 (list 'calcFunc-sqrt
(math-sub 1 (math-sqr u
))))))
1749 (math-defintegral calcFunc-arctanh
1750 (and (equal u math-integ-var
)
1751 (math-sub (math-mul u
(list 'calcFunc-arctan u
))
1752 (math-div (list 'calcFunc-ln
1753 (math-add 1 (math-sqr u
)))
1756 ;;; (Ax + B) / (ax^2 + bx + c)^n forms.
1757 (math-defintegral-2 /
1758 (math-integral-rational-funcs u v
))
1760 (defun math-integral-rational-funcs (u v
)
1761 (let ((pu (math-is-polynomial u math-integ-var
1))
1765 (if (and (eq (car-safe v
) '^
) (natnump (nth 2 v
)))
1766 (setq vpow
(nth 2 v
)
1768 (and (setq pv
(math-is-polynomial v math-integ-var
2))
1769 (let ((int (math-mul-thru
1771 (math-integral-q02 (car pv
) (nth 1 pv
)
1772 (nth 2 pv
) v vpow
))))
1774 (setq int
(math-add int
1779 (nth 2 pv
) v vpow
)))))
1782 (defun math-integral-q12 (a b c v vpow
)
1786 (math-sub (math-div math-integ-var b
)
1787 (math-mul (math-div a
(math-sqr b
))
1788 (list 'calcFunc-ln v
))))
1790 (math-div (math-add (list 'calcFunc-ln v
)
1794 (let ((nm1 (math-sub vpow
1))
1795 (nm2 (math-sub vpow
2)))
1797 (math-div a
(math-mul nm1
(math-pow v nm1
)))
1798 (math-div 1 (math-mul nm2
(math-pow v nm2
))))
1801 (setq q
(math-sub (math-mul 4 (math-mul a c
)) (math-sqr b
))))
1802 (let ((part (math-div b
(math-mul 2 c
))))
1803 (math-mul-thru (math-pow c vpow
)
1804 (math-integral-q12 part
1 nil
1805 (math-add math-integ-var part
)
1808 (and (math-ratp q
) (math-negp q
)
1809 (let ((calc-symbolic-mode t
))
1810 (math-ratp (math-sqrt (math-neg q
))))
1811 (throw 'int-rat nil
)) ; should have used calcFunc-apart first
1812 (math-sub (math-div (list 'calcFunc-ln v
) (math-mul 2 c
))
1813 (math-mul-thru (math-div b
(math-mul 2 c
))
1814 (math-integral-q02 a b c v
1))))
1816 (let ((n (1- vpow
)))
1817 (math-sub (math-neg (math-div
1818 (math-add (math-mul b math-integ-var
)
1820 (math-mul n
(math-mul q
(math-pow v n
)))))
1821 (math-mul-thru (math-div (math-mul b
(1- (* 2 n
)))
1823 (math-integral-q02 a b c v n
))))))))
1825 (defun math-integral-q02 (a b c v vpow
)
1829 (math-div (list 'calcFunc-ln v
) b
))
1831 (math-div (math-pow v
(- 1 vpow
))
1832 (math-mul (- 1 vpow
) b
)))))
1834 (setq q
(math-sub (math-mul 4 (math-mul a c
)) (math-sqr b
))))
1835 (let ((part (math-div b
(math-mul 2 c
))))
1836 (math-mul-thru (math-pow c vpow
)
1837 (math-integral-q02 part
1 nil
1838 (math-add math-integ-var part
)
1841 (setq part
(math-add (math-mul 2 (math-mul c math-integ-var
)) b
))
1843 (let ((n (1- vpow
)))
1844 (math-add (math-div part
(math-mul n
(math-mul q
(math-pow v n
))))
1845 (math-mul-thru (math-div (math-mul (- (* 4 n
) 2) c
)
1847 (math-integral-q02 a b c v n
)))))
1848 ((math-guess-if-neg q
)
1849 (setq rq
(list 'calcFunc-sqrt
(math-neg q
)))
1850 ;;(math-div-thru (list 'calcFunc-ln
1851 ;; (math-div (math-sub part rq)
1852 ;; (math-add part rq)))
1854 (math-div (math-mul -
2 (list 'calcFunc-arctanh
1855 (math-div part rq
)))
1858 (setq rq
(list 'calcFunc-sqrt q
))
1859 (math-div (math-mul 2 (math-to-radians-2
1860 (list 'calcFunc-arctan
1861 (math-div part rq
))))
1865 (math-defintegral calcFunc-erf
1866 (and (equal u math-integ-var
)
1867 (math-add (math-mul u
(list 'calcFunc-erf u
))
1868 (math-div 1 (math-mul (list 'calcFunc-exp
(math-sqr u
))
1869 (list 'calcFunc-sqrt
1870 '(var pi var-pi
)))))))
1872 (math-defintegral calcFunc-erfc
1873 (and (equal u math-integ-var
)
1874 (math-sub (math-mul u
(list 'calcFunc-erfc u
))
1875 (math-div 1 (math-mul (list 'calcFunc-exp
(math-sqr u
))
1876 (list 'calcFunc-sqrt
1877 '(var pi var-pi
)))))))
1882 (defvar math-tabulate-initial nil
)
1883 (defvar math-tabulate-function nil
)
1885 ;; These variables are local to calcFunc-table, but are used by
1886 ;; math-scan-for-limits.
1891 (defun calcFunc-table (expr math-var
&optional calc-low calc-high step
)
1893 (setq calc-low
'(neg (var inf var-inf
)) calc-high
'(var inf var-inf
)))
1894 (or calc-high
(setq calc-high calc-low calc-low
1))
1895 (and (or (math-infinitep calc-low
) (math-infinitep calc-high
))
1897 (math-scan-for-limits expr
))
1898 (and step
(math-zerop step
) (math-reject-arg step
'nonzerop
))
1899 (let ((known (+ (if (Math-objectp calc-low
) 1 0)
1900 (if (Math-objectp calc-high
) 1 0)
1901 (if (or (null step
) (Math-objectp step
)) 1 0)))
1902 (count '(var inf var-inf
))
1904 (or (= known
2) ; handy optimization
1905 (equal calc-high
'(var inf var-inf
))
1907 (setq count
(math-div (math-sub calc-high calc-low
) (or step
1)))
1908 (or (Math-objectp count
)
1909 (setq count
(math-simplify count
)))
1910 (if (Math-messy-integerp count
)
1911 (setq count
(math-trunc count
)))))
1912 (if (Math-negp count
)
1914 (if (integerp count
)
1915 (let ((var-DUMMY nil
)
1916 (vec math-tabulate-initial
)
1917 (math-working-step-2 (1+ count
))
1918 (math-working-step 0))
1919 (setq expr
(math-evaluate-expr
1920 (math-expr-subst expr math-var
'(var DUMMY var-DUMMY
))))
1922 (setq math-working-step
(1+ math-working-step
)
1924 vec
(cond ((eq math-tabulate-function
'calcFunc-sum
)
1925 (math-add vec
(math-evaluate-expr expr
)))
1926 ((eq math-tabulate-function
'calcFunc-prod
)
1927 (math-mul vec
(math-evaluate-expr expr
)))
1929 (cons (math-evaluate-expr expr
) vec
)))
1930 calc-low
(math-add calc-low
(or step
1))
1932 (if math-tabulate-function
1934 (cons 'vec
(nreverse vec
))))
1935 (if (Math-integerp count
)
1936 (calc-record-why 'fixnump calc-high
)
1937 (if (Math-num-integerp calc-low
)
1938 (if (Math-num-integerp calc-high
)
1939 (calc-record-why 'integerp step
)
1940 (calc-record-why 'integerp calc-high
))
1941 (calc-record-why 'integerp calc-low
)))
1942 (append (list (or math-tabulate-function
'calcFunc-table
)
1944 (and (not (and (equal calc-low
'(neg (var inf var-inf
)))
1945 (equal calc-high
'(var inf var-inf
))))
1946 (list calc-low calc-high
))
1947 (and step
(list step
))))))
1949 (defun math-scan-for-limits (x)
1950 (cond ((Math-primp x
))
1951 ((and (eq (car x
) 'calcFunc-subscr
)
1952 (Math-vectorp (nth 1 x
))
1953 (math-expr-contains (nth 2 x
) math-var
))
1954 (let* ((calc-next-why nil
)
1955 (low-val (math-solve-for (nth 2 x
) 1 math-var nil
))
1956 (high-val (math-solve-for (nth 2 x
) (1- (length (nth 1 x
)))
1959 (and low-val
(math-realp low-val
)
1960 high-val
(math-realp high-val
))
1961 (and (Math-lessp high-val low-val
)
1962 (setq temp low-val low-val high-val high-val temp
))
1963 (setq calc-low
(math-max calc-low
(math-ceiling low-val
))
1964 calc-high
(math-min calc-high
(math-floor high-val
)))))
1966 (while (setq x
(cdr x
))
1967 (math-scan-for-limits (car x
))))))
1970 (defvar math-disable-sums nil
)
1971 (defun calcFunc-sum (expr var
&optional low high step
)
1972 (if math-disable-sums
(math-reject-arg))
1973 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec
2)))
1974 (math-sum-rec expr var low high step
)))
1975 (math-disable-sums t
))
1976 (math-normalize res
)))
1978 (defun math-sum-rec (expr var
&optional low high step
)
1979 (or low
(setq low
'(neg (var inf var-inf
)) high
'(var inf var-inf
)))
1980 (and low
(not high
) (setq high low low
1))
1984 ((not (math-expr-contains expr var
))
1985 (math-mul expr
(math-add (math-div (math-sub high low
) (or step
1))
1987 ((and step
(not (math-equal-int step
1)))
1988 (if (math-negp step
)
1989 (math-sum-rec expr var high low
(math-neg step
))
1990 (let ((lo (math-simplify (math-div low step
))))
1991 (if (math-known-num-integerp lo
)
1992 (math-sum-rec (math-normalize
1993 (math-expr-subst expr var
1994 (math-mul step var
)))
1995 var lo
(math-simplify (math-div high step
)))
1996 (math-sum-rec (math-normalize
1997 (math-expr-subst expr var
1998 (math-add (math-mul step var
)
2001 (math-simplify (math-div (math-sub high low
)
2003 ((memq (setq t1
(math-compare low high
)) '(0 1))
2005 (math-expr-subst expr var low
)
2007 ((setq t1
(math-is-polynomial expr var
20))
2011 (setq poly
(math-poly-mix poly
1
2012 (math-sum-integer-power n
) (car t1
))
2015 (setq n
(math-build-polynomial-expr poly high
))
2018 (math-sub n
(math-build-polynomial-expr poly
2019 (math-sub low
1))))))
2020 ((and (memq (car expr
) '(+ -
))
2021 (setq t1
(math-sum-rec (nth 1 expr
) var low high
)
2022 t2
(math-sum-rec (nth 2 expr
) var low high
))
2023 (not (and (math-expr-calls t1
'(calcFunc-sum))
2024 (math-expr-calls t2
'(calcFunc-sum)))))
2025 (list (car expr
) t1 t2
))
2026 ((and (eq (car expr
) '*)
2027 (setq t1
(math-sum-const-factors expr var
)))
2028 (math-mul (car t1
) (math-sum-rec (cdr t1
) var low high
)))
2029 ((and (eq (car expr
) '*) (memq (car-safe (nth 1 expr
)) '(+ -
)))
2030 (math-sum-rec (math-add-or-sub (math-mul (nth 1 (nth 1 expr
))
2032 (math-mul (nth 2 (nth 1 expr
))
2034 nil
(eq (car (nth 1 expr
)) '-
))
2036 ((and (eq (car expr
) '*) (memq (car-safe (nth 2 expr
)) '(+ -
)))
2037 (math-sum-rec (math-add-or-sub (math-mul (nth 1 expr
)
2038 (nth 1 (nth 2 expr
)))
2039 (math-mul (nth 1 expr
)
2040 (nth 2 (nth 2 expr
)))
2041 nil
(eq (car (nth 2 expr
)) '-
))
2043 ((and (eq (car expr
) '/)
2044 (not (math-primp (nth 1 expr
)))
2045 (setq t1
(math-sum-const-factors (nth 1 expr
) var
)))
2047 (math-sum-rec (math-div (cdr t1
) (nth 2 expr
))
2049 ((and (eq (car expr
) '/)
2050 (setq t1
(math-sum-const-factors (nth 2 expr
) var
)))
2051 (math-div (math-sum-rec (math-div (nth 1 expr
) (cdr t1
))
2054 ((eq (car expr
) 'neg
)
2055 (math-neg (math-sum-rec (nth 1 expr
) var low high
)))
2056 ((and (eq (car expr
) '^
)
2057 (not (math-expr-contains (nth 1 expr
) var
))
2058 (setq t1
(math-is-polynomial (nth 2 expr
) var
1)))
2059 (let ((x (math-pow (nth 1 expr
) (nth 1 t1
))))
2060 (math-div (math-mul (math-sub (math-pow x
(math-add 1 high
))
2062 (math-pow (nth 1 expr
) (car t1
)))
2064 ((and (setq t1
(math-to-exponentials expr
))
2065 (setq t1
(math-sum-rec t1 var low high
))
2066 (not (math-expr-calls t1
'(calcFunc-sum))))
2068 ((memq (car expr
) '(calcFunc-ln calcFunc-log10
))
2069 (list (car expr
) (calcFunc-prod (nth 1 expr
) var low high
)))
2070 ((and (eq (car expr
) 'calcFunc-log
)
2072 (not (math-expr-contains (nth 2 expr
) var
)))
2074 (calcFunc-prod (nth 1 expr
) var low high
)
2076 (if (equal val
'(var nan var-nan
)) (setq val nil
))
2078 (let* ((math-tabulate-initial 0)
2079 (math-tabulate-function 'calcFunc-sum
))
2080 (calcFunc-table expr var low high
)))))
2082 (defun calcFunc-asum (expr var low
&optional high step no-mul-flag
)
2083 (or high
(setq high low low
1))
2084 (if (and step
(not (math-equal-int step
1)))
2085 (if (math-negp step
)
2086 (math-mul (math-pow -
1 low
)
2087 (calcFunc-asum expr var high low
(math-neg step
) t
))
2088 (let ((lo (math-simplify (math-div low step
))))
2089 (if (math-num-integerp lo
)
2090 (calcFunc-asum (math-normalize
2091 (math-expr-subst expr var
2092 (math-mul step var
)))
2093 var lo
(math-simplify (math-div high step
)))
2094 (calcFunc-asum (math-normalize
2095 (math-expr-subst expr var
2096 (math-add (math-mul step var
)
2099 (math-simplify (math-div (math-sub high low
)
2101 (math-mul (if no-mul-flag
1 (math-pow -
1 low
))
2102 (calcFunc-sum (math-mul (math-pow -
1 var
) expr
) var low high
))))
2104 (defun math-sum-const-factors (expr var
)
2108 (while (eq (car-safe p
) '*)
2109 (if (math-expr-contains (nth 1 p
) var
)
2110 (setq not-const
(cons (nth 1 p
) not-const
))
2111 (setq const
(cons (nth 1 p
) const
)))
2113 (if (math-expr-contains p var
)
2114 (setq not-const
(cons p not-const
))
2115 (setq const
(cons p const
)))
2117 (cons (let ((temp (car const
)))
2118 (while (setq const
(cdr const
))
2119 (setq temp
(list '* (car const
) temp
)))
2121 (let ((temp (or (car not-const
) 1)))
2122 (while (setq not-const
(cdr not-const
))
2123 (setq temp
(list '* (car not-const
) temp
)))
2126 (defvar math-sum-int-pow-cache
(list '(0 1)))
2127 ;; Following is from CRC Math Tables, 27th ed, pp. 52-53.
2128 (defun math-sum-integer-power (pow)
2129 (let ((calc-prefer-frac t
)
2130 (n (length math-sum-int-pow-cache
)))
2132 (let* ((new (list 0 0))
2134 (pp (cdr (nth (1- n
) math-sum-int-pow-cache
)))
2139 (setq q
(math-div (car pp
) p
)
2140 new
(cons (math-mul q n
) new
)
2141 sum
(math-add sum q
)
2144 (setcar lin
(math-sub 1 (math-mul n sum
)))
2145 (setq math-sum-int-pow-cache
2146 (nconc math-sum-int-pow-cache
(list (nreverse new
)))
2148 (nth pow math-sum-int-pow-cache
)))
2150 (defun math-to-exponentials (expr)
2153 (let ((x (nth 1 expr
))
2154 (pi (if calc-symbolic-mode
'(var pi var-pi
) (math-pi)))
2155 (i (if calc-symbolic-mode
'(var i var-i
) '(cplx 0 1))))
2156 (cond ((eq (car expr
) 'calcFunc-exp
)
2157 (list '^
'(var e var-e
) x
))
2158 ((eq (car expr
) 'calcFunc-sin
)
2159 (or (eq calc-angle-mode
'rad
)
2160 (setq x
(list '/ (list '* x pi
) 180)))
2162 (list '^
'(var e var-e
) (list '* x i
))
2163 (list '^
'(var e var-e
)
2164 (list 'neg
(list '* x i
))))
2166 ((eq (car expr
) 'calcFunc-cos
)
2167 (or (eq calc-angle-mode
'rad
)
2168 (setq x
(list '/ (list '* x pi
) 180)))
2170 (list '^
'(var e var-e
)
2172 (list '^
'(var e var-e
)
2173 (list 'neg
(list '* x i
))))
2175 ((eq (car expr
) 'calcFunc-sinh
)
2177 (list '^
'(var e var-e
) x
)
2178 (list '^
'(var e var-e
) (list 'neg x
)))
2180 ((eq (car expr
) 'calcFunc-cosh
)
2182 (list '^
'(var e var-e
) x
)
2183 (list '^
'(var e var-e
) (list 'neg x
)))
2187 (defun math-to-exps (expr)
2188 (cond (calc-symbolic-mode expr
)
2190 (if (equal expr
'(var e var-e
)) (math-e) expr
))
2191 ((and (eq (car expr
) '^
)
2192 (equal (nth 1 expr
) '(var e var-e
)))
2193 (list 'calcFunc-exp
(nth 2 expr
)))
2195 (cons (car expr
) (mapcar 'math-to-exps
(cdr expr
))))))
2198 (defvar math-disable-prods nil
)
2199 (defun calcFunc-prod (expr var
&optional low high step
)
2200 (if math-disable-prods
(math-reject-arg))
2201 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec
2)))
2202 (math-prod-rec expr var low high step
)))
2203 (math-disable-prods t
))
2204 (math-normalize res
)))
2206 (defun math-prod-rec (expr var
&optional low high step
)
2207 (or low
(setq low
'(neg (var inf var-inf
)) high
'(var inf var-inf
)))
2208 (and low
(not high
) (setq high
'(var inf var-inf
)))
2212 ((not (math-expr-contains expr var
))
2213 (math-pow expr
(math-add (math-div (math-sub high low
) (or step
1))
2215 ((and step
(not (math-equal-int step
1)))
2216 (if (math-negp step
)
2217 (math-prod-rec expr var high low
(math-neg step
))
2218 (let ((lo (math-simplify (math-div low step
))))
2219 (if (math-known-num-integerp lo
)
2220 (math-prod-rec (math-normalize
2221 (math-expr-subst expr var
2222 (math-mul step var
)))
2223 var lo
(math-simplify (math-div high step
)))
2224 (math-prod-rec (math-normalize
2225 (math-expr-subst expr var
2226 (math-add (math-mul step
2230 (math-simplify (math-div (math-sub high low
)
2232 ((and (memq (car expr
) '(* /))
2233 (setq t1
(math-prod-rec (nth 1 expr
) var low high
)
2234 t2
(math-prod-rec (nth 2 expr
) var low high
))
2235 (not (and (math-expr-calls t1
'(calcFunc-prod))
2236 (math-expr-calls t2
'(calcFunc-prod)))))
2237 (list (car expr
) t1 t2
))
2238 ((and (eq (car expr
) '^
)
2239 (not (math-expr-contains (nth 2 expr
) var
)))
2240 (math-pow (math-prod-rec (nth 1 expr
) var low high
)
2242 ((and (eq (car expr
) '^
)
2243 (not (math-expr-contains (nth 1 expr
) var
)))
2244 (math-pow (nth 1 expr
)
2245 (calcFunc-sum (nth 2 expr
) var low high
)))
2246 ((eq (car expr
) 'sqrt
)
2247 (math-normalize (list 'calcFunc-sqrt
2248 (list 'calcFunc-prod
(nth 1 expr
)
2250 ((eq (car expr
) 'neg
)
2251 (math-mul (math-pow -
1 (math-add (math-sub high low
) 1))
2252 (math-prod-rec (nth 1 expr
) var low high
)))
2253 ((eq (car expr
) 'calcFunc-exp
)
2254 (list 'calcFunc-exp
(calcFunc-sum (nth 1 expr
) var low high
)))
2255 ((and (setq t1
(math-is-polynomial expr var
1))
2258 ((or (and (math-equal-int (nth 1 t1
) 1)
2259 (setq low
(math-simplify
2260 (math-add low
(car t1
)))
2262 (math-add high
(car t1
)))))
2263 (and (math-equal-int (nth 1 t1
) -
1)
2266 (math-sub (car t1
) high
))
2268 (math-sub (car t1
) t2
)))))
2269 (if (or (math-zerop low
) (math-zerop high
))
2271 (if (and (or (math-negp low
) (math-negp high
))
2272 (or (math-num-integerp low
)
2273 (math-num-integerp high
)))
2274 (if (math-posp high
)
2276 (math-mul (math-pow -
1
2278 (math-add low high
) 1))
2280 (list 'calcFunc-fact
2282 (list 'calcFunc-fact
2283 (math-sub -
1 high
)))))
2285 (list 'calcFunc-fact high
)
2286 (list 'calcFunc-fact
(math-sub low
1))))))
2287 ((and (or (and (math-equal-int (nth 1 t1
) 2)
2288 (setq t2
(math-simplify
2289 (math-add (math-mul low
2)
2292 (math-add (math-mul high
2)
2294 (and (math-equal-int (nth 1 t1
) -
2)
2295 (setq t2
(math-simplify
2302 (or (math-integerp t2
)
2303 (and (math-messy-integerp t2
)
2304 (setq t2
(math-trunc t2
)))
2306 (and (math-messy-integerp t3
)
2307 (setq t3
(math-trunc t3
)))))
2308 (if (or (math-zerop t2
) (math-zerop t3
))
2310 (if (or (math-evenp t2
) (math-evenp t3
))
2311 (if (or (math-negp t2
) (math-negp t3
))
2312 (if (math-posp high
)
2315 (list 'calcFunc-dfact
2317 (list 'calcFunc-dfact
2320 (list 'calcFunc-dfact t3
)
2321 (list 'calcFunc-dfact
2326 (list '/ (list '-
(list '- t2 t3
)
2330 (list 'calcFunc-dfact
2332 (list 'calcFunc-dfact
2336 (list 'calcFunc-dfact t3
)
2337 (list 'calcFunc-dfact
2341 (if (equal val
'(var nan var-nan
)) (setq val nil
))
2343 (let* ((math-tabulate-initial 1)
2344 (math-tabulate-function 'calcFunc-prod
))
2345 (calcFunc-table expr var low high
)))))
2350 (defvar math-solve-ranges nil
)
2351 (defvar math-solve-sign
)
2352 ;;; Attempt to reduce math-solve-lhs = math-solve-rhs to
2353 ;;; math-solve-var = math-solve-rhs', where math-solve-var appears
2354 ;;; in math-solve-lhs but not in math-solve-rhs or math-solve-rhs';
2355 ;;; return math-solve-rhs'.
2356 ;;; Uses global values: math-solve-var, math-solve-full.
2357 (defvar math-solve-var
)
2358 (defvar math-solve-full
)
2360 ;; The variables math-solve-lhs, math-solve-rhs and math-try-solve-sign
2361 ;; are local to math-try-solve-for, but are used by math-try-solve-prod.
2362 ;; (math-solve-lhs and math-solve-rhs are is also local to
2363 ;; math-decompose-poly, but used by math-solve-poly-funny-powers.)
2364 (defvar math-solve-lhs
)
2365 (defvar math-solve-rhs
)
2366 (defvar math-try-solve-sign
)
2368 (defun math-try-solve-for
2369 (math-solve-lhs math-solve-rhs
&optional math-try-solve-sign no-poly
)
2370 (let (math-t1 math-t2 math-t3
)
2371 (cond ((equal math-solve-lhs math-solve-var
)
2372 (setq math-solve-sign math-try-solve-sign
)
2373 (if (eq math-solve-full
'all
)
2374 (let ((vec (list 'vec
(math-evaluate-expr math-solve-rhs
)))
2376 (while math-solve-ranges
2377 (setq p
(car math-solve-ranges
)
2380 (while (setq p
(cdr p
))
2381 (setq newvec
(nconc newvec
2382 (cdr (math-expr-subst
2383 vec var
(car p
))))))
2385 math-solve-ranges
(cdr math-solve-ranges
)))
2386 (math-normalize vec
))
2388 ((Math-primp math-solve-lhs
)
2390 ((and (eq (car math-solve-lhs
) '-
)
2391 (eq (car-safe (nth 1 math-solve-lhs
)) (car-safe (nth 2 math-solve-lhs
)))
2392 (Math-zerop math-solve-rhs
)
2393 (= (length (nth 1 math-solve-lhs
)) 2)
2394 (= (length (nth 2 math-solve-lhs
)) 2)
2395 (setq math-t1
(get (car (nth 1 math-solve-lhs
)) 'math-inverse
))
2396 (setq math-t2
(funcall math-t1
'(var SOLVEDUM SOLVEDUM
)))
2397 (eq (math-expr-contains-count math-t2
'(var SOLVEDUM SOLVEDUM
)) 1)
2398 (setq math-t3
(math-solve-above-dummy math-t2
))
2399 (setq math-t1
(math-try-solve-for
2400 (math-sub (nth 1 (nth 1 math-solve-lhs
))
2403 (nth 1 (nth 2 math-solve-lhs
))))
2406 ((eq (car math-solve-lhs
) 'neg
)
2407 (math-try-solve-for (nth 1 math-solve-lhs
) (math-neg math-solve-rhs
)
2408 (and math-try-solve-sign
(- math-try-solve-sign
))))
2409 ((and (not (eq math-solve-full
't
)) (math-try-solve-prod)))
2412 (math-decompose-poly math-solve-lhs
2413 math-solve-var
15 math-solve-rhs
)))
2414 (setq math-t1
(cdr (nth 1 math-t2
))
2415 math-t1
(let ((math-solve-ranges math-solve-ranges
))
2416 (cond ((= (length math-t1
) 5)
2417 (apply 'math-solve-quartic
(car math-t2
) math-t1
))
2418 ((= (length math-t1
) 4)
2419 (apply 'math-solve-cubic
(car math-t2
) math-t1
))
2420 ((= (length math-t1
) 3)
2421 (apply 'math-solve-quadratic
(car math-t2
) math-t1
))
2422 ((= (length math-t1
) 2)
2423 (apply 'math-solve-linear
2424 (car math-t2
) math-try-solve-sign math-t1
))
2426 (math-poly-all-roots (car math-t2
) math-t1
))
2427 (calc-symbolic-mode nil
)
2431 (math-poly-any-root (reverse math-t1
) 0 t
)
2434 (if (eq (nth 2 math-t2
) 1)
2436 (math-solve-prod math-t1
(math-try-solve-for (nth 2 math-t2
) 0 nil t
)))
2437 (calc-record-why "*Unable to find a symbolic solution")
2439 ((and (math-solve-find-root-term math-solve-lhs nil
)
2440 (eq (math-expr-contains-count math-solve-lhs math-t1
) 1)) ; just in case
2441 (math-try-solve-for (math-simplify
2442 (math-sub (if (or math-t3
(math-evenp math-t2
))
2443 (math-pow math-t1 math-t2
)
2444 (math-neg (math-pow math-t1 math-t2
)))
2446 (math-sub (math-normalize
2448 math-solve-lhs math-t1
0))
2450 math-t2 math-solve-var
)))
2452 ((eq (car math-solve-lhs
) '+)
2453 (cond ((not (math-expr-contains (nth 1 math-solve-lhs
) math-solve-var
))
2454 (math-try-solve-for (nth 2 math-solve-lhs
)
2455 (math-sub math-solve-rhs
(nth 1 math-solve-lhs
))
2456 math-try-solve-sign
))
2457 ((not (math-expr-contains (nth 2 math-solve-lhs
) math-solve-var
))
2458 (math-try-solve-for (nth 1 math-solve-lhs
)
2459 (math-sub math-solve-rhs
(nth 2 math-solve-lhs
))
2460 math-try-solve-sign
))))
2461 ((eq (car math-solve-lhs
) 'calcFunc-eq
)
2462 (math-try-solve-for (math-sub (nth 1 math-solve-lhs
) (nth 2 math-solve-lhs
))
2463 math-solve-rhs math-try-solve-sign no-poly
))
2464 ((eq (car math-solve-lhs
) '-
)
2465 (cond ((or (and (eq (car-safe (nth 1 math-solve-lhs
)) 'calcFunc-sin
)
2466 (eq (car-safe (nth 2 math-solve-lhs
)) 'calcFunc-cos
))
2467 (and (eq (car-safe (nth 1 math-solve-lhs
)) 'calcFunc-cos
)
2468 (eq (car-safe (nth 2 math-solve-lhs
)) 'calcFunc-sin
)))
2469 (math-try-solve-for (math-sub (nth 1 math-solve-lhs
)
2470 (list (car (nth 1 math-solve-lhs
))
2472 (math-quarter-circle t
)
2473 (nth 1 (nth 2 math-solve-lhs
)))))
2475 ((not (math-expr-contains (nth 1 math-solve-lhs
) math-solve-var
))
2476 (math-try-solve-for (nth 2 math-solve-lhs
)
2477 (math-sub (nth 1 math-solve-lhs
) math-solve-rhs
)
2478 (and math-try-solve-sign
2479 (- math-try-solve-sign
))))
2480 ((not (math-expr-contains (nth 2 math-solve-lhs
) math-solve-var
))
2481 (math-try-solve-for (nth 1 math-solve-lhs
)
2482 (math-add math-solve-rhs
(nth 2 math-solve-lhs
))
2483 math-try-solve-sign
))))
2484 ((and (eq math-solve-full
't
) (math-try-solve-prod)))
2485 ((and (eq (car math-solve-lhs
) '%
)
2486 (not (math-expr-contains (nth 2 math-solve-lhs
) math-solve-var
)))
2487 (math-try-solve-for (nth 1 math-solve-lhs
) (math-add math-solve-rhs
2489 (nth 2 math-solve-lhs
)))))
2490 ((eq (car math-solve-lhs
) 'calcFunc-log
)
2491 (cond ((not (math-expr-contains (nth 2 math-solve-lhs
) math-solve-var
))
2492 (math-try-solve-for (nth 1 math-solve-lhs
)
2493 (math-pow (nth 2 math-solve-lhs
) math-solve-rhs
)))
2494 ((not (math-expr-contains (nth 1 math-solve-lhs
) math-solve-var
))
2495 (math-try-solve-for (nth 2 math-solve-lhs
) (math-pow
2496 (nth 1 math-solve-lhs
)
2497 (math-div 1 math-solve-rhs
))))))
2498 ((and (= (length math-solve-lhs
) 2)
2499 (symbolp (car math-solve-lhs
))
2500 (setq math-t1
(get (car math-solve-lhs
) 'math-inverse
))
2501 (setq math-t2
(funcall math-t1 math-solve-rhs
)))
2502 (setq math-t1
(get (car math-solve-lhs
) 'math-inverse-sign
))
2503 (math-try-solve-for (nth 1 math-solve-lhs
) (math-normalize math-t2
)
2504 (and math-try-solve-sign math-t1
2505 (if (integerp math-t1
)
2506 (* math-t1 math-try-solve-sign
)
2507 (funcall math-t1 math-solve-lhs
2508 math-try-solve-sign
)))))
2509 ((and (symbolp (car math-solve-lhs
))
2510 (setq math-t1
(get (car math-solve-lhs
) 'math-inverse-n
))
2511 (setq math-t2
(funcall math-t1 math-solve-lhs math-solve-rhs
)))
2513 ((setq math-t1
(math-expand-formula math-solve-lhs
))
2514 (math-try-solve-for math-t1 math-solve-rhs math-try-solve-sign
))
2516 (calc-record-why "*No inverse known" math-solve-lhs
)
2520 (defun math-try-solve-prod ()
2521 (cond ((eq (car math-solve-lhs
) '*)
2522 (cond ((not (math-expr-contains (nth 1 math-solve-lhs
) math-solve-var
))
2523 (math-try-solve-for (nth 2 math-solve-lhs
)
2524 (math-div math-solve-rhs
(nth 1 math-solve-lhs
))
2525 (math-solve-sign math-try-solve-sign
2526 (nth 1 math-solve-lhs
))))
2527 ((not (math-expr-contains (nth 2 math-solve-lhs
) math-solve-var
))
2528 (math-try-solve-for (nth 1 math-solve-lhs
)
2529 (math-div math-solve-rhs
(nth 2 math-solve-lhs
))
2530 (math-solve-sign math-try-solve-sign
2531 (nth 2 math-solve-lhs
))))
2532 ((Math-zerop math-solve-rhs
)
2533 (math-solve-prod (let ((math-solve-ranges math-solve-ranges
))
2534 (math-try-solve-for (nth 2 math-solve-lhs
) 0))
2535 (math-try-solve-for (nth 1 math-solve-lhs
) 0)))))
2536 ((eq (car math-solve-lhs
) '/)
2537 (cond ((not (math-expr-contains (nth 1 math-solve-lhs
) math-solve-var
))
2538 (math-try-solve-for (nth 2 math-solve-lhs
)
2539 (math-div (nth 1 math-solve-lhs
) math-solve-rhs
)
2540 (math-solve-sign math-try-solve-sign
2541 (nth 1 math-solve-lhs
))))
2542 ((not (math-expr-contains (nth 2 math-solve-lhs
) math-solve-var
))
2543 (math-try-solve-for (nth 1 math-solve-lhs
)
2544 (math-mul math-solve-rhs
(nth 2 math-solve-lhs
))
2545 (math-solve-sign math-try-solve-sign
2546 (nth 2 math-solve-lhs
))))
2547 ((setq math-t1
(math-try-solve-for (math-sub (nth 1 math-solve-lhs
)
2548 (math-mul (nth 2 math-solve-lhs
)
2552 ((eq (car math-solve-lhs
) '^
)
2553 (cond ((not (math-expr-contains (nth 1 math-solve-lhs
) math-solve-var
))
2555 (nth 2 math-solve-lhs
)
2556 (math-add (math-normalize
2557 (list 'calcFunc-log math-solve-rhs
(nth 1 math-solve-lhs
)))
2560 (math-mul '(var pi var-pi
)
2564 (list 'calcFunc-ln
(nth 1 math-solve-lhs
)))))))
2565 ((not (math-expr-contains (nth 2 math-solve-lhs
) math-solve-var
))
2566 (cond ((and (integerp (nth 2 math-solve-lhs
))
2567 (>= (nth 2 math-solve-lhs
) 2)
2568 (setq math-t1
(math-integer-log2 (nth 2 math-solve-lhs
))))
2569 (setq math-t2 math-solve-rhs
)
2570 (if (and (eq math-solve-full t
)
2571 (math-known-realp (nth 1 math-solve-lhs
)))
2573 (while (>= (setq math-t1
(1- math-t1
)) 0)
2574 (setq math-t2
(list 'calcFunc-sqrt math-t2
)))
2575 (setq math-t2
(math-solve-get-sign math-t2
)))
2576 (while (>= (setq math-t1
(1- math-t1
)) 0)
2577 (setq math-t2
(math-solve-get-sign
2579 (list 'calcFunc-sqrt math-t2
))))))
2581 (nth 1 math-solve-lhs
)
2582 (math-normalize math-t2
)))
2583 ((math-looks-negp (nth 2 math-solve-lhs
))
2585 (list '^
(nth 1 math-solve-lhs
)
2586 (math-neg (nth 2 math-solve-lhs
)))
2587 (math-div 1 math-solve-rhs
)))
2588 ((and (eq math-solve-full t
)
2589 (Math-integerp (nth 2 math-solve-lhs
))
2590 (math-known-realp (nth 1 math-solve-lhs
)))
2591 (setq math-t1
(math-normalize
2592 (list 'calcFunc-nroot math-solve-rhs
2593 (nth 2 math-solve-lhs
))))
2594 (if (math-evenp (nth 2 math-solve-lhs
))
2595 (setq math-t1
(math-solve-get-sign math-t1
)))
2597 (nth 1 math-solve-lhs
) math-t1
2598 (and math-try-solve-sign
2599 (math-oddp (nth 2 math-solve-lhs
))
2600 (math-solve-sign math-try-solve-sign
2601 (nth 2 math-solve-lhs
)))))
2602 (t (math-try-solve-for
2603 (nth 1 math-solve-lhs
)
2607 (if (Math-realp (nth 2 math-solve-lhs
))
2612 (and (integerp (nth 2 math-solve-lhs
))
2614 (nth 2 math-solve-lhs
)))))
2615 (math-div (nth 2 math-solve-lhs
) 2))
2622 (and (integerp (nth 2 math-solve-lhs
))
2624 (nth 2 math-solve-lhs
))))))
2625 (nth 2 math-solve-lhs
)))))
2627 (list 'calcFunc-nroot
2629 (nth 2 math-solve-lhs
))))
2630 (and math-try-solve-sign
2631 (math-oddp (nth 2 math-solve-lhs
))
2632 (math-solve-sign math-try-solve-sign
2633 (nth 2 math-solve-lhs
)))))))))
2636 (defun math-solve-prod (lsoln rsoln
)
2641 ((eq math-solve-full
'all
)
2642 (cons 'vec
(append (cdr lsoln
) (cdr rsoln
))))
2645 (list 'calcFunc-gt
(math-solve-get-sign 1) 0)
2650 ;;; This deals with negative, fractional, and symbolic powers of "x".
2651 ;; The variable math-solve-b is local to math-decompose-poly,
2652 ;; but is used by math-solve-poly-funny-powers.
2653 (defvar math-solve-b
)
2655 (defun math-solve-poly-funny-powers (sub-rhs) ; uses "t1", "t2"
2656 (setq math-t1 math-solve-lhs
)
2657 (let ((pp math-poly-neg-powers
)
2660 (setq fac
(math-pow (car pp
) (or math-poly-mult-powers
1))
2661 math-t1
(math-mul math-t1 fac
)
2662 math-solve-rhs
(math-mul math-solve-rhs fac
)
2664 (if sub-rhs
(setq math-t1
(math-sub math-t1 math-solve-rhs
)))
2665 (let ((math-poly-neg-powers nil
))
2666 (setq math-t2
(math-mul (or math-poly-mult-powers
1)
2667 (let ((calc-prefer-frac t
))
2668 (math-div 1 math-poly-frac-powers
)))
2669 math-t1
(math-is-polynomial
2670 (math-simplify (calcFunc-expand math-t1
)) math-solve-b
50))))
2672 ;;; This converts "a x^8 + b x^5 + c x^2" to "(a (x^3)^2 + b (x^3) + c) * x^2".
2673 (defun math-solve-crunch-poly (max-degree) ; uses "t1", "t3"
2675 (while (and math-t1
(Math-zerop (car math-t1
)))
2676 (setq math-t1
(cdr math-t1
)
2679 (let* ((degree (1- (length math-t1
)))
2681 (while (and (> scale
1) (= (car math-t3
) 1))
2682 (and (= (% degree scale
) 0)
2688 (if (= (% n scale
) 0)
2689 (setq new-t1
(nconc new-t1
(list (car p
))))
2690 (or (Math-zerop (car p
))
2695 (setq math-t3
(cons scale
(cdr math-t3
))
2697 (setq scale
(1- scale
)))
2698 (setq math-t3
(list (math-mul (car math-t3
) math-t2
)
2699 (math-mul count math-t2
)))
2700 (<= (1- (length math-t1
)) max-degree
)))))
2702 (defun calcFunc-poly (expr var
&optional degree
)
2704 (or (natnump degree
) (math-reject-arg degree
'fixnatnump
))
2706 (let ((p (math-is-polynomial expr var degree
'gen
)))
2711 (math-reject-arg expr
"Expected a polynomial"))))
2713 (defun calcFunc-gpoly (expr var
&optional degree
)
2715 (or (natnump degree
) (math-reject-arg degree
'fixnatnump
))
2717 (let* ((math-poly-base-variable var
)
2718 (d (math-decompose-poly expr var degree nil
)))
2721 (math-reject-arg expr
"Expected a polynomial"))))
2723 (defun math-decompose-poly (math-solve-lhs math-solve-var degree sub-rhs
)
2724 (let ((math-solve-rhs (or sub-rhs
1))
2725 math-t1 math-t2 math-t3
)
2726 (setq math-t2
(math-polynomial-base
2729 (lambda (math-solve-b)
2730 (let ((math-poly-neg-powers '(1))
2731 (math-poly-mult-powers nil
)
2732 (math-poly-frac-powers 1)
2733 (math-poly-exp-base t
))
2734 (and (not (equal math-solve-b math-solve-lhs
))
2735 (or (not (memq (car-safe math-solve-b
) '(+ -
))) sub-rhs
)
2736 (setq math-t3
'(1 0) math-t2
1
2737 math-t1
(math-is-polynomial math-solve-lhs
2739 (if (and (equal math-poly-neg-powers
'(1))
2740 (memq math-poly-mult-powers
'(nil 1))
2741 (eq math-poly-frac-powers
1)
2743 (setq math-t1
(cons (math-sub (car math-t1
) math-solve-rhs
)
2745 (math-solve-poly-funny-powers sub-rhs
))
2746 (math-solve-crunch-poly degree
)
2747 (or (math-expr-contains math-solve-b math-solve-var
)
2748 (math-expr-contains (car math-t3
) math-solve-var
))))))))
2750 (list (math-pow math-t2
(car math-t3
))
2753 (math-pow math-t2
(nth 1 math-t3
))
2754 (math-div (math-pow math-t2
(nth 1 math-t3
)) math-solve-rhs
))))))
2756 (defun math-solve-linear (var sign b a
)
2757 (math-try-solve-for var
2758 (math-div (math-neg b
) a
)
2759 (math-solve-sign sign a
)
2762 (defun math-solve-quadratic (var c b a
)
2765 (if (math-looks-evenp b
)
2766 (let ((halfb (math-div b
2)))
2770 (math-solve-get-sign
2772 (list 'calcFunc-sqrt
2773 (math-add (math-sqr halfb
)
2774 (math-mul (math-neg c
) a
))))))
2779 (math-solve-get-sign
2781 (list 'calcFunc-sqrt
2782 (math-add (math-sqr b
)
2783 (math-mul 4 (math-mul (math-neg c
) a
)))))))
2787 (defun math-solve-cubic (var d c b a
)
2788 (let* ((p (math-div b a
))
2792 (aa (math-sub q
(math-div psqr
3)))
2794 (math-div (math-sub (math-mul 2 (math-mul psqr p
))
2795 (math-mul 9 (math-mul p q
)))
2799 (math-try-solve-for (math-pow (math-add var
(math-div p
3)) 3)
2800 (math-neg bb
) nil t
)
2803 (math-mul (math-add var
(math-div p
3))
2804 (math-add (math-sqr (math-add var
(math-div p
3)))
2807 (setq m
(math-mul 2 (list 'calcFunc-sqrt
(math-div aa -
3))))
2816 (math-sub (list 'calcFunc-arccos
2817 (math-div (math-mul 3 bb
)
2821 (math-add 1 (math-solve-get-int
2824 calc-symbolic-mode
))))
2829 (defun math-solve-quartic (var d c b a aa
)
2830 (setq a
(math-div a aa
))
2831 (setq b
(math-div b aa
))
2832 (setq c
(math-div c aa
))
2833 (setq d
(math-div d aa
))
2836 (let* ((asqr (math-sqr a
))
2837 (asqr4 (math-div asqr
4))
2838 (y (let ((math-solve-full nil
)
2840 (math-solve-cubic math-solve-var
2842 (math-mul 4 (math-mul b d
))
2845 (math-sub (math-mul a c
)
2849 (rsqr (math-add (math-sub asqr4 b
) y
))
2850 (r (list 'calcFunc-sqrt rsqr
))
2851 (sign1 (math-solve-get-sign 1))
2852 (de (list 'calcFunc-sqrt
2854 (math-sub (math-mul 3 asqr4
)
2856 (if (Math-zerop rsqr
)
2860 (list 'calcFunc-sqrt
2861 (math-sub (math-sqr y
)
2867 (math-mul 4 (math-mul a b
))
2873 (math-sub (math-add (math-mul sign1
(math-div r
2))
2874 (math-solve-get-sign (math-div de
2)))
2878 (defvar math-symbolic-solve nil
)
2879 (defvar math-int-coefs nil
)
2881 ;; The variable math-int-threshold is local to math-poly-all-roots,
2882 ;; but is used by math-poly-newton-root.
2883 (defvar math-int-threshold
)
2884 ;; The variables math-int-scale, math-int-factors and math-double-roots
2885 ;; are local to math-poly-all-roots, but are used by math-poly-integer-root.
2886 (defvar math-int-scale
)
2887 (defvar math-int-factors
)
2888 (defvar math-double-roots
)
2890 (defun math-poly-all-roots (var p
&optional math-factoring
)
2892 (let* ((math-symbolic-solve calc-symbolic-mode
)
2894 (deg (1- (length p
)))
2895 (orig-p (reverse p
))
2896 (math-int-coefs nil
)
2897 (math-int-scale nil
)
2898 (math-double-roots nil
)
2899 (math-int-factors nil
)
2900 (math-int-threshold nil
)
2902 ;; If rational coefficients, look for exact rational factors.
2903 (while (and pp
(Math-ratp (car pp
)))
2906 (if (or math-factoring math-symbolic-solve
)
2908 (let ((lead (car orig-p
))
2909 (calc-prefer-frac t
)
2910 (scale (apply 'math-lcm-denoms p
)))
2911 (setq math-int-scale
(math-abs (math-mul scale lead
))
2912 math-int-threshold
(math-div '(float 5 -
2) math-int-scale
)
2913 math-int-coefs
(cdr (math-div (cons 'vec orig-p
) lead
)))))
2915 (let ((calc-prefer-frac nil
)
2916 (calc-symbolic-mode nil
)
2918 (def-p (copy-sequence orig-p
)))
2920 (if (Math-numberp (car pp
))
2923 (while (> deg
(if math-symbolic-solve
2 4))
2924 (let* ((x (math-poly-any-root def-p
'(float 0 0) nil
))
2926 (if (and (eq (car-safe x
) 'cplx
)
2927 (math-nearly-zerop (nth 2 x
) (nth 1 x
)))
2928 (setq x
(calcFunc-re x
)))
2930 (setq roots
(cons x roots
)))
2931 (or (math-numberp x
)
2932 (setq x
(math-evaluate-expr x
)))
2935 (while (setq pp
(cdr pp
))
2938 (setq b
(math-add (math-mul x b
) c
)))
2939 (setq def-p
(cdr def-p
)
2941 (setq p
(reverse def-p
))))
2943 (let ((math-solve-var '(var DUMMY var-DUMMY
))
2944 (math-solve-sign nil
)
2945 (math-solve-ranges nil
)
2946 (math-solve-full 'all
))
2947 (if (= (length p
) (length math-int-coefs
))
2948 (setq p
(reverse math-int-coefs
)))
2949 (setq roots
(append (cdr (apply (cond ((= deg
2)
2950 'math-solve-quadratic
)
2954 'math-solve-quartic
))
2958 (setq roots
(cons (math-div (math-neg (car p
)) (nth 1 p
))
2963 (math-poly-integer-root (car roots
))
2964 (setq roots
(cdr roots
)))
2965 (list math-int-factors
(nreverse math-int-coefs
) math-int-scale
))
2966 (let ((vec nil
) res
)
2968 (let ((root (car roots
))
2969 (math-solve-full (and math-solve-full
'all
)))
2970 (if (math-floatp root
)
2971 (setq root
(math-poly-any-root orig-p root t
)))
2972 (setq vec
(append vec
2973 (cdr (or (math-try-solve-for var root nil t
)
2974 (throw 'ouch nil
))))))
2975 (setq roots
(cdr roots
)))
2976 (setq vec
(cons 'vec
(nreverse vec
)))
2977 (if math-symbolic-solve
2978 (setq vec
(math-normalize vec
)))
2979 (if (eq math-solve-full t
)
2980 (list 'calcFunc-subscr
2982 (math-solve-get-int 1 (1- (length orig-p
)) 1))
2985 (defun math-lcm-denoms (&rest fracs
)
2988 (if (eq (car-safe (car fracs
)) 'frac
)
2989 (setq den
(calcFunc-lcm den
(nth 2 (car fracs
)))))
2990 (setq fracs
(cdr fracs
)))
2993 (defun math-poly-any-root (p x polish
) ; p is a reverse poly coeff list
2994 (let* ((newt (if (math-zerop x
)
2995 (math-poly-newton-root
2996 p
'(cplx (float 123 -
6) (float 1 -
4)) 4)
2997 (math-poly-newton-root p x
4)))
2998 (res (if (math-zerop (cdr newt
))
3000 (if (and (math-lessp (cdr newt
) '(float 1 -
3)) (not polish
))
3001 (setq newt
(math-poly-newton-root p
(car newt
) 30)))
3002 (if (math-zerop (cdr newt
))
3004 (math-poly-laguerre-root p x polish
)))))
3005 (and math-symbolic-solve
(math-floatp res
)
3009 (defun math-poly-newton-root (p x iters
)
3010 (let* ((calc-prefer-frac nil
)
3011 (calc-symbolic-mode nil
)
3012 (try-integer math-int-coefs
)
3014 (while (and (> (setq iters
(1- iters
)) 0)
3016 (math-working "newton" x
)
3019 (while (setq pp
(cdr pp
))
3020 (setq d
(math-add (math-mul x d
) b
)
3021 b
(math-add (math-mul x b
) (car pp
))))
3022 (not (math-zerop d
)))
3024 (setq dx
(math-div b d
)
3027 (let ((adx (math-abs-approx dx
)))
3028 (and (math-lessp adx math-int-threshold
)
3029 (let ((iroot (math-poly-integer-root x
)))
3032 (setq try-integer nil
))))))
3033 (or (not (or (eq dx
0)
3034 (math-nearly-zerop dx
(math-abs-approx x
))))
3035 (progn (setq dx
0) nil
)))))
3036 (cons x
(if (math-zerop x
)
3037 1 (math-div (math-abs-approx dx
) (math-abs-approx x
))))))
3039 (defun math-poly-integer-root (x)
3040 (and (math-lessp (calcFunc-xpon (math-abs-approx x
)) calc-internal-prec
)
3042 (let* ((calc-prefer-frac t
)
3043 (xre (calcFunc-re x
))
3044 (xim (calcFunc-im x
))
3045 (xresq (math-sqr xre
))
3046 (ximsq (math-sqr xim
)))
3047 (if (math-lessp ximsq
(calcFunc-scf xresq -
1))
3048 ;; Look for linear factor
3049 (let* ((rnd (math-div (math-round (math-mul xre math-int-scale
))
3051 (icp math-int-coefs
)
3054 (while (setq icp
(cdr icp
))
3055 (setq newcoef
(cons rem newcoef
)
3056 rem
(math-add (car icp
)
3057 (math-mul rem rnd
))))
3058 (and (math-zerop rem
)
3060 (setq math-int-coefs
(nreverse newcoef
)
3061 math-int-factors
(cons (list (math-neg rnd
))
3064 ;; Look for irreducible quadratic factor
3065 (let* ((rnd1 (math-div (math-round
3066 (math-mul xre
(math-mul -
2 math-int-scale
)))
3068 (sqscale (math-sqr math-int-scale
))
3069 (rnd0 (math-div (math-round (math-mul (math-add xresq ximsq
)
3072 (rem1 (car math-int-coefs
))
3073 (icp (cdr math-int-coefs
))
3076 (found (assoc (list rnd0 rnd1
(math-posp xim
))
3080 (setq math-double-roots
(delq found math-double-roots
)
3082 (while (setq icp
(cdr icp
))
3084 newcoef
(cons rem1 newcoef
)
3085 rem1
(math-sub rem0
(math-mul this rnd1
))
3086 rem0
(math-sub (car icp
) (math-mul this rnd0
)))))
3087 (and (math-zerop rem0
)
3089 (let ((aa (math-div rnd1 -
2)))
3090 (or found
(setq math-int-coefs
(reverse newcoef
)
3091 math-double-roots
(cons (list
3096 math-int-factors
(cons (cons rnd0 rnd1
)
3099 (let ((calc-symbolic-mode math-symbolic-solve
))
3100 (math-mul (math-sqrt (math-sub (math-sqr aa
)
3102 (if (math-negp xim
) -
1 1)))))))))))
3104 ;;; The following routine is from Numerical Recipes, section 9.5.
3105 (defun math-poly-laguerre-root (p x polish
)
3106 (let* ((calc-prefer-frac nil
)
3107 (calc-symbolic-mode nil
)
3110 (try-newt (not polish
))
3114 (and (or (< (setq iters
(1+ iters
)) 50)
3115 (math-reject-arg x
"*Laguerre's method failed to converge"))
3116 (let ((err (math-abs-approx (car p
)))
3117 (abx (math-abs-approx x
))
3121 (while (setq pp
(cdr pp
))
3122 (setq f
(math-add (math-mul x f
) d
)
3123 d
(math-add (math-mul x d
) b
)
3124 b
(math-add (math-mul x b
) (car pp
))
3125 err
(math-add (math-abs-approx b
) (math-mul abx err
))))
3126 (math-lessp (calcFunc-scf err
(- -
2 calc-internal-prec
))
3127 (math-abs-approx b
)))
3128 (or (not (math-zerop d
))
3129 (not (math-zerop f
))
3131 (setq x
(math-pow (math-neg b
) (list 'frac
1 m
)))
3133 (let* ((g (math-div d b
))
3135 (h (math-sub g2
(math-mul 2 (math-div f b
))))
3137 (math-mul (1- m
) (math-sub (math-mul m h
) g2
))))
3138 (gp (math-add g sq
))
3139 (gm (math-sub g sq
)))
3140 (if (math-lessp (calcFunc-abssqr gp
) (calcFunc-abssqr gm
))
3142 (setq dx
(math-div m gp
)
3145 (math-lessp (math-abs-approx dx
)
3146 (calcFunc-scf (math-abs-approx x
) -
3)))
3147 (let ((newt (math-poly-newton-root p x1
7)))
3150 (if (math-zerop (cdr newt
))
3151 (setq x
(car newt
) x1 x
)
3152 (if (math-lessp (cdr newt
) '(float 1 -
6))
3153 (let ((newt2 (math-poly-newton-root
3155 (if (math-zerop (cdr newt2
))
3156 (setq x
(car newt2
) x1 x
)
3157 (setq x
(car newt
))))))))
3159 (math-nearly-equal x x1
))))
3160 (let ((cdx (math-abs-approx dx
)))
3165 (math-lessp cdx dxold
)
3168 (let ((digs (calcFunc-xpon
3169 (math-div (math-abs-approx x
) cdx
))))
3171 "*Could not attain full precision")
3173 (let ((calc-internal-prec (max 3 digs
)))
3174 (setq x
(math-normalize x
))))))
3178 (math-lessp (calcFunc-scf (math-abs-approx x
)
3179 (- calc-internal-prec
))
3181 (or (and (math-floatp x
)
3182 (math-poly-integer-root x
))
3185 (defun math-solve-above-dummy (x)
3186 (and (not (Math-primp x
))
3187 (if (and (equal (nth 1 x
) '(var SOLVEDUM SOLVEDUM
))
3191 (while (and (setq x
(cdr x
))
3192 (not (setq res
(math-solve-above-dummy (car x
))))))
3195 (defun math-solve-find-root-term (x neg
) ; sets "t2", "t3"
3196 (if (math-solve-find-root-in-prod x
)
3199 (and (memq (car-safe x
) '(+ -
))
3200 (or (math-solve-find-root-term (nth 1 x
) neg
)
3201 (math-solve-find-root-term (nth 2 x
)
3202 (if (eq (car x
) '-
) (not neg
) neg
))))))
3204 (defun math-solve-find-root-in-prod (x)
3206 (math-expr-contains x math-solve-var
)
3207 (or (and (eq (car x
) 'calcFunc-sqrt
)
3209 (and (eq (car x
) '^
)
3210 (or (and (memq (math-quarter-integer (nth 2 x
)) '(1 2 3))
3212 (and (eq (car-safe (nth 2 x
)) 'frac
)
3213 (eq (nth 2 (nth 2 x
)) 3)
3215 (and (memq (car x
) '(* /))
3216 (or (and (not (math-expr-contains (nth 1 x
) math-solve-var
))
3217 (math-solve-find-root-in-prod (nth 2 x
)))
3218 (and (not (math-expr-contains (nth 2 x
) math-solve-var
))
3219 (math-solve-find-root-in-prod (nth 1 x
))))))))
3221 ;; The variable math-solve-vars is local to math-solve-system,
3222 ;; but is used by math-solve-system-rec.
3223 (defvar math-solve-vars
)
3225 ;; The variable math-solve-simplifying is local to math-solve-system
3226 ;; and math-solve-system-rec, but is used by math-solve-system-subst.
3227 (defvar math-solve-simplifying
)
3229 (defun math-solve-system (exprs math-solve-vars math-solve-full
)
3230 (setq exprs
(mapcar 'list
(if (Math-vectorp exprs
)
3233 math-solve-vars
(if (Math-vectorp math-solve-vars
)
3234 (cdr math-solve-vars
)
3235 (list math-solve-vars
)))
3236 (or (let ((math-solve-simplifying nil
))
3237 (math-solve-system-rec exprs math-solve-vars nil
))
3238 (let ((math-solve-simplifying t
))
3239 (math-solve-system-rec exprs math-solve-vars nil
))))
3241 ;;; The following backtracking solver works by choosing a variable
3242 ;;; and equation, and trying to solve the equation for the variable.
3243 ;;; If it succeeds it calls itself recursively with that variable and
3244 ;;; equation removed from their respective lists, and with the solution
3245 ;;; added to solns as well as being substituted into all existing
3246 ;;; equations. The algorithm terminates when any solution path
3247 ;;; manages to remove all the variables from var-list.
3249 ;;; To support calcFunc-roots, entries in eqn-list and solns are
3250 ;;; actually lists of equations.
3252 ;; The variables math-solve-system-res and math-solve-system-vv are
3253 ;; local to math-solve-system-rec, but are used by math-solve-system-subst.
3254 (defvar math-solve-system-vv
)
3255 (defvar math-solve-system-res
)
3258 (defun math-solve-system-rec (eqn-list var-list solns
)
3261 (math-solve-system-res nil
))
3263 ;; Try each variable in turn.
3267 (let* ((math-solve-system-vv (car v
))
3269 (elim (eq (car-safe math-solve-system-vv
) 'calcFunc-elim
)))
3271 (setq math-solve-system-vv
(nth 1 math-solve-system-vv
)))
3273 ;; Try each equation in turn.
3280 (setq math-solve-system-res nil
)
3282 ;; Try to solve for math-solve-system-vv the list of equations e2.
3284 (setq res2
(or (and (eq (car e2
) eprev
)
3286 (math-solve-for (car e2
) 0
3287 math-solve-system-vv
3289 (setq eprev
(car e2
)
3290 math-solve-system-res
(cons (if (eq math-solve-full
'all
)
3293 math-solve-system-res
)
3296 (setq math-solve-system-res nil
)
3298 ;; Found a solution. Now try other variables.
3299 (setq math-solve-system-res
(nreverse math-solve-system-res
)
3300 math-solve-system-res
(math-solve-system-rec
3302 'math-solve-system-subst
3304 (copy-sequence eqn-list
)))
3305 (delq (car v
) (copy-sequence var-list
))
3306 (let ((math-solve-simplifying nil
)
3312 (math-solve-system-subst
3318 math-solve-system-vv
3319 (apply 'append math-solve-system-res
))
3321 (not math-solve-system-res
))))
3323 (not math-solve-system-res
)))
3325 math-solve-system-res
)
3327 ;; Eliminated all variables, so now put solution into the proper format.
3328 (setq solns
(sort solns
3331 (not (memq (car x
) (memq (car y
) math-solve-vars
)))))))
3332 (if (eq math-solve-full
'all
)
3337 (mapcar (function (lambda (x) (cons 'vec
(cdr x
)))) solns
)
3338 (mapcar (function (lambda (x) (cons 'vec x
))) eqn-list
)))))
3342 (mapcar (function (lambda (x) (cons 'calcFunc-eq x
))) solns
)
3343 (mapcar 'car eqn-list
)))))))
3345 (defun math-solve-system-subst (x) ; uses "res" and "v"
3347 (res2 math-solve-system-res
))
3349 (setq accum
(nconc accum
3352 (if math-solve-simplifying
3355 (car x
) math-solve-system-vv r
))
3357 (car x
) math-solve-system-vv r
))))
3364 ;; calc-command-flags is declared in calc.el
3365 (defvar calc-command-flags
)
3367 (defun math-get-from-counter (name)
3368 (let ((ctr (assq name calc-command-flags
)))
3370 (setcdr ctr
(1+ (cdr ctr
)))
3371 (setq ctr
(cons name
1)
3372 calc-command-flags
(cons ctr calc-command-flags
)))
3375 (defvar var-GenCount
)
3377 (defun math-solve-get-sign (val)
3378 (setq val
(math-simplify val
))
3379 (if (and (eq (car-safe val
) '*)
3380 (Math-numberp (nth 1 val
)))
3381 (list '* (nth 1 val
) (math-solve-get-sign (nth 2 val
)))
3382 (and (eq (car-safe val
) 'calcFunc-sqrt
)
3383 (eq (car-safe (nth 1 val
)) '^
)
3384 (setq val
(math-normalize (list '^
3386 (math-div (nth 2 (nth 1 val
)) 2)))))
3388 (if (and (calc-var-value 'var-GenCount
)
3389 (Math-natnump var-GenCount
)
3390 (not (eq math-solve-full
'all
)))
3392 (math-mul (list 'calcFunc-as var-GenCount
) val
)
3393 (setq var-GenCount
(math-add var-GenCount
1))
3394 (calc-refresh-evaltos 'var-GenCount
))
3395 (let* ((var (concat "s" (int-to-string (math-get-from-counter 'solve-sign
))))
3396 (var2 (list 'var
(intern var
) (intern (concat "var-" var
)))))
3397 (if (eq math-solve-full
'all
)
3398 (setq math-solve-ranges
(cons (list var2
1 -
1)
3399 math-solve-ranges
)))
3400 (math-mul var2 val
)))
3401 (calc-record-why "*Choosing positive solution")
3404 (defun math-solve-get-int (val &optional range first
)
3406 (if (and (calc-var-value 'var-GenCount
)
3407 (Math-natnump var-GenCount
)
3408 (not (eq math-solve-full
'all
)))
3410 (math-mul val
(list 'calcFunc-an var-GenCount
))
3411 (setq var-GenCount
(math-add var-GenCount
1))
3412 (calc-refresh-evaltos 'var-GenCount
))
3413 (let* ((var (concat "n" (int-to-string
3414 (math-get-from-counter 'solve-int
))))
3415 (var2 (list 'var
(intern var
) (intern (concat "var-" var
)))))
3416 (if (and range
(eq math-solve-full
'all
))
3417 (setq math-solve-ranges
(cons (cons var2
3418 (cdr (calcFunc-index
3419 range
(or first
0))))
3420 math-solve-ranges
)))
3421 (math-mul val var2
)))
3422 (calc-record-why "*Choosing 0 for arbitrary integer in solution")
3425 (defun math-solve-sign (sign expr
)
3427 (let ((s1 (math-possible-signs expr
)))
3428 (cond ((memq s1
'(4 6))
3433 (defun math-looks-evenp (expr)
3434 (if (Math-integerp expr
)
3436 (if (memq (car expr
) '(* /))
3437 (math-looks-evenp (nth 1 expr
)))))
3439 (defun math-solve-for (lhs rhs math-solve-var math-solve-full
&optional sign
)
3440 (if (math-expr-contains rhs math-solve-var
)
3441 (math-solve-for (math-sub lhs rhs
) 0 math-solve-var math-solve-full
)
3442 (and (math-expr-contains lhs math-solve-var
)
3443 (math-with-extra-prec 1
3444 (let* ((math-poly-base-variable math-solve-var
)
3445 (res (math-try-solve-for lhs rhs sign
)))
3446 (if (and (eq math-solve-full
'all
)
3447 (math-known-realp math-solve-var
))
3448 (let ((old-len (length res
))
3453 (and (not (memq (car-safe x
)
3457 new-len
(length res
))
3458 (if (< new-len old-len
)
3459 (calc-record-why (if (= new-len
1)
3460 "*All solutions were complex"
3462 "*Omitted %d complex solutions"
3463 (- old-len new-len
)))))))
3466 (defun math-solve-eqn (expr var full
)
3467 (if (memq (car-safe expr
) '(calcFunc-neq calcFunc-lt calcFunc-gt
3468 calcFunc-leq calcFunc-geq
))
3469 (let ((res (math-solve-for (cons '-
(cdr expr
))
3471 (if (eq (car expr
) 'calcFunc-neq
) nil
1))))
3473 (if (eq math-solve-sign
1)
3474 (list (car expr
) var res
)
3475 (if (eq math-solve-sign -
1)
3476 (list (car expr
) res var
)
3477 (or (eq (car expr
) 'calcFunc-neq
)
3479 "*Can't determine direction of inequality"))
3480 (and (memq (car expr
) '(calcFunc-neq calcFunc-lt calcFunc-gt
))
3481 (list 'calcFunc-neq var res
))))))
3482 (let ((res (math-solve-for expr
0 var full
)))
3484 (list 'calcFunc-eq var res
)))))
3486 (defun math-reject-solution (expr var func
)
3487 (if (math-expr-contains expr var
)
3488 (or (equal (car calc-next-why
) '(* "Unable to find a symbolic solution"))
3489 (calc-record-why "*Unable to find a solution")))
3490 (list func expr var
))
3492 (defun calcFunc-solve (expr var
)
3493 (or (if (or (Math-vectorp expr
) (Math-vectorp var
))
3494 (math-solve-system expr var nil
)
3495 (math-solve-eqn expr var nil
))
3496 (math-reject-solution expr var
'calcFunc-solve
)))
3498 (defun calcFunc-fsolve (expr var
)
3499 (or (if (or (Math-vectorp expr
) (Math-vectorp var
))
3500 (math-solve-system expr var t
)
3501 (math-solve-eqn expr var t
))
3502 (math-reject-solution expr var
'calcFunc-fsolve
)))
3504 (defun calcFunc-roots (expr var
)
3505 (let ((math-solve-ranges nil
))
3506 (or (if (or (Math-vectorp expr
) (Math-vectorp var
))
3507 (math-solve-system expr var
'all
)
3508 (math-solve-for expr
0 var
'all
))
3509 (math-reject-solution expr var
'calcFunc-roots
))))
3511 (defun calcFunc-finv (expr var
)
3512 (let ((res (math-solve-for expr math-integ-var var nil
)))
3514 (math-normalize (math-expr-subst res math-integ-var var
))
3515 (math-reject-solution expr var
'calcFunc-finv
))))
3517 (defun calcFunc-ffinv (expr var
)
3518 (let ((res (math-solve-for expr math-integ-var var t
)))
3520 (math-normalize (math-expr-subst res math-integ-var var
))
3521 (math-reject-solution expr var
'calcFunc-finv
))))
3524 (put 'calcFunc-inv
'math-inverse
3525 (function (lambda (x) (math-div 1 x
))))
3526 (put 'calcFunc-inv
'math-inverse-sign -
1)
3528 (put 'calcFunc-sqrt
'math-inverse
3529 (function (lambda (x) (math-sqr x
))))
3531 (put 'calcFunc-conj
'math-inverse
3532 (function (lambda (x) (list 'calcFunc-conj x
))))
3534 (put 'calcFunc-abs
'math-inverse
3535 (function (lambda (x) (math-solve-get-sign x
))))
3537 (put 'calcFunc-deg
'math-inverse
3538 (function (lambda (x) (list 'calcFunc-rad x
))))
3539 (put 'calcFunc-deg
'math-inverse-sign
1)
3541 (put 'calcFunc-rad
'math-inverse
3542 (function (lambda (x) (list 'calcFunc-deg x
))))
3543 (put 'calcFunc-rad
'math-inverse-sign
1)
3545 (put 'calcFunc-ln
'math-inverse
3546 (function (lambda (x) (list 'calcFunc-exp x
))))
3547 (put 'calcFunc-ln
'math-inverse-sign
1)
3549 (put 'calcFunc-log10
'math-inverse
3550 (function (lambda (x) (list 'calcFunc-exp10 x
))))
3551 (put 'calcFunc-log10
'math-inverse-sign
1)
3553 (put 'calcFunc-lnp1
'math-inverse
3554 (function (lambda (x) (list 'calcFunc-expm1 x
))))
3555 (put 'calcFunc-lnp1
'math-inverse-sign
1)
3557 (put 'calcFunc-exp
'math-inverse
3558 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-ln x
))
3560 (math-mul '(var pi var-pi
)
3562 '(var i var-i
))))))))
3563 (put 'calcFunc-exp
'math-inverse-sign
1)
3565 (put 'calcFunc-expm1
'math-inverse
3566 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-lnp1 x
))
3568 (math-mul '(var pi var-pi
)
3570 '(var i var-i
))))))))
3571 (put 'calcFunc-expm1
'math-inverse-sign
1)
3573 (put 'calcFunc-sin
'math-inverse
3574 (function (lambda (x) (let ((n (math-solve-get-int 1)))
3575 (math-add (math-mul (math-normalize
3576 (list 'calcFunc-arcsin x
))
3578 (math-mul (math-half-circle t
)
3581 (put 'calcFunc-cos
'math-inverse
3582 (function (lambda (x) (math-add (math-solve-get-sign
3584 (list 'calcFunc-arccos x
)))
3586 (math-full-circle t
))))))
3588 (put 'calcFunc-tan
'math-inverse
3589 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-arctan x
))
3591 (math-half-circle t
))))))
3593 (put 'calcFunc-arcsin
'math-inverse
3594 (function (lambda (x) (math-normalize (list 'calcFunc-sin x
)))))
3596 (put 'calcFunc-arccos
'math-inverse
3597 (function (lambda (x) (math-normalize (list 'calcFunc-cos x
)))))
3599 (put 'calcFunc-arctan
'math-inverse
3600 (function (lambda (x) (math-normalize (list 'calcFunc-tan x
)))))
3602 (put 'calcFunc-sinh
'math-inverse
3603 (function (lambda (x) (let ((n (math-solve-get-int 1)))
3604 (math-add (math-mul (math-normalize
3605 (list 'calcFunc-arcsinh x
))
3607 (math-mul (math-half-circle t
)
3611 (put 'calcFunc-sinh
'math-inverse-sign
1)
3613 (put 'calcFunc-cosh
'math-inverse
3614 (function (lambda (x) (math-add (math-solve-get-sign
3616 (list 'calcFunc-arccosh x
)))
3617 (math-mul (math-full-circle t
)
3619 '(var i var-i
)))))))
3621 (put 'calcFunc-tanh
'math-inverse
3622 (function (lambda (x) (math-add (math-normalize
3623 (list 'calcFunc-arctanh x
))
3624 (math-mul (math-half-circle t
)
3626 '(var i var-i
)))))))
3627 (put 'calcFunc-tanh
'math-inverse-sign
1)
3629 (put 'calcFunc-arcsinh
'math-inverse
3630 (function (lambda (x) (math-normalize (list 'calcFunc-sinh x
)))))
3631 (put 'calcFunc-arcsinh
'math-inverse-sign
1)
3633 (put 'calcFunc-arccosh
'math-inverse
3634 (function (lambda (x) (math-normalize (list 'calcFunc-cosh x
)))))
3636 (put 'calcFunc-arctanh
'math-inverse
3637 (function (lambda (x) (math-normalize (list 'calcFunc-tanh x
)))))
3638 (put 'calcFunc-arctanh
'math-inverse-sign
1)
3642 (defun calcFunc-taylor (expr var num
)
3643 (let ((x0 0) (v var
))
3644 (if (memq (car-safe var
) '(+ - calcFunc-eq
))
3645 (setq x0
(if (eq (car var
) '+) (math-neg (nth 2 var
)) (nth 2 var
))
3647 (or (and (eq (car-safe v
) 'var
)
3648 (math-expr-contains expr v
)
3650 (let ((accum (math-expr-subst expr v x0
))
3651 (var2 (if (eq (car var
) 'calcFunc-eq
)
3657 (while (and (<= (setq n
(1+ n
)) num
)
3658 (setq fprime
(calcFunc-deriv fprime v nil t
)))
3659 (setq fprime
(math-simplify fprime
)
3660 nfac
(math-mul nfac n
)
3661 accum
(math-add accum
3662 (math-div (math-mul (math-pow var2 n
)
3667 (math-normalize accum
))))
3668 (list 'calcFunc-taylor expr var num
))))
3672 ;; arch-tag: f2932ec8-dd63-418b-a542-11a644b9d4c4
3673 ;;; calcalg2.el ends here