Merge from origin/emacs-24
[emacs.git] / lisp / calc / calcalg2.el
blob1b86f79149a8e8b19ca1a37a6196366401de24f4
1 ;;; calcalg2.el --- more algebraic functions for Calc
3 ;; Copyright (C) 1990-1993, 2001-2014 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
6 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
8 ;; This file is part of GNU Emacs.
10 ;; GNU Emacs is free software: you can redistribute it and/or modify
11 ;; it under the terms of the GNU General Public License as published by
12 ;; the Free Software Foundation, either version 3 of the License, or
13 ;; (at your option) any later version.
15 ;; GNU Emacs is distributed in the hope that it will be useful,
16 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 ;; GNU General Public License for more details.
20 ;; You should have received a copy of the GNU General Public License
21 ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
23 ;;; Commentary:
25 ;;; Code:
27 ;; This file is autoloaded from calc-ext.el.
29 (require 'calc-ext)
30 (require 'calc-macs)
32 (defun calc-derivative (var num)
33 (interactive "sDifferentiate with respect to: \np")
34 (calc-slow-wrapper
35 (when (< num 0)
36 (error "Order of derivative must be positive"))
37 (let ((func (if (calc-is-hyperbolic) 'calcFunc-tderiv 'calcFunc-deriv))
38 n expr)
39 (if (or (equal var "") (equal var "$"))
40 (setq n 2
41 expr (calc-top-n 2)
42 var (calc-top-n 1))
43 (setq var (math-read-expr var))
44 (when (eq (car-safe var) 'error)
45 (error "Bad format in expression: %s" (nth 1 var)))
46 (setq n 1
47 expr (calc-top-n 1)))
48 (while (>= (setq num (1- num)) 0)
49 (setq expr (list func expr var)))
50 (calc-enter-result n "derv" expr))))
52 (defun calc-integral (var &optional arg)
53 (interactive "sIntegration variable: \nP")
54 (if arg
55 (calc-tabular-command 'calcFunc-integ "Integration" "intg" nil var nil nil)
56 (calc-slow-wrapper
57 (if (or (equal var "") (equal var "$"))
58 (calc-enter-result 2 "intg" (list 'calcFunc-integ
59 (calc-top-n 2)
60 (calc-top-n 1)))
61 (let ((var (math-read-expr var)))
62 (if (eq (car-safe var) 'error)
63 (error "Bad format in expression: %s" (nth 1 var)))
64 (calc-enter-result 1 "intg" (list 'calcFunc-integ
65 (calc-top-n 1)
66 var)))))))
68 (defun calc-num-integral (&optional varname lowname highname)
69 (interactive "sIntegration variable: ")
70 (calc-tabular-command 'calcFunc-ninteg "Integration" "nint"
71 nil varname lowname highname))
73 (defun calc-summation (arg &optional varname lowname highname)
74 (interactive "P\nsSummation variable: ")
75 (calc-tabular-command 'calcFunc-sum "Summation" "sum"
76 arg varname lowname highname))
78 (defun calc-alt-summation (arg &optional varname lowname highname)
79 (interactive "P\nsSummation variable: ")
80 (calc-tabular-command 'calcFunc-asum "Summation" "asum"
81 arg varname lowname highname))
83 (defun calc-product (arg &optional varname lowname highname)
84 (interactive "P\nsIndex variable: ")
85 (calc-tabular-command 'calcFunc-prod "Index" "prod"
86 arg varname lowname highname))
88 (defun calc-tabulate (arg &optional varname lowname highname)
89 (interactive "P\nsIndex variable: ")
90 (calc-tabular-command 'calcFunc-table "Index" "tabl"
91 arg varname lowname highname))
93 (defun calc-tabular-command (func prompt prefix arg varname lowname highname)
94 (calc-slow-wrapper
95 (let (var (low nil) (high nil) (step nil) stepname stepnum (num 1) expr)
96 (if (consp arg)
97 (setq stepnum 1)
98 (setq stepnum 0))
99 (if (or (equal varname "") (equal varname "$") (null varname))
100 (setq high (calc-top-n (+ stepnum 1))
101 low (calc-top-n (+ stepnum 2))
102 var (calc-top-n (+ stepnum 3))
103 num (+ stepnum 4))
104 (setq var (if (stringp varname) (math-read-expr varname) varname))
105 (if (eq (car-safe var) 'error)
106 (error "Bad format in expression: %s" (nth 1 var)))
107 (or lowname
108 (setq lowname (read-string (concat prompt " variable: " varname
109 ", from: "))))
110 (if (or (equal lowname "") (equal lowname "$"))
111 (setq high (calc-top-n (+ stepnum 1))
112 low (calc-top-n (+ stepnum 2))
113 num (+ stepnum 3))
114 (setq low (if (stringp lowname) (math-read-expr lowname) lowname))
115 (if (eq (car-safe low) 'error)
116 (error "Bad format in expression: %s" (nth 1 low)))
117 (or highname
118 (setq highname (read-string (concat prompt " variable: " varname
119 ", from: " lowname
120 ", to: "))))
121 (if (or (equal highname "") (equal highname "$"))
122 (setq high (calc-top-n (+ stepnum 1))
123 num (+ stepnum 2))
124 (setq high (if (stringp highname) (math-read-expr highname)
125 highname))
126 (if (eq (car-safe high) 'error)
127 (error "Bad format in expression: %s" (nth 1 high)))
128 (if (consp arg)
129 (progn
130 (setq stepname (read-string (concat prompt " variable: "
131 varname
132 ", from: " lowname
133 ", to: " highname
134 ", step: ")))
135 (if (or (equal stepname "") (equal stepname "$"))
136 (setq step (calc-top-n 1)
137 num 2)
138 (setq step (math-read-expr stepname))
139 (if (eq (car-safe step) 'error)
140 (error "Bad format in expression: %s"
141 (nth 1 step)))))))))
142 (or step
143 (if (consp arg)
144 (setq step (calc-top-n 1))
145 (if arg
146 (setq step (prefix-numeric-value arg)))))
147 (setq expr (calc-top-n num))
148 (calc-enter-result num prefix (append (list func expr var low high)
149 (and step (list step)))))))
151 (defun calc-solve-for (var)
152 (interactive "sVariable(s) to solve for: ")
153 (calc-slow-wrapper
154 (let ((func (if (calc-is-inverse)
155 (if (calc-is-hyperbolic) 'calcFunc-ffinv 'calcFunc-finv)
156 (if (calc-is-hyperbolic) 'calcFunc-fsolve 'calcFunc-solve))))
157 (if (or (equal var "") (equal var "$"))
158 (calc-enter-result 2 "solv" (list func
159 (calc-top-n 2)
160 (calc-top-n 1)))
161 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
162 (not (string-match "\\[" var)))
163 (math-read-expr (concat "[" var "]"))
164 (math-read-expr var))))
165 (if (eq (car-safe var) 'error)
166 (error "Bad format in expression: %s" (nth 1 var)))
167 (calc-enter-result 1 "solv" (list func
168 (calc-top-n 1)
169 var)))))))
171 (defun calc-poly-roots (var)
172 (interactive "sVariable to solve for: ")
173 (calc-slow-wrapper
174 (if (or (equal var "") (equal var "$"))
175 (calc-enter-result 2 "prts" (list 'calcFunc-roots
176 (calc-top-n 2)
177 (calc-top-n 1)))
178 (let ((var (if (and (string-match ",\\|[^ ] +[^ ]" var)
179 (not (string-match "\\[" var)))
180 (math-read-expr (concat "[" var "]"))
181 (math-read-expr var))))
182 (if (eq (car-safe var) 'error)
183 (error "Bad format in expression: %s" (nth 1 var)))
184 (calc-enter-result 1 "prts" (list 'calcFunc-roots
185 (calc-top-n 1)
186 var))))))
188 (defun calc-taylor (var nterms)
189 (interactive "sTaylor expansion variable: \nNNumber of terms: ")
190 (calc-slow-wrapper
191 (let ((var (math-read-expr var)))
192 (if (eq (car-safe var) 'error)
193 (error "Bad format in expression: %s" (nth 1 var)))
194 (calc-enter-result 1 "tylr" (list 'calcFunc-taylor
195 (calc-top-n 1)
197 (prefix-numeric-value nterms))))))
200 ;; The following are global variables used by math-derivative and some
201 ;; related functions
202 (defvar math-deriv-var)
203 (defvar math-deriv-total)
204 (defvar math-deriv-symb)
205 (defvar math-decls-cache)
206 (defvar math-decls-all)
208 (defun math-derivative (expr)
209 (cond ((equal expr math-deriv-var)
211 ((or (Math-scalarp expr)
212 (eq (car expr) 'sdev)
213 (and (eq (car expr) 'var)
214 (or (not math-deriv-total)
215 (math-const-var expr)
216 (progn
217 (math-setup-declarations)
218 (memq 'const (nth 1 (or (assq (nth 2 expr)
219 math-decls-cache)
220 math-decls-all)))))))
222 ((eq (car expr) '+)
223 (math-add (math-derivative (nth 1 expr))
224 (math-derivative (nth 2 expr))))
225 ((eq (car expr) '-)
226 (math-sub (math-derivative (nth 1 expr))
227 (math-derivative (nth 2 expr))))
228 ((memq (car expr) '(calcFunc-eq calcFunc-neq calcFunc-lt
229 calcFunc-gt calcFunc-leq calcFunc-geq))
230 (list (car expr)
231 (math-derivative (nth 1 expr))
232 (math-derivative (nth 2 expr))))
233 ((eq (car expr) 'neg)
234 (math-neg (math-derivative (nth 1 expr))))
235 ((eq (car expr) '*)
236 (math-add (math-mul (nth 2 expr)
237 (math-derivative (nth 1 expr)))
238 (math-mul (nth 1 expr)
239 (math-derivative (nth 2 expr)))))
240 ((eq (car expr) '/)
241 (math-sub (math-div (math-derivative (nth 1 expr))
242 (nth 2 expr))
243 (math-div (math-mul (nth 1 expr)
244 (math-derivative (nth 2 expr)))
245 (math-sqr (nth 2 expr)))))
246 ((eq (car expr) '^)
247 (let ((du (math-derivative (nth 1 expr)))
248 (dv (math-derivative (nth 2 expr))))
249 (or (Math-zerop du)
250 (setq du (math-mul (nth 2 expr)
251 (math-mul (math-normalize
252 (list '^
253 (nth 1 expr)
254 (math-add (nth 2 expr) -1)))
255 du))))
256 (or (Math-zerop dv)
257 (setq dv (math-mul (math-normalize
258 (list 'calcFunc-ln (nth 1 expr)))
259 (math-mul expr dv))))
260 (math-add du dv)))
261 ((eq (car expr) '%)
262 (math-derivative (nth 1 expr))) ; a reasonable definition
263 ((eq (car expr) 'vec)
264 (math-map-vec 'math-derivative expr))
265 ((and (memq (car expr) '(calcFunc-conj calcFunc-re calcFunc-im))
266 (= (length expr) 2))
267 (list (car expr) (math-derivative (nth 1 expr))))
268 ((and (memq (car expr) '(calcFunc-subscr calcFunc-mrow calcFunc-mcol))
269 (= (length expr) 3))
270 (let ((d (math-derivative (nth 1 expr))))
271 (if (math-numberp d)
272 0 ; assume x and x_1 are independent vars
273 (list (car expr) d (nth 2 expr)))))
274 (t (or (and (symbolp (car expr))
275 (if (= (length expr) 2)
276 (let ((handler (get (car expr) 'math-derivative)))
277 (and handler
278 (let ((deriv (math-derivative (nth 1 expr))))
279 (if (Math-zerop deriv)
280 deriv
281 (math-mul (funcall handler (nth 1 expr))
282 deriv)))))
283 (let ((handler (get (car expr) 'math-derivative-n)))
284 (and handler
285 (funcall handler expr)))))
286 (and (not (eq math-deriv-symb 'pre-expand))
287 (let ((exp (math-expand-formula expr)))
288 (and exp
289 (or (let ((math-deriv-symb 'pre-expand))
290 (catch 'math-deriv (math-derivative expr)))
291 (math-derivative exp)))))
292 (if (or (Math-objvecp expr)
293 (eq (car expr) 'var)
294 (not (symbolp (car expr))))
295 (if math-deriv-symb
296 (throw 'math-deriv nil)
297 (list (if math-deriv-total 'calcFunc-tderiv 'calcFunc-deriv)
298 expr
299 math-deriv-var))
300 (let ((accum 0)
301 (arg expr)
302 (n 1)
303 derv)
304 (while (setq arg (cdr arg))
305 (or (Math-zerop (setq derv (math-derivative (car arg))))
306 (let ((func (intern (concat (symbol-name (car expr))
308 (if (> n 1)
309 (int-to-string n)
310 ""))))
311 (prop (cond ((= (length expr) 2)
312 'math-derivative-1)
313 ((= (length expr) 3)
314 'math-derivative-2)
315 ((= (length expr) 4)
316 'math-derivative-3)
317 ((= (length expr) 5)
318 'math-derivative-4)
319 ((= (length expr) 6)
320 'math-derivative-5))))
321 (setq accum
322 (math-add
323 accum
324 (math-mul
325 derv
326 (let ((handler (get func prop)))
327 (or (and prop handler
328 (apply handler (cdr expr)))
329 (if (and math-deriv-symb
330 (not (get func
331 'calc-user-defn)))
332 (throw 'math-deriv nil)
333 (cons func (cdr expr))))))))))
334 (setq n (1+ n)))
335 accum))))))
337 (defun calcFunc-deriv (expr math-deriv-var &optional deriv-value math-deriv-symb)
338 (let* ((math-deriv-total nil)
339 (res (catch 'math-deriv (math-derivative expr))))
340 (or (eq (car-safe res) 'calcFunc-deriv)
341 (null res)
342 (setq res (math-normalize res)))
343 (and res
344 (if deriv-value
345 (math-expr-subst res math-deriv-var deriv-value)
346 res))))
348 (defun calcFunc-tderiv (expr math-deriv-var &optional deriv-value math-deriv-symb)
349 (math-setup-declarations)
350 (let* ((math-deriv-total t)
351 (res (catch 'math-deriv (math-derivative expr))))
352 (or (eq (car-safe res) 'calcFunc-tderiv)
353 (null res)
354 (setq res (math-normalize res)))
355 (and res
356 (if deriv-value
357 (math-expr-subst res math-deriv-var deriv-value)
358 res))))
360 (put 'calcFunc-inv\' 'math-derivative-1
361 (function (lambda (u) (math-neg (math-div 1 (math-sqr u))))))
363 (put 'calcFunc-sqrt\' 'math-derivative-1
364 (function (lambda (u) (math-div 1 (math-mul 2 (list 'calcFunc-sqrt u))))))
366 (put 'calcFunc-deg\' 'math-derivative-1
367 (function (lambda (u) (math-div-float '(float 18 1) (math-pi)))))
369 (put 'calcFunc-rad\' 'math-derivative-1
370 (function (lambda (u) (math-pi-over-180))))
372 (put 'calcFunc-ln\' 'math-derivative-1
373 (function (lambda (u) (math-div 1 u))))
375 (put 'calcFunc-log10\' 'math-derivative-1
376 (function (lambda (u)
377 (math-div (math-div 1 (math-normalize '(calcFunc-ln 10)))
378 u))))
380 (put 'calcFunc-lnp1\' 'math-derivative-1
381 (function (lambda (u) (math-div 1 (math-add u 1)))))
383 (put 'calcFunc-log\' 'math-derivative-2
384 (function (lambda (x b)
385 (and (not (Math-zerop b))
386 (let ((lnv (math-normalize
387 (list 'calcFunc-ln b))))
388 (math-div 1 (math-mul lnv x)))))))
390 (put 'calcFunc-log\'2 'math-derivative-2
391 (function (lambda (x b)
392 (let ((lnv (list 'calcFunc-ln b)))
393 (math-neg (math-div (list 'calcFunc-log x b)
394 (math-mul lnv b)))))))
396 (put 'calcFunc-exp\' 'math-derivative-1
397 (function (lambda (u) (math-normalize (list 'calcFunc-exp u)))))
399 (put 'calcFunc-expm1\' 'math-derivative-1
400 (function (lambda (u) (math-normalize (list 'calcFunc-expm1 u)))))
402 (put 'calcFunc-sin\' 'math-derivative-1
403 (function (lambda (u) (math-to-radians-2 (math-normalize
404 (list 'calcFunc-cos u)) t))))
406 (put 'calcFunc-cos\' 'math-derivative-1
407 (function (lambda (u) (math-neg (math-to-radians-2
408 (math-normalize
409 (list 'calcFunc-sin u)) t)))))
411 (put 'calcFunc-tan\' 'math-derivative-1
412 (function (lambda (u) (math-to-radians-2
413 (math-sqr
414 (math-normalize
415 (list 'calcFunc-sec u))) t))))
417 (put 'calcFunc-sec\' 'math-derivative-1
418 (function (lambda (u) (math-to-radians-2
419 (math-mul
420 (math-normalize
421 (list 'calcFunc-sec u))
422 (math-normalize
423 (list 'calcFunc-tan u))) t))))
425 (put 'calcFunc-csc\' 'math-derivative-1
426 (function (lambda (u) (math-neg
427 (math-to-radians-2
428 (math-mul
429 (math-normalize
430 (list 'calcFunc-csc u))
431 (math-normalize
432 (list 'calcFunc-cot u))) t)))))
434 (put 'calcFunc-cot\' 'math-derivative-1
435 (function (lambda (u) (math-neg
436 (math-to-radians-2
437 (math-sqr
438 (math-normalize
439 (list 'calcFunc-csc u))) t)))))
441 (put 'calcFunc-arcsin\' 'math-derivative-1
442 (function (lambda (u)
443 (math-from-radians-2
444 (math-div 1 (math-normalize
445 (list 'calcFunc-sqrt
446 (math-sub 1 (math-sqr u))))) t))))
448 (put 'calcFunc-arccos\' 'math-derivative-1
449 (function (lambda (u)
450 (math-from-radians-2
451 (math-div -1 (math-normalize
452 (list 'calcFunc-sqrt
453 (math-sub 1 (math-sqr u))))) t))))
455 (put 'calcFunc-arctan\' 'math-derivative-1
456 (function (lambda (u) (math-from-radians-2
457 (math-div 1 (math-add 1 (math-sqr u))) t))))
459 (put 'calcFunc-sinh\' 'math-derivative-1
460 (function (lambda (u) (math-normalize (list 'calcFunc-cosh u)))))
462 (put 'calcFunc-cosh\' 'math-derivative-1
463 (function (lambda (u) (math-normalize (list 'calcFunc-sinh u)))))
465 (put 'calcFunc-tanh\' 'math-derivative-1
466 (function (lambda (u) (math-sqr
467 (math-normalize
468 (list 'calcFunc-sech u))))))
470 (put 'calcFunc-sech\' 'math-derivative-1
471 (function (lambda (u) (math-neg
472 (math-mul
473 (math-normalize (list 'calcFunc-sech u))
474 (math-normalize (list 'calcFunc-tanh u)))))))
476 (put 'calcFunc-csch\' 'math-derivative-1
477 (function (lambda (u) (math-neg
478 (math-mul
479 (math-normalize (list 'calcFunc-csch u))
480 (math-normalize (list 'calcFunc-coth u)))))))
482 (put 'calcFunc-coth\' 'math-derivative-1
483 (function (lambda (u) (math-neg
484 (math-sqr
485 (math-normalize
486 (list 'calcFunc-csch u)))))))
488 (put 'calcFunc-arcsinh\' 'math-derivative-1
489 (function (lambda (u)
490 (math-div 1 (math-normalize
491 (list 'calcFunc-sqrt
492 (math-add (math-sqr u) 1)))))))
494 (put 'calcFunc-arccosh\' 'math-derivative-1
495 (function (lambda (u)
496 (math-div 1 (math-normalize
497 (list 'calcFunc-sqrt
498 (math-add (math-sqr u) -1)))))))
500 (put 'calcFunc-arctanh\' 'math-derivative-1
501 (function (lambda (u) (math-div 1 (math-sub 1 (math-sqr u))))))
503 (put 'calcFunc-bern\'2 'math-derivative-2
504 (function (lambda (n x)
505 (math-mul n (list 'calcFunc-bern (math-add n -1) x)))))
507 (put 'calcFunc-euler\'2 'math-derivative-2
508 (function (lambda (n x)
509 (math-mul n (list 'calcFunc-euler (math-add n -1) x)))))
511 (put 'calcFunc-gammag\'2 'math-derivative-2
512 (function (lambda (a x) (math-deriv-gamma a x 1))))
514 (put 'calcFunc-gammaG\'2 'math-derivative-2
515 (function (lambda (a x) (math-deriv-gamma a x -1))))
517 (put 'calcFunc-gammaP\'2 'math-derivative-2
518 (function (lambda (a x) (math-deriv-gamma a x
519 (math-div
520 1 (math-normalize
521 (list 'calcFunc-gamma
522 a)))))))
524 (put 'calcFunc-gammaQ\'2 'math-derivative-2
525 (function (lambda (a x) (math-deriv-gamma a x
526 (math-div
527 -1 (math-normalize
528 (list 'calcFunc-gamma
529 a)))))))
531 (defun math-deriv-gamma (a x scale)
532 (math-mul scale
533 (math-mul (math-pow x (math-add a -1))
534 (list 'calcFunc-exp (math-neg x)))))
536 (put 'calcFunc-betaB\' 'math-derivative-3
537 (function (lambda (x a b) (math-deriv-beta x a b 1))))
539 (put 'calcFunc-betaI\' 'math-derivative-3
540 (function (lambda (x a b) (math-deriv-beta x a b
541 (math-div
542 1 (list 'calcFunc-beta
543 a b))))))
545 (defun math-deriv-beta (x a b scale)
546 (math-mul (math-mul (math-pow x (math-add a -1))
547 (math-pow (math-sub 1 x) (math-add b -1)))
548 scale))
550 (put 'calcFunc-erf\' 'math-derivative-1
551 (function (lambda (x) (math-div 2
552 (math-mul (list 'calcFunc-exp
553 (math-sqr x))
554 (if calc-symbolic-mode
555 '(calcFunc-sqrt
556 (var pi var-pi))
557 (math-sqrt-pi)))))))
559 (put 'calcFunc-erfc\' 'math-derivative-1
560 (function (lambda (x) (math-div -2
561 (math-mul (list 'calcFunc-exp
562 (math-sqr x))
563 (if calc-symbolic-mode
564 '(calcFunc-sqrt
565 (var pi var-pi))
566 (math-sqrt-pi)))))))
568 (put 'calcFunc-besJ\'2 'math-derivative-2
569 (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besJ
570 (math-add v -1)
572 (list 'calcFunc-besJ
573 (math-add v 1)
575 2))))
577 (put 'calcFunc-besY\'2 'math-derivative-2
578 (function (lambda (v z) (math-div (math-sub (list 'calcFunc-besY
579 (math-add v -1)
581 (list 'calcFunc-besY
582 (math-add v 1)
584 2))))
586 (put 'calcFunc-sum 'math-derivative-n
587 (function
588 (lambda (expr)
589 (if (math-expr-contains (cons 'vec (cdr (cdr expr))) math-deriv-var)
590 (throw 'math-deriv nil)
591 (cons 'calcFunc-sum
592 (cons (math-derivative (nth 1 expr))
593 (cdr (cdr expr))))))))
595 (put 'calcFunc-prod 'math-derivative-n
596 (function
597 (lambda (expr)
598 (if (math-expr-contains (cons 'vec (cdr (cdr expr))) math-deriv-var)
599 (throw 'math-deriv nil)
600 (math-mul expr
601 (cons 'calcFunc-sum
602 (cons (math-div (math-derivative (nth 1 expr))
603 (nth 1 expr))
604 (cdr (cdr expr)))))))))
606 (put 'calcFunc-integ 'math-derivative-n
607 (function
608 (lambda (expr)
609 (if (= (length expr) 3)
610 (if (equal (nth 2 expr) math-deriv-var)
611 (nth 1 expr)
612 (math-normalize
613 (list 'calcFunc-integ
614 (math-derivative (nth 1 expr))
615 (nth 2 expr))))
616 (if (= (length expr) 5)
617 (let ((lower (math-expr-subst (nth 1 expr) (nth 2 expr)
618 (nth 3 expr)))
619 (upper (math-expr-subst (nth 1 expr) (nth 2 expr)
620 (nth 4 expr))))
621 (math-add (math-sub (math-mul upper
622 (math-derivative (nth 4 expr)))
623 (math-mul lower
624 (math-derivative (nth 3 expr))))
625 (if (equal (nth 2 expr) math-deriv-var)
627 (math-normalize
628 (list 'calcFunc-integ
629 (math-derivative (nth 1 expr)) (nth 2 expr)
630 (nth 3 expr) (nth 4 expr)))))))))))
632 (put 'calcFunc-if 'math-derivative-n
633 (function
634 (lambda (expr)
635 (and (= (length expr) 4)
636 (list 'calcFunc-if (nth 1 expr)
637 (math-derivative (nth 2 expr))
638 (math-derivative (nth 3 expr)))))))
640 (put 'calcFunc-subscr 'math-derivative-n
641 (function
642 (lambda (expr)
643 (and (= (length expr) 3)
644 (list 'calcFunc-subscr (nth 1 expr)
645 (math-derivative (nth 2 expr)))))))
648 (defvar math-integ-var '(var X ---))
649 (defvar math-integ-var-2 '(var Y ---))
650 (defvar math-integ-vars (list 'f math-integ-var math-integ-var-2))
651 (defvar math-integ-var-list (list math-integ-var))
652 (defvar math-integ-var-list-list (list math-integ-var-list))
654 ;; math-integ-depth is a local variable for math-try-integral, but is used
655 ;; by math-integral and math-tracing-integral
656 ;; which are called (directly or indirectly) by math-try-integral.
657 (defvar math-integ-depth)
658 ;; math-integ-level is a local variable for math-try-integral, but is used
659 ;; by math-integral, math-do-integral, math-tracing-integral,
660 ;; math-sub-integration, math-integrate-by-parts and
661 ;; math-integrate-by-substitution, which are called (directly or
662 ;; indirectly) by math-try-integral.
663 (defvar math-integ-level)
664 ;; math-integral-limit is a local variable for calcFunc-integ, but is
665 ;; used by math-tracing-integral, math-sub-integration and
666 ;; math-try-integration.
667 (defvar math-integral-limit)
669 (defmacro math-tracing-integral (&rest parts)
670 `(and trace-buffer
671 (with-current-buffer trace-buffer
672 (goto-char (point-max))
673 (and (bolp)
674 (insert (make-string (- math-integral-limit
675 math-integ-level) 32)
676 (format "%2d " math-integ-depth)
677 (make-string math-integ-level 32)))
678 ;;(condition-case err
679 (insert ,@parts)
680 ;; (error (insert (prin1-to-string err))))
681 (sit-for 0))))
683 ;;; The following wrapper caches results and avoids infinite recursion.
684 ;;; Each cache entry is: ( A B ) Integral of A is B;
685 ;;; ( A N ) Integral of A failed at level N;
686 ;;; ( A busy ) Currently working on integral of A;
687 ;;; ( A parts ) Currently working, integ-by-parts;
688 ;;; ( A parts2 ) Currently working, integ-by-parts;
689 ;;; ( A cancelled ) Ignore this cache entry;
690 ;;; ( A [B] ) Same result as for math-cur-record = B.
692 ;; math-cur-record is a local variable for math-try-integral, but is used
693 ;; by math-integral, math-replace-integral-parts and math-integrate-by-parts
694 ;; which are called (directly or indirectly) by math-try-integral, as well as
695 ;; by calc-dump-integral-cache
696 (defvar math-cur-record)
697 ;; math-enable-subst and math-any-substs are local variables for
698 ;; calcFunc-integ, but are used by math-integral and math-try-integral.
699 (defvar math-enable-subst)
700 (defvar math-any-substs)
702 ;; math-integ-msg is a local variable for math-try-integral, but is
703 ;; used (both locally and non-locally) by math-integral.
704 (defvar math-integ-msg)
706 (defvar math-integral-cache nil)
707 (defvar math-integral-cache-state nil)
709 (defun math-integral (expr &optional simplify same-as-above)
710 (let* ((simp math-cur-record)
711 (math-cur-record (assoc expr math-integral-cache))
712 (math-integ-depth (1+ math-integ-depth))
713 (val 'cancelled))
714 (math-tracing-integral "Integrating "
715 (math-format-value expr 1000)
716 "...\n")
717 (and math-cur-record
718 (progn
719 (math-tracing-integral "Found "
720 (math-format-value (nth 1 math-cur-record) 1000))
721 (and (consp (nth 1 math-cur-record))
722 (math-replace-integral-parts math-cur-record))
723 (math-tracing-integral " => "
724 (math-format-value (nth 1 math-cur-record) 1000)
725 "\n")))
726 (or (and math-cur-record
727 (not (eq (nth 1 math-cur-record) 'cancelled))
728 (or (not (integerp (nth 1 math-cur-record)))
729 (>= (nth 1 math-cur-record) math-integ-level)))
730 (and (math-integral-contains-parts expr)
731 (progn
732 (setq val nil)
734 (unwind-protect
735 (progn
736 (let (math-integ-msg)
737 (if (eq calc-display-working-message 'lots)
738 (progn
739 (calc-set-command-flag 'clear-message)
740 (setq math-integ-msg (format
741 "Working... Integrating %s"
742 (math-format-flat-expr expr 0)))
743 (message "%s" math-integ-msg)))
744 (if math-cur-record
745 (setcar (cdr math-cur-record)
746 (if same-as-above (vector simp) 'busy))
747 (setq math-cur-record
748 (list expr (if same-as-above (vector simp) 'busy))
749 math-integral-cache (cons math-cur-record
750 math-integral-cache)))
751 (if (eq simplify 'yes)
752 (progn
753 (math-tracing-integral "Simplifying...")
754 (setq simp (math-simplify expr))
755 (setq val (if (equal simp expr)
756 (progn
757 (math-tracing-integral " no change\n")
758 (math-do-integral expr))
759 (math-tracing-integral " simplified\n")
760 (math-integral simp 'no t))))
761 (or (setq val (math-do-integral expr))
762 (eq simplify 'no)
763 (let ((simp (math-simplify expr)))
764 (or (equal simp expr)
765 (progn
766 (math-tracing-integral "Trying again after "
767 "simplification...\n")
768 (setq val (math-integral simp 'no t))))))))
769 (if (eq calc-display-working-message 'lots)
770 (message "%s" math-integ-msg)))
771 (setcar (cdr math-cur-record) (or val
772 (if (or math-enable-subst
773 (not math-any-substs))
774 math-integ-level
775 'cancelled)))))
776 (setq val math-cur-record)
777 (while (vectorp (nth 1 val))
778 (setq val (aref (nth 1 val) 0)))
779 (setq val (if (memq (nth 1 val) '(parts parts2))
780 (progn
781 (setcar (cdr val) 'parts2)
782 (list 'var 'PARTS val))
783 (and (consp (nth 1 val))
784 (nth 1 val))))
785 (math-tracing-integral "Integral of "
786 (math-format-value expr 1000)
787 " is "
788 (math-format-value val 1000)
789 "\n")
790 val))
792 (defun math-integral-contains-parts (expr)
793 (if (Math-primp expr)
794 (and (eq (car-safe expr) 'var)
795 (eq (nth 1 expr) 'PARTS)
796 (listp (nth 2 expr)))
797 (while (and (setq expr (cdr expr))
798 (not (math-integral-contains-parts (car expr)))))
799 expr))
801 (defun math-replace-integral-parts (expr)
802 (or (Math-primp expr)
803 (while (setq expr (cdr expr))
804 (and (consp (car expr))
805 (if (eq (car (car expr)) 'var)
806 (and (eq (nth 1 (car expr)) 'PARTS)
807 (consp (nth 2 (car expr)))
808 (if (listp (nth 1 (nth 2 (car expr))))
809 (progn
810 (setcar expr (nth 1 (nth 2 (car expr))))
811 (math-replace-integral-parts (cons 'foo expr)))
812 (setcar (cdr math-cur-record) 'cancelled)))
813 (math-replace-integral-parts (car expr)))))))
815 (defvar math-linear-subst-tried t
816 "Non-nil means that a linear substitution has been tried.")
818 ;; The variable math-has-rules is a local variable for math-try-integral,
819 ;; but is used by math-do-integral, which is called (non-directly) by
820 ;; math-try-integral.
821 (defvar math-has-rules)
823 ;; math-old-integ is a local variable for math-do-integral, but is
824 ;; used by math-sub-integration.
825 (defvar math-old-integ)
827 ;; The variables math-t1, math-t2 and math-t3 are local to
828 ;; math-do-integral, math-try-solve-for and math-decompose-poly, but
829 ;; are used by functions they call (directly or indirectly);
830 ;; math-do-integral calls math-do-integral-methods;
831 ;; math-try-solve-for calls math-try-solve-prod,
832 ;; math-solve-find-root-term and math-solve-find-root-in-prod;
833 ;; math-decompose-poly calls math-solve-poly-funny-powers and
834 ;; math-solve-crunch-poly.
835 (defvar math-t1)
836 (defvar math-t2)
837 (defvar math-t3)
839 (defun math-do-integral (expr)
840 (let ((math-linear-subst-tried nil)
841 math-t1 math-t2)
842 (or (cond ((not (math-expr-contains expr math-integ-var))
843 (math-mul expr math-integ-var))
844 ((equal expr math-integ-var)
845 (math-div (math-sqr expr) 2))
846 ((eq (car expr) '+)
847 (and (setq math-t1 (math-integral (nth 1 expr)))
848 (setq math-t2 (math-integral (nth 2 expr)))
849 (math-add math-t1 math-t2)))
850 ((eq (car expr) '-)
851 (and (setq math-t1 (math-integral (nth 1 expr)))
852 (setq math-t2 (math-integral (nth 2 expr)))
853 (math-sub math-t1 math-t2)))
854 ((eq (car expr) 'neg)
855 (and (setq math-t1 (math-integral (nth 1 expr)))
856 (math-neg math-t1)))
857 ((eq (car expr) '*)
858 (cond ((not (math-expr-contains (nth 1 expr) math-integ-var))
859 (and (setq math-t1 (math-integral (nth 2 expr)))
860 (math-mul (nth 1 expr) math-t1)))
861 ((not (math-expr-contains (nth 2 expr) math-integ-var))
862 (and (setq math-t1 (math-integral (nth 1 expr)))
863 (math-mul math-t1 (nth 2 expr))))
864 ((memq (car-safe (nth 1 expr)) '(+ -))
865 (math-integral (list (car (nth 1 expr))
866 (math-mul (nth 1 (nth 1 expr))
867 (nth 2 expr))
868 (math-mul (nth 2 (nth 1 expr))
869 (nth 2 expr)))
870 'yes t))
871 ((memq (car-safe (nth 2 expr)) '(+ -))
872 (math-integral (list (car (nth 2 expr))
873 (math-mul (nth 1 (nth 2 expr))
874 (nth 1 expr))
875 (math-mul (nth 2 (nth 2 expr))
876 (nth 1 expr)))
877 'yes t))))
878 ((eq (car expr) '/)
879 (cond ((and (not (math-expr-contains (nth 1 expr)
880 math-integ-var))
881 (not (math-equal-int (nth 1 expr) 1)))
882 (and (setq math-t1 (math-integral (math-div 1 (nth 2 expr))))
883 (math-mul (nth 1 expr) math-t1)))
884 ((not (math-expr-contains (nth 2 expr) math-integ-var))
885 (and (setq math-t1 (math-integral (nth 1 expr)))
886 (math-div math-t1 (nth 2 expr))))
887 ((and (eq (car-safe (nth 1 expr)) '*)
888 (not (math-expr-contains (nth 1 (nth 1 expr))
889 math-integ-var)))
890 (and (setq math-t1 (math-integral
891 (math-div (nth 2 (nth 1 expr))
892 (nth 2 expr))))
893 (math-mul math-t1 (nth 1 (nth 1 expr)))))
894 ((and (eq (car-safe (nth 1 expr)) '*)
895 (not (math-expr-contains (nth 2 (nth 1 expr))
896 math-integ-var)))
897 (and (setq math-t1 (math-integral
898 (math-div (nth 1 (nth 1 expr))
899 (nth 2 expr))))
900 (math-mul math-t1 (nth 2 (nth 1 expr)))))
901 ((and (eq (car-safe (nth 2 expr)) '*)
902 (not (math-expr-contains (nth 1 (nth 2 expr))
903 math-integ-var)))
904 (and (setq math-t1 (math-integral
905 (math-div (nth 1 expr)
906 (nth 2 (nth 2 expr)))))
907 (math-div math-t1 (nth 1 (nth 2 expr)))))
908 ((and (eq (car-safe (nth 2 expr)) '*)
909 (not (math-expr-contains (nth 2 (nth 2 expr))
910 math-integ-var)))
911 (and (setq math-t1 (math-integral
912 (math-div (nth 1 expr)
913 (nth 1 (nth 2 expr)))))
914 (math-div math-t1 (nth 2 (nth 2 expr)))))
915 ((eq (car-safe (nth 2 expr)) 'calcFunc-exp)
916 (math-integral
917 (math-mul (nth 1 expr)
918 (list 'calcFunc-exp
919 (math-neg (nth 1 (nth 2 expr)))))))))
920 ((eq (car expr) '^)
921 (cond ((not (math-expr-contains (nth 1 expr) math-integ-var))
922 (or (and (setq math-t1 (math-is-polynomial (nth 2 expr)
923 math-integ-var 1))
924 (math-div expr
925 (math-mul (nth 1 math-t1)
926 (math-normalize
927 (list 'calcFunc-ln
928 (nth 1 expr))))))
929 (math-integral
930 (list 'calcFunc-exp
931 (math-mul (nth 2 expr)
932 (math-normalize
933 (list 'calcFunc-ln
934 (nth 1 expr)))))
935 'yes t)))
936 ((not (math-expr-contains (nth 2 expr) math-integ-var))
937 (if (and (integerp (nth 2 expr)) (< (nth 2 expr) 0))
938 (math-integral
939 (list '/ 1 (math-pow (nth 1 expr) (- (nth 2 expr))))
940 nil t)
941 (or (and (setq math-t1 (math-is-polynomial (nth 1 expr)
942 math-integ-var
944 (setq math-t2 (math-add (nth 2 expr) 1))
945 (math-div (math-pow (nth 1 expr) math-t2)
946 (math-mul math-t2 (nth 1 math-t1))))
947 (and (Math-negp (nth 2 expr))
948 (math-integral
949 (math-div 1
950 (math-pow (nth 1 expr)
951 (math-neg
952 (nth 2 expr))))
953 nil t))
954 nil))))))
956 ;; Integral of a polynomial.
957 (and (setq math-t1 (math-is-polynomial expr math-integ-var 20))
958 (let ((accum 0)
959 (n 1))
960 (while math-t1
961 (if (setq accum (math-add accum
962 (math-div (math-mul (car math-t1)
963 (math-pow
964 math-integ-var
967 math-t1 (cdr math-t1))
968 (setq n (1+ n))))
969 accum))
971 ;; Try looking it up!
972 (cond ((= (length expr) 2)
973 (and (symbolp (car expr))
974 (setq math-t1 (get (car expr) 'math-integral))
975 (progn
976 (while (and math-t1
977 (not (setq math-t2 (funcall (car math-t1)
978 (nth 1 expr)))))
979 (setq math-t1 (cdr math-t1)))
980 (and math-t2 (math-normalize math-t2)))))
981 ((= (length expr) 3)
982 (and (symbolp (car expr))
983 (setq math-t1 (get (car expr) 'math-integral-2))
984 (progn
985 (while (and math-t1
986 (not (setq math-t2 (funcall (car math-t1)
987 (nth 1 expr)
988 (nth 2 expr)))))
989 (setq math-t1 (cdr math-t1)))
990 (and math-t2 (math-normalize math-t2))))))
992 ;; Integral of a rational function.
993 (and (math-ratpoly-p expr math-integ-var)
994 (setq math-t1 (calcFunc-apart expr math-integ-var))
995 (not (equal math-t1 expr))
996 (math-integral math-t1))
998 ;; Try user-defined integration rules.
999 (and math-has-rules
1000 (let ((math-old-integ (symbol-function 'calcFunc-integ))
1001 (input (list 'calcFunc-integtry expr math-integ-var))
1002 res part)
1003 (unwind-protect
1004 (progn
1005 (fset 'calcFunc-integ 'math-sub-integration)
1006 (setq res (math-rewrite input
1007 '(var IntegRules var-IntegRules)
1009 (fset 'calcFunc-integ math-old-integ)
1010 (and (not (equal res input))
1011 (if (setq part (math-expr-calls
1012 res '(calcFunc-integsubst)))
1013 (and (memq (length part) '(3 4 5))
1014 (let ((parts (mapcar
1015 (function
1016 (lambda (x)
1017 (math-expr-subst
1018 x (nth 2 part)
1019 math-integ-var)))
1020 (cdr part))))
1021 (math-integrate-by-substitution
1022 expr (car parts) t
1023 (or (nth 2 parts)
1024 (list 'calcFunc-integfailed
1025 math-integ-var))
1026 (nth 3 parts))))
1027 (if (not (math-expr-calls res
1028 '(calcFunc-integtry
1029 calcFunc-integfailed)))
1030 res))))
1031 (fset 'calcFunc-integ math-old-integ))))
1033 ;; See if the function is a symbolic derivative.
1034 (and (string-match "'" (symbol-name (car expr)))
1035 (let ((name (symbol-name (car expr)))
1036 (p expr) (n 0) (which nil) (bad nil))
1037 (while (setq n (1+ n) p (cdr p))
1038 (if (equal (car p) math-integ-var)
1039 (if which (setq bad t) (setq which n))
1040 (if (math-expr-contains (car p) math-integ-var)
1041 (setq bad t))))
1042 (and which (not bad)
1043 (let ((prime (if (= which 1) "'" (format "'%d" which))))
1044 (and (string-match (concat prime "\\('['0-9]*\\|$\\)")
1045 name)
1046 (cons (intern
1047 (concat
1048 (substring name 0 (match-beginning 0))
1049 (substring name (+ (match-beginning 0)
1050 (length prime)))))
1051 (cdr expr)))))))
1053 ;; Try transformation methods (parts, substitutions).
1054 (and (> math-integ-level 0)
1055 (math-do-integral-methods expr))
1057 ;; Try expanding the function's definition.
1058 (let ((res (math-expand-formula expr)))
1059 (and res
1060 (math-integral res))))))
1062 (defun math-sub-integration (expr &rest rest)
1063 (or (if (or (not rest)
1064 (and (< math-integ-level math-integral-limit)
1065 (eq (car rest) math-integ-var)))
1066 (math-integral expr)
1067 (let ((res (apply math-old-integ expr rest)))
1068 (and (or (= math-integ-level math-integral-limit)
1069 (not (math-expr-calls res 'calcFunc-integ)))
1070 res)))
1071 (list 'calcFunc-integfailed expr)))
1073 ;; math-so-far is a local variable for math-do-integral-methods, but
1074 ;; is used by math-integ-try-linear-substitutions and
1075 ;; math-integ-try-substitutions.
1076 (defvar math-so-far)
1078 ;; math-integ-expr is a local variable for math-do-integral-methods,
1079 ;; but is used by math-integ-try-linear-substitutions and
1080 ;; math-integ-try-substitutions.
1081 (defvar math-integ-expr)
1083 (defun math-do-integral-methods (math-integ-expr)
1084 (let ((math-so-far math-integ-var-list-list)
1085 rat-in)
1087 ;; Integration by substitution, for various likely sub-expressions.
1088 ;; (In first pass, we look only for sub-exprs that are linear in X.)
1089 (or (math-integ-try-linear-substitutions math-integ-expr)
1090 (math-integ-try-substitutions math-integ-expr)
1092 ;; If function has sines and cosines, try tan(x/2) substitution.
1093 (and (let ((p (setq rat-in (math-expr-rational-in math-integ-expr))))
1094 (while (and p
1095 (memq (car (car p)) '(calcFunc-sin
1096 calcFunc-cos
1097 calcFunc-tan
1098 calcFunc-sec
1099 calcFunc-csc
1100 calcFunc-cot))
1101 (equal (nth 1 (car p)) math-integ-var))
1102 (setq p (cdr p)))
1103 (null p))
1104 (or (and (math-integ-parts-easy math-integ-expr)
1105 (math-integ-try-parts math-integ-expr t))
1106 (math-integrate-by-good-substitution
1107 math-integ-expr (list 'calcFunc-tan (math-div math-integ-var 2)))))
1109 ;; If function has sinh and cosh, try tanh(x/2) substitution.
1110 (and (let ((p rat-in))
1111 (while (and p
1112 (memq (car (car p)) '(calcFunc-sinh
1113 calcFunc-cosh
1114 calcFunc-tanh
1115 calcFunc-sech
1116 calcFunc-csch
1117 calcFunc-coth
1118 calcFunc-exp))
1119 (equal (nth 1 (car p)) math-integ-var))
1120 (setq p (cdr p)))
1121 (null p))
1122 (or (and (math-integ-parts-easy math-integ-expr)
1123 (math-integ-try-parts math-integ-expr t))
1124 (math-integrate-by-good-substitution
1125 math-integ-expr (list 'calcFunc-tanh (math-div math-integ-var 2)))))
1127 ;; If function has square roots, try sin, tan, or sec substitution.
1128 (and (let ((p rat-in))
1129 (setq math-t1 nil)
1130 (while (and p
1131 (or (equal (car p) math-integ-var)
1132 (and (eq (car (car p)) 'calcFunc-sqrt)
1133 (setq math-t1 (math-is-polynomial
1134 (nth 1 (setq math-t2 (car p)))
1135 math-integ-var 2)))))
1136 (setq p (cdr p)))
1137 (and (null p) math-t1))
1138 (if (cdr (cdr math-t1))
1139 (if (math-guess-if-neg (nth 2 math-t1))
1140 (let* ((c (math-sqrt (math-neg (nth 2 math-t1))))
1141 (d (math-div (nth 1 math-t1) (math-mul -2 c)))
1142 (a (math-sqrt (math-add (car math-t1) (math-sqr d)))))
1143 (math-integrate-by-good-substitution
1144 math-integ-expr (list 'calcFunc-arcsin
1145 (math-div-thru
1146 (math-add (math-mul c math-integ-var) d)
1147 a))))
1148 (let* ((c (math-sqrt (nth 2 math-t1)))
1149 (d (math-div (nth 1 math-t1) (math-mul 2 c)))
1150 (aa (math-sub (car math-t1) (math-sqr d))))
1151 (if (and nil (not (and (eq d 0) (eq c 1))))
1152 (math-integrate-by-good-substitution
1153 math-integ-expr (math-add (math-mul c math-integ-var) d))
1154 (if (math-guess-if-neg aa)
1155 (math-integrate-by-good-substitution
1156 math-integ-expr (list 'calcFunc-arccosh
1157 (math-div-thru
1158 (math-add (math-mul c math-integ-var)
1160 (math-sqrt (math-neg aa)))))
1161 (math-integrate-by-good-substitution
1162 math-integ-expr (list 'calcFunc-arcsinh
1163 (math-div-thru
1164 (math-add (math-mul c math-integ-var)
1166 (math-sqrt aa))))))))
1167 (math-integrate-by-good-substitution math-integ-expr math-t2)) )
1169 ;; Try integration by parts.
1170 (math-integ-try-parts math-integ-expr)
1172 ;; Give up.
1173 nil)))
1175 (defun math-integ-parts-easy (expr)
1176 (cond ((Math-primp expr) t)
1177 ((memq (car expr) '(+ - *))
1178 (and (math-integ-parts-easy (nth 1 expr))
1179 (math-integ-parts-easy (nth 2 expr))))
1180 ((eq (car expr) '/)
1181 (and (math-integ-parts-easy (nth 1 expr))
1182 (math-atomic-factorp (nth 2 expr))))
1183 ((eq (car expr) '^)
1184 (and (natnump (nth 2 expr))
1185 (math-integ-parts-easy (nth 1 expr))))
1186 ((eq (car expr) 'neg)
1187 (math-integ-parts-easy (nth 1 expr)))
1188 (t t)))
1190 ;; math-prev-parts-v is local to calcFunc-integ (as well as
1191 ;; math-integrate-by-parts), but is used by math-integ-try-parts.
1192 (defvar math-prev-parts-v)
1194 ;; math-good-parts is local to calcFunc-integ (as well as
1195 ;; math-integ-try-parts), but is used by math-integrate-by-parts.
1196 (defvar math-good-parts)
1199 (defun math-integ-try-parts (expr &optional math-good-parts)
1200 ;; Integration by parts:
1201 ;; integ(f(x) g(x),x) = f(x) h(x) - integ(h(x) f'(x),x)
1202 ;; where h(x) = integ(g(x),x).
1203 (or (let ((exp (calcFunc-expand expr)))
1204 (and (not (equal exp expr))
1205 (math-integral exp)))
1206 (and (eq (car expr) '*)
1207 (let ((first-bad (or (math-polynomial-p (nth 1 expr)
1208 math-integ-var)
1209 (equal (nth 2 expr) math-prev-parts-v))))
1210 (or (and first-bad ; so try this one first
1211 (math-integrate-by-parts (nth 1 expr) (nth 2 expr)))
1212 (math-integrate-by-parts (nth 2 expr) (nth 1 expr))
1213 (and (not first-bad)
1214 (math-integrate-by-parts (nth 1 expr) (nth 2 expr))))))
1215 (and (eq (car expr) '/)
1216 (math-expr-contains (nth 1 expr) math-integ-var)
1217 (let ((recip (math-div 1 (nth 2 expr))))
1218 (or (math-integrate-by-parts (nth 1 expr) recip)
1219 (math-integrate-by-parts recip (nth 1 expr)))))
1220 (and (eq (car expr) '^)
1221 (math-integrate-by-parts (math-pow (nth 1 expr)
1222 (math-sub (nth 2 expr) 1))
1223 (nth 1 expr)))))
1225 (defun math-integrate-by-parts (u vprime)
1226 (let ((math-integ-level (if (or math-good-parts
1227 (math-polynomial-p u math-integ-var))
1228 math-integ-level
1229 (1- math-integ-level)))
1230 (math-doing-parts t)
1231 v temp)
1232 (and (>= math-integ-level 0)
1233 (unwind-protect
1234 (progn
1235 (setcar (cdr math-cur-record) 'parts)
1236 (math-tracing-integral "Integrating by parts, u = "
1237 (math-format-value u 1000)
1238 ", v' = "
1239 (math-format-value vprime 1000)
1240 "\n")
1241 (and (setq v (math-integral vprime))
1242 (setq temp (calcFunc-deriv u math-integ-var nil t))
1243 (setq temp (let ((math-prev-parts-v v))
1244 (math-integral (math-mul v temp) 'yes)))
1245 (setq temp (math-sub (math-mul u v) temp))
1246 (if (eq (nth 1 math-cur-record) 'parts)
1247 (calcFunc-expand temp)
1248 (setq v (list 'var 'PARTS math-cur-record)
1249 temp (let (calc-next-why)
1250 (math-simplify-extended
1251 (math-solve-for (math-sub v temp) 0 v nil)))
1252 temp (if (and (eq (car-safe temp) '/)
1253 (math-zerop (nth 2 temp)))
1254 nil temp)))))
1255 (setcar (cdr math-cur-record) 'busy)))))
1257 ;;; This tries two different formulations, hoping the algebraic simplifier
1258 ;;; will be strong enough to handle at least one.
1259 (defun math-integrate-by-substitution (expr u &optional user uinv uinvprime)
1260 (and (> math-integ-level 0)
1261 (let ((math-integ-level (max (- math-integ-level 2) 0)))
1262 (math-integrate-by-good-substitution expr u user uinv uinvprime))))
1264 (defun math-integrate-by-good-substitution (expr u &optional user
1265 uinv uinvprime)
1266 (let ((math-living-dangerously t)
1267 deriv temp)
1268 (and (setq uinv (if uinv
1269 (math-expr-subst uinv math-integ-var
1270 math-integ-var-2)
1271 (let (calc-next-why)
1272 (math-solve-for u
1273 math-integ-var-2
1274 math-integ-var nil))))
1275 (progn
1276 (math-tracing-integral "Integrating by substitution, u = "
1277 (math-format-value u 1000)
1278 "\n")
1279 (or (and (setq deriv (calcFunc-deriv u
1280 math-integ-var nil
1281 (not user)))
1282 (setq temp (math-integral (math-expr-subst
1283 (math-expr-subst
1284 (math-expr-subst
1285 (math-div expr deriv)
1287 math-integ-var-2)
1288 math-integ-var
1289 uinv)
1290 math-integ-var-2
1291 math-integ-var)
1292 'yes)))
1293 (and (setq deriv (or uinvprime
1294 (calcFunc-deriv uinv
1295 math-integ-var-2
1296 math-integ-var
1297 (not user))))
1298 (setq temp (math-integral (math-mul
1299 (math-expr-subst
1300 (math-expr-subst
1301 (math-expr-subst
1302 expr
1304 math-integ-var-2)
1305 math-integ-var
1306 uinv)
1307 math-integ-var-2
1308 math-integ-var)
1309 deriv)
1310 'yes)))))
1311 (math-simplify-extended
1312 (math-expr-subst temp math-integ-var u)))))
1314 ;;; Look for substitutions of the form u = a x + b.
1315 (defun math-integ-try-linear-substitutions (sub-expr)
1316 (setq math-linear-subst-tried t)
1317 (and (not (Math-primp sub-expr))
1318 (or (and (not (memq (car sub-expr) '(+ - * / neg)))
1319 (not (and (eq (car sub-expr) '^)
1320 (integerp (nth 2 sub-expr))))
1321 (math-expr-contains sub-expr math-integ-var)
1322 (let ((res nil))
1323 (while (and (setq sub-expr (cdr sub-expr))
1324 (or (not (math-linear-in (car sub-expr)
1325 math-integ-var))
1326 (assoc (car sub-expr) math-so-far)
1327 (progn
1328 (setq math-so-far (cons (list (car sub-expr))
1329 math-so-far))
1330 (not (setq res
1331 (math-integrate-by-substitution
1332 math-integ-expr (car sub-expr))))))))
1333 res))
1334 (let ((res nil))
1335 (while (and (setq sub-expr (cdr sub-expr))
1336 (not (setq res (math-integ-try-linear-substitutions
1337 (car sub-expr))))))
1338 res))))
1340 ;;; Recursively try different substitutions based on various sub-expressions.
1341 (defun math-integ-try-substitutions (sub-expr &optional allow-rat)
1342 (and (not (Math-primp sub-expr))
1343 (not (assoc sub-expr math-so-far))
1344 (math-expr-contains sub-expr math-integ-var)
1345 (or (and (if (and (not (memq (car sub-expr) '(+ - * / neg)))
1346 (not (and (eq (car sub-expr) '^)
1347 (integerp (nth 2 sub-expr)))))
1348 (setq allow-rat t)
1349 (prog1 allow-rat (setq allow-rat nil)))
1350 (not (eq sub-expr math-integ-expr))
1351 (or (math-integrate-by-substitution math-integ-expr sub-expr)
1352 (and (eq (car sub-expr) '^)
1353 (integerp (nth 2 sub-expr))
1354 (< (nth 2 sub-expr) 0)
1355 (math-integ-try-substitutions
1356 (math-pow (nth 1 sub-expr) (- (nth 2 sub-expr)))
1357 t))))
1358 (let ((res nil))
1359 (setq math-so-far (cons (list sub-expr) math-so-far))
1360 (while (and (setq sub-expr (cdr sub-expr))
1361 (not (setq res (math-integ-try-substitutions
1362 (car sub-expr) allow-rat)))))
1363 res))))
1365 ;; The variable math-expr-parts is local to math-expr-rational-in,
1366 ;; but is used by math-expr-rational-in-rec
1367 (defvar math-expr-parts)
1369 (defun math-expr-rational-in (expr)
1370 (let ((math-expr-parts nil))
1371 (math-expr-rational-in-rec expr)
1372 (mapcar 'car math-expr-parts)))
1374 (defun math-expr-rational-in-rec (expr)
1375 (cond ((Math-primp expr)
1376 (and (equal expr math-integ-var)
1377 (not (assoc expr math-expr-parts))
1378 (setq math-expr-parts (cons (list expr) math-expr-parts))))
1379 ((or (memq (car expr) '(+ - * / neg))
1380 (and (eq (car expr) '^) (integerp (nth 2 expr))))
1381 (math-expr-rational-in-rec (nth 1 expr))
1382 (and (nth 2 expr) (math-expr-rational-in-rec (nth 2 expr))))
1383 ((and (eq (car expr) '^)
1384 (eq (math-quarter-integer (nth 2 expr)) 2))
1385 (math-expr-rational-in-rec (list 'calcFunc-sqrt (nth 1 expr))))
1387 (and (not (assoc expr math-expr-parts))
1388 (math-expr-contains expr math-integ-var)
1389 (setq math-expr-parts (cons (list expr) math-expr-parts))))))
1391 (defun math-expr-calls (expr funcs &optional arg-contains)
1392 (if (consp expr)
1393 (if (or (memq (car expr) funcs)
1394 (and (eq (car expr) '^) (eq (car funcs) 'calcFunc-sqrt)
1395 (eq (math-quarter-integer (nth 2 expr)) 2)))
1396 (and (or (not arg-contains)
1397 (math-expr-contains expr arg-contains))
1398 expr)
1399 (and (not (Math-primp expr))
1400 (let ((res nil))
1401 (while (and (setq expr (cdr expr))
1402 (not (setq res (math-expr-calls
1403 (car expr) funcs arg-contains)))))
1404 res)))))
1406 (defun math-fix-const-terms (expr except-vars)
1407 (cond ((not (math-expr-depends expr except-vars)) 0)
1408 ((Math-primp expr) expr)
1409 ((eq (car expr) '+)
1410 (math-add (math-fix-const-terms (nth 1 expr) except-vars)
1411 (math-fix-const-terms (nth 2 expr) except-vars)))
1412 ((eq (car expr) '-)
1413 (math-sub (math-fix-const-terms (nth 1 expr) except-vars)
1414 (math-fix-const-terms (nth 2 expr) except-vars)))
1415 (t expr)))
1417 ;; Command for debugging the Calculator's symbolic integrator.
1418 (defun calc-dump-integral-cache (&optional arg)
1419 (interactive "P")
1420 (let ((buf (current-buffer)))
1421 (unwind-protect
1422 (let ((p math-integral-cache)
1423 math-cur-record)
1424 (display-buffer (get-buffer-create "*Integral Cache*"))
1425 (set-buffer (get-buffer "*Integral Cache*"))
1426 (erase-buffer)
1427 (while p
1428 (setq math-cur-record (car p))
1429 (or arg (math-replace-integral-parts math-cur-record))
1430 (insert (math-format-flat-expr (car math-cur-record) 0)
1431 " --> "
1432 (if (symbolp (nth 1 math-cur-record))
1433 (concat "(" (symbol-name (nth 1 math-cur-record)) ")")
1434 (math-format-flat-expr (nth 1 math-cur-record) 0))
1435 "\n")
1436 (setq p (cdr p)))
1437 (goto-char (point-min)))
1438 (set-buffer buf))))
1440 ;; The variable math-max-integral-limit is local to calcFunc-integ,
1441 ;; but is used by math-try-integral.
1442 (defvar math-max-integral-limit)
1444 (defun math-try-integral (expr)
1445 (let ((math-integ-level math-integral-limit)
1446 (math-integ-depth 0)
1447 (math-integ-msg "Working...done")
1448 (math-cur-record nil) ; a technicality
1449 (math-integrating t)
1450 (calc-prefer-frac t)
1451 (calc-symbolic-mode t)
1452 (math-has-rules (calc-has-rules 'var-IntegRules)))
1453 (or (math-integral expr 'yes)
1454 (and math-any-substs
1455 (setq math-enable-subst t)
1456 (math-integral expr 'yes))
1457 (and (> math-max-integral-limit math-integral-limit)
1458 (setq math-integral-limit math-max-integral-limit
1459 math-integ-level math-integral-limit)
1460 (math-integral expr 'yes)))))
1462 (defvar var-IntegLimit nil)
1464 (defun calcFunc-integ (expr var &optional low high)
1465 (cond
1466 ;; Do these even if the parts turn out not to be integrable.
1467 ((eq (car-safe expr) '+)
1468 (math-add (calcFunc-integ (nth 1 expr) var low high)
1469 (calcFunc-integ (nth 2 expr) var low high)))
1470 ((eq (car-safe expr) '-)
1471 (math-sub (calcFunc-integ (nth 1 expr) var low high)
1472 (calcFunc-integ (nth 2 expr) var low high)))
1473 ((eq (car-safe expr) 'neg)
1474 (math-neg (calcFunc-integ (nth 1 expr) var low high)))
1475 ((and (eq (car-safe expr) '*)
1476 (not (math-expr-contains (nth 1 expr) var)))
1477 (math-mul (nth 1 expr) (calcFunc-integ (nth 2 expr) var low high)))
1478 ((and (eq (car-safe expr) '*)
1479 (not (math-expr-contains (nth 2 expr) var)))
1480 (math-mul (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr)))
1481 ((and (eq (car-safe expr) '/)
1482 (not (math-expr-contains (nth 1 expr) var))
1483 (not (math-equal-int (nth 1 expr) 1)))
1484 (math-mul (nth 1 expr)
1485 (calcFunc-integ (math-div 1 (nth 2 expr)) var low high)))
1486 ((and (eq (car-safe expr) '/)
1487 (not (math-expr-contains (nth 2 expr) var)))
1488 (math-div (calcFunc-integ (nth 1 expr) var low high) (nth 2 expr)))
1489 ((and (eq (car-safe expr) '/)
1490 (eq (car-safe (nth 1 expr)) '*)
1491 (not (math-expr-contains (nth 1 (nth 1 expr)) var)))
1492 (math-mul (nth 1 (nth 1 expr))
1493 (calcFunc-integ (math-div (nth 2 (nth 1 expr)) (nth 2 expr))
1494 var low high)))
1495 ((and (eq (car-safe expr) '/)
1496 (eq (car-safe (nth 1 expr)) '*)
1497 (not (math-expr-contains (nth 2 (nth 1 expr)) var)))
1498 (math-mul (nth 2 (nth 1 expr))
1499 (calcFunc-integ (math-div (nth 1 (nth 1 expr)) (nth 2 expr))
1500 var low high)))
1501 ((and (eq (car-safe expr) '/)
1502 (eq (car-safe (nth 2 expr)) '*)
1503 (not (math-expr-contains (nth 1 (nth 2 expr)) var)))
1504 (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 2 (nth 2 expr)))
1505 var low high)
1506 (nth 1 (nth 2 expr))))
1507 ((and (eq (car-safe expr) '/)
1508 (eq (car-safe (nth 2 expr)) '*)
1509 (not (math-expr-contains (nth 2 (nth 2 expr)) var)))
1510 (math-div (calcFunc-integ (math-div (nth 1 expr) (nth 1 (nth 2 expr)))
1511 var low high)
1512 (nth 2 (nth 2 expr))))
1513 ((eq (car-safe expr) 'vec)
1514 (cons 'vec (mapcar (function (lambda (x) (calcFunc-integ x var low high)))
1515 (cdr expr))))
1517 (let ((state (list calc-angle-mode
1518 ;;calc-symbolic-mode
1519 ;;calc-prefer-frac
1520 calc-internal-prec
1521 (calc-var-value 'var-IntegRules)
1522 (calc-var-value 'var-IntegSimpRules))))
1523 (or (equal state math-integral-cache-state)
1524 (setq math-integral-cache-state state
1525 math-integral-cache nil)))
1526 (let* ((math-max-integral-limit (or (and (natnump var-IntegLimit)
1527 var-IntegLimit)
1529 (math-integral-limit 1)
1530 (sexpr (math-expr-subst expr var math-integ-var))
1531 (trace-buffer (get-buffer "*Trace*"))
1532 (calc-language (if (eq calc-language 'big) nil calc-language))
1533 (math-any-substs t)
1534 (math-enable-subst nil)
1535 (math-prev-parts-v nil)
1536 (math-doing-parts nil)
1537 (math-good-parts nil)
1538 (res
1539 (if trace-buffer
1540 (let ((calcbuf (current-buffer))
1541 (calcwin (selected-window)))
1542 (unwind-protect
1543 (progn
1544 (if (get-buffer-window trace-buffer)
1545 (select-window (get-buffer-window trace-buffer)))
1546 (set-buffer trace-buffer)
1547 (goto-char (point-max))
1548 (or (assq 'scroll-stop (buffer-local-variables))
1549 (progn
1550 (make-local-variable 'scroll-step)
1551 (setq scroll-step 3)))
1552 (insert "\n\n\n")
1553 (set-buffer calcbuf)
1554 (math-try-integral sexpr))
1555 (select-window calcwin)
1556 (set-buffer calcbuf)))
1557 (math-try-integral sexpr))))
1558 (if res
1559 (progn
1560 (if (calc-has-rules 'var-IntegAfterRules)
1561 (setq res (math-rewrite res '(var IntegAfterRules
1562 var-IntegAfterRules))))
1563 (math-simplify
1564 (if (and low high)
1565 (math-sub (math-expr-subst res math-integ-var high)
1566 (math-expr-subst res math-integ-var low))
1567 (setq res (math-fix-const-terms res math-integ-vars))
1568 (if low
1569 (math-expr-subst res math-integ-var low)
1570 (math-expr-subst res math-integ-var var)))))
1571 (append (list 'calcFunc-integ expr var)
1572 (and low (list low))
1573 (and high (list high))))))))
1576 (math-defintegral calcFunc-inv
1577 (math-integral (math-div 1 u)))
1579 (math-defintegral calcFunc-conj
1580 (let ((int (math-integral u)))
1581 (and int
1582 (list 'calcFunc-conj int))))
1584 (math-defintegral calcFunc-deg
1585 (let ((int (math-integral u)))
1586 (and int
1587 (list 'calcFunc-deg int))))
1589 (math-defintegral calcFunc-rad
1590 (let ((int (math-integral u)))
1591 (and int
1592 (list 'calcFunc-rad int))))
1594 (math-defintegral calcFunc-re
1595 (let ((int (math-integral u)))
1596 (and int
1597 (list 'calcFunc-re int))))
1599 (math-defintegral calcFunc-im
1600 (let ((int (math-integral u)))
1601 (and int
1602 (list 'calcFunc-im int))))
1604 (math-defintegral calcFunc-sqrt
1605 (and (equal u math-integ-var)
1606 (math-mul '(frac 2 3)
1607 (list 'calcFunc-sqrt (math-pow u 3)))))
1609 (math-defintegral calcFunc-exp
1610 (or (and (equal u math-integ-var)
1611 (list 'calcFunc-exp u))
1612 (let ((p (math-is-polynomial u math-integ-var 2)))
1613 (and (nth 2 p)
1614 (let ((sqa (math-sqrt (math-neg (nth 2 p)))))
1615 (math-div
1616 (math-mul
1617 (math-mul (math-div (list 'calcFunc-sqrt '(var pi var-pi))
1618 sqa)
1619 (math-normalize
1620 (list 'calcFunc-exp
1621 (math-div (math-sub (math-mul (car p)
1622 (nth 2 p))
1623 (math-div
1624 (math-sqr (nth 1 p))
1626 (nth 2 p)))))
1627 (list 'calcFunc-erf
1628 (math-sub (math-mul sqa math-integ-var)
1629 (math-div (nth 1 p) (math-mul 2 sqa)))))
1630 2))))))
1632 (math-defintegral calcFunc-ln
1633 (or (and (equal u math-integ-var)
1634 (math-sub (math-mul u (list 'calcFunc-ln u)) u))
1635 (and (eq (car u) '*)
1636 (math-integral (math-add (list 'calcFunc-ln (nth 1 u))
1637 (list 'calcFunc-ln (nth 2 u)))))
1638 (and (eq (car u) '/)
1639 (math-integral (math-sub (list 'calcFunc-ln (nth 1 u))
1640 (list 'calcFunc-ln (nth 2 u)))))
1641 (and (eq (car u) '^)
1642 (math-integral (math-mul (nth 2 u)
1643 (list 'calcFunc-ln (nth 1 u)))))))
1645 (math-defintegral calcFunc-log10
1646 (and (equal u math-integ-var)
1647 (math-sub (math-mul u (list 'calcFunc-ln u))
1648 (math-div u (list 'calcFunc-ln 10)))))
1650 (math-defintegral-2 calcFunc-log
1651 (math-integral (math-div (list 'calcFunc-ln u)
1652 (list 'calcFunc-ln v))))
1654 (math-defintegral calcFunc-sin
1655 (or (and (equal u math-integ-var)
1656 (math-neg (math-from-radians-2 (list 'calcFunc-cos u))))
1657 (and (nth 2 (math-is-polynomial u math-integ-var 2))
1658 (math-integral (math-to-exponentials (list 'calcFunc-sin u))))))
1660 (math-defintegral calcFunc-cos
1661 (or (and (equal u math-integ-var)
1662 (math-from-radians-2 (list 'calcFunc-sin u)))
1663 (and (nth 2 (math-is-polynomial u math-integ-var 2))
1664 (math-integral (math-to-exponentials (list 'calcFunc-cos u))))))
1666 (math-defintegral calcFunc-tan
1667 (and (equal u math-integ-var)
1668 (math-from-radians-2
1669 (list 'calcFunc-ln (list 'calcFunc-sec u)))))
1671 (math-defintegral calcFunc-sec
1672 (and (equal u math-integ-var)
1673 (math-from-radians-2
1674 (list 'calcFunc-ln
1675 (math-add
1676 (list 'calcFunc-sec u)
1677 (list 'calcFunc-tan u))))))
1679 (math-defintegral calcFunc-csc
1680 (and (equal u math-integ-var)
1681 (math-from-radians-2
1682 (list 'calcFunc-ln
1683 (math-sub
1684 (list 'calcFunc-csc u)
1685 (list 'calcFunc-cot u))))))
1687 (math-defintegral calcFunc-cot
1688 (and (equal u math-integ-var)
1689 (math-from-radians-2
1690 (list 'calcFunc-ln (list 'calcFunc-sin u)))))
1692 (math-defintegral calcFunc-arcsin
1693 (and (equal u math-integ-var)
1694 (math-add (math-mul u (list 'calcFunc-arcsin u))
1695 (math-from-radians-2
1696 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))))
1698 (math-defintegral calcFunc-arccos
1699 (and (equal u math-integ-var)
1700 (math-sub (math-mul u (list 'calcFunc-arccos u))
1701 (math-from-radians-2
1702 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u)))))))
1704 (math-defintegral calcFunc-arctan
1705 (and (equal u math-integ-var)
1706 (math-sub (math-mul u (list 'calcFunc-arctan u))
1707 (math-from-radians-2
1708 (math-div (list 'calcFunc-ln (math-add 1 (math-sqr u)))
1709 2)))))
1711 (math-defintegral calcFunc-sinh
1712 (and (equal u math-integ-var)
1713 (list 'calcFunc-cosh u)))
1715 (math-defintegral calcFunc-cosh
1716 (and (equal u math-integ-var)
1717 (list 'calcFunc-sinh u)))
1719 (math-defintegral calcFunc-tanh
1720 (and (equal u math-integ-var)
1721 (list 'calcFunc-ln (list 'calcFunc-cosh u))))
1723 (math-defintegral calcFunc-sech
1724 (and (equal u math-integ-var)
1725 (list 'calcFunc-arctan (list 'calcFunc-sinh u))))
1727 (math-defintegral calcFunc-csch
1728 (and (equal u math-integ-var)
1729 (list 'calcFunc-ln (list 'calcFunc-tanh (math-div u 2)))))
1731 (math-defintegral calcFunc-coth
1732 (and (equal u math-integ-var)
1733 (list 'calcFunc-ln (list 'calcFunc-sinh u))))
1735 (math-defintegral calcFunc-arcsinh
1736 (and (equal u math-integ-var)
1737 (math-sub (math-mul u (list 'calcFunc-arcsinh u))
1738 (list 'calcFunc-sqrt (math-add (math-sqr u) 1)))))
1740 (math-defintegral calcFunc-arccosh
1741 (and (equal u math-integ-var)
1742 (math-sub (math-mul u (list 'calcFunc-arccosh u))
1743 (list 'calcFunc-sqrt (math-sub 1 (math-sqr u))))))
1745 (math-defintegral calcFunc-arctanh
1746 (and (equal u math-integ-var)
1747 (math-sub (math-mul u (list 'calcFunc-arctan u))
1748 (math-div (list 'calcFunc-ln
1749 (math-add 1 (math-sqr u)))
1750 2))))
1752 ;;; (Ax + B) / (ax^2 + bx + c)^n forms.
1753 (math-defintegral-2 /
1754 (math-integral-rational-funcs u v))
1756 (defun math-integral-rational-funcs (u v)
1757 (let ((pu (math-is-polynomial u math-integ-var 1))
1758 (vpow 1) pv)
1759 (and pu
1760 (catch 'int-rat
1761 (if (and (eq (car-safe v) '^) (natnump (nth 2 v)))
1762 (setq vpow (nth 2 v)
1763 v (nth 1 v)))
1764 (and (setq pv (math-is-polynomial v math-integ-var 2))
1765 (let ((int (math-mul-thru
1766 (car pu)
1767 (math-integral-q02 (car pv) (nth 1 pv)
1768 (nth 2 pv) v vpow))))
1769 (if (cdr pu)
1770 (setq int (math-add int
1771 (math-mul-thru
1772 (nth 1 pu)
1773 (math-integral-q12
1774 (car pv) (nth 1 pv)
1775 (nth 2 pv) v vpow)))))
1776 int))))))
1778 (defun math-integral-q12 (a b c v vpow)
1779 (let (q)
1780 (cond ((not c)
1781 (cond ((= vpow 1)
1782 (math-sub (math-div math-integ-var b)
1783 (math-mul (math-div a (math-sqr b))
1784 (list 'calcFunc-ln v))))
1785 ((= vpow 2)
1786 (math-div (math-add (list 'calcFunc-ln v)
1787 (math-div a v))
1788 (math-sqr b)))
1790 (let ((nm1 (math-sub vpow 1))
1791 (nm2 (math-sub vpow 2)))
1792 (math-div (math-sub
1793 (math-div a (math-mul nm1 (math-pow v nm1)))
1794 (math-div 1 (math-mul nm2 (math-pow v nm2))))
1795 (math-sqr b))))))
1796 ((math-zerop
1797 (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b))))
1798 (let ((part (math-div b (math-mul 2 c))))
1799 (math-mul-thru (math-pow c vpow)
1800 (math-integral-q12 part 1 nil
1801 (math-add math-integ-var part)
1802 (* vpow 2)))))
1803 ((= vpow 1)
1804 (and (math-ratp q) (math-negp q)
1805 (let ((calc-symbolic-mode t))
1806 (math-ratp (math-sqrt (math-neg q))))
1807 (throw 'int-rat nil)) ; should have used calcFunc-apart first
1808 (math-sub (math-div (list 'calcFunc-ln v) (math-mul 2 c))
1809 (math-mul-thru (math-div b (math-mul 2 c))
1810 (math-integral-q02 a b c v 1))))
1812 (let ((n (1- vpow)))
1813 (math-sub (math-neg (math-div
1814 (math-add (math-mul b math-integ-var)
1815 (math-mul 2 a))
1816 (math-mul n (math-mul q (math-pow v n)))))
1817 (math-mul-thru (math-div (math-mul b (1- (* 2 n)))
1818 (math-mul n q))
1819 (math-integral-q02 a b c v n))))))))
1821 (defun math-integral-q02 (a b c v vpow)
1822 (let (q rq part)
1823 (cond ((not c)
1824 (cond ((= vpow 1)
1825 (math-div (list 'calcFunc-ln v) b))
1827 (math-div (math-pow v (- 1 vpow))
1828 (math-mul (- 1 vpow) b)))))
1829 ((math-zerop
1830 (setq q (math-sub (math-mul 4 (math-mul a c)) (math-sqr b))))
1831 (let ((part (math-div b (math-mul 2 c))))
1832 (math-mul-thru (math-pow c vpow)
1833 (math-integral-q02 part 1 nil
1834 (math-add math-integ-var part)
1835 (* vpow 2)))))
1836 ((progn
1837 (setq part (math-add (math-mul 2 (math-mul c math-integ-var)) b))
1838 (> vpow 1))
1839 (let ((n (1- vpow)))
1840 (math-add (math-div part (math-mul n (math-mul q (math-pow v n))))
1841 (math-mul-thru (math-div (math-mul (- (* 4 n) 2) c)
1842 (math-mul n q))
1843 (math-integral-q02 a b c v n)))))
1844 ((math-guess-if-neg q)
1845 (setq rq (list 'calcFunc-sqrt (math-neg q)))
1846 ;;(math-div-thru (list 'calcFunc-ln
1847 ;; (math-div (math-sub part rq)
1848 ;; (math-add part rq)))
1849 ;; rq)
1850 (math-div (math-mul -2 (list 'calcFunc-arctanh
1851 (math-div part rq)))
1852 rq))
1854 (setq rq (list 'calcFunc-sqrt q))
1855 (math-div (math-mul 2 (math-to-radians-2
1856 (list 'calcFunc-arctan
1857 (math-div part rq))))
1858 rq)))))
1861 (math-defintegral calcFunc-erf
1862 (and (equal u math-integ-var)
1863 (math-add (math-mul u (list 'calcFunc-erf u))
1864 (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u))
1865 (list 'calcFunc-sqrt
1866 '(var pi var-pi)))))))
1868 (math-defintegral calcFunc-erfc
1869 (and (equal u math-integ-var)
1870 (math-sub (math-mul u (list 'calcFunc-erfc u))
1871 (math-div 1 (math-mul (list 'calcFunc-exp (math-sqr u))
1872 (list 'calcFunc-sqrt
1873 '(var pi var-pi)))))))
1878 (defvar math-tabulate-initial nil)
1879 (defvar math-tabulate-function nil)
1881 ;; These variables are local to calcFunc-table, but are used by
1882 ;; math-scan-for-limits.
1883 (defvar calc-low)
1884 (defvar calc-high)
1885 (defvar math-var)
1887 (defun calcFunc-table (expr math-var &optional calc-low calc-high step)
1888 (or calc-low
1889 (setq calc-low '(neg (var inf var-inf)) calc-high '(var inf var-inf)))
1890 (or calc-high (setq calc-high calc-low calc-low 1))
1891 (and (or (math-infinitep calc-low) (math-infinitep calc-high))
1892 (not step)
1893 (math-scan-for-limits expr))
1894 (and step (math-zerop step) (math-reject-arg step 'nonzerop))
1895 (let ((known (+ (if (Math-objectp calc-low) 1 0)
1896 (if (Math-objectp calc-high) 1 0)
1897 (if (or (null step) (Math-objectp step)) 1 0)))
1898 (count '(var inf var-inf))
1899 vec)
1900 (or (= known 2) ; handy optimization
1901 (equal calc-high '(var inf var-inf))
1902 (progn
1903 (setq count (math-div (math-sub calc-high calc-low) (or step 1)))
1904 (or (Math-objectp count)
1905 (setq count (math-simplify count)))
1906 (if (Math-messy-integerp count)
1907 (setq count (math-trunc count)))))
1908 (if (Math-negp count)
1909 (setq count -1))
1910 (if (integerp count)
1911 (let ((var-DUMMY nil)
1912 (vec math-tabulate-initial)
1913 (math-working-step-2 (1+ count))
1914 (math-working-step 0))
1915 (setq expr (math-evaluate-expr
1916 (math-expr-subst expr math-var '(var DUMMY var-DUMMY))))
1917 (while (>= count 0)
1918 (setq math-working-step (1+ math-working-step)
1919 var-DUMMY calc-low
1920 vec (cond ((eq math-tabulate-function 'calcFunc-sum)
1921 (math-add vec (math-evaluate-expr expr)))
1922 ((eq math-tabulate-function 'calcFunc-prod)
1923 (math-mul vec (math-evaluate-expr expr)))
1925 (cons (math-evaluate-expr expr) vec)))
1926 calc-low (math-add calc-low (or step 1))
1927 count (1- count)))
1928 (if math-tabulate-function
1930 (cons 'vec (nreverse vec))))
1931 (if (Math-integerp count)
1932 (calc-record-why 'fixnump calc-high)
1933 (if (Math-num-integerp calc-low)
1934 (if (Math-num-integerp calc-high)
1935 (calc-record-why 'integerp step)
1936 (calc-record-why 'integerp calc-high))
1937 (calc-record-why 'integerp calc-low)))
1938 (append (list (or math-tabulate-function 'calcFunc-table)
1939 expr math-var)
1940 (and (not (and (equal calc-low '(neg (var inf var-inf)))
1941 (equal calc-high '(var inf var-inf))))
1942 (list calc-low calc-high))
1943 (and step (list step))))))
1945 (defun math-scan-for-limits (x)
1946 (cond ((Math-primp x))
1947 ((and (eq (car x) 'calcFunc-subscr)
1948 (Math-vectorp (nth 1 x))
1949 (math-expr-contains (nth 2 x) math-var))
1950 (let* ((calc-next-why nil)
1951 (low-val (math-solve-for (nth 2 x) 1 math-var nil))
1952 (high-val (math-solve-for (nth 2 x) (1- (length (nth 1 x)))
1953 math-var nil))
1954 temp)
1955 (and low-val (math-realp low-val)
1956 high-val (math-realp high-val))
1957 (and (Math-lessp high-val low-val)
1958 (setq temp low-val low-val high-val high-val temp))
1959 (setq calc-low (math-max calc-low (math-ceiling low-val))
1960 calc-high (math-min calc-high (math-floor high-val)))))
1962 (while (setq x (cdr x))
1963 (math-scan-for-limits (car x))))))
1966 (defvar math-disable-sums nil)
1967 (defun calcFunc-sum (expr var &optional low high step)
1968 (if math-disable-sums (math-reject-arg))
1969 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2)))
1970 (math-sum-rec expr var low high step)))
1971 (math-disable-sums t))
1972 (math-normalize res)))
1974 (defun math-sum-rec (expr var &optional low high step)
1975 (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
1976 (and low (not high) (setq high low low 1))
1977 (let (t1 t2 val)
1978 (setq val
1979 (cond
1980 ((not (math-expr-contains expr var))
1981 (math-mul expr (math-add (math-div (math-sub high low) (or step 1))
1982 1)))
1983 ((and step (not (math-equal-int step 1)))
1984 (if (math-negp step)
1985 (math-sum-rec expr var high low (math-neg step))
1986 (let ((lo (math-simplify (math-div low step))))
1987 (if (math-known-num-integerp lo)
1988 (math-sum-rec (math-normalize
1989 (math-expr-subst expr var
1990 (math-mul step var)))
1991 var lo (math-simplify (math-div high step)))
1992 (math-sum-rec (math-normalize
1993 (math-expr-subst expr var
1994 (math-add (math-mul step var)
1995 low)))
1996 var 0
1997 (math-simplify (math-div (math-sub high low)
1998 step)))))))
1999 ((memq (setq t1 (math-compare low high)) '(0 1))
2000 (if (eq t1 0)
2001 (math-expr-subst expr var low)
2003 ((setq t1 (math-is-polynomial expr var 20))
2004 (let ((poly nil)
2005 (n 0))
2006 (while t1
2007 (setq poly (math-poly-mix poly 1
2008 (math-sum-integer-power n) (car t1))
2009 n (1+ n)
2010 t1 (cdr t1)))
2011 (setq n (math-build-polynomial-expr poly high))
2012 (if (= low 1)
2014 (math-sub n (math-build-polynomial-expr poly
2015 (math-sub low 1))))))
2016 ((and (memq (car expr) '(+ -))
2017 (setq t1 (math-sum-rec (nth 1 expr) var low high)
2018 t2 (math-sum-rec (nth 2 expr) var low high))
2019 (not (and (math-expr-calls t1 '(calcFunc-sum))
2020 (math-expr-calls t2 '(calcFunc-sum)))))
2021 (list (car expr) t1 t2))
2022 ((and (eq (car expr) '*)
2023 (setq t1 (math-sum-const-factors expr var)))
2024 (math-mul (car t1) (math-sum-rec (cdr t1) var low high)))
2025 ((and (eq (car expr) '*) (memq (car-safe (nth 1 expr)) '(+ -)))
2026 (math-sum-rec (math-add-or-sub (math-mul (nth 1 (nth 1 expr))
2027 (nth 2 expr))
2028 (math-mul (nth 2 (nth 1 expr))
2029 (nth 2 expr))
2030 nil (eq (car (nth 1 expr)) '-))
2031 var low high))
2032 ((and (eq (car expr) '*) (memq (car-safe (nth 2 expr)) '(+ -)))
2033 (math-sum-rec (math-add-or-sub (math-mul (nth 1 expr)
2034 (nth 1 (nth 2 expr)))
2035 (math-mul (nth 1 expr)
2036 (nth 2 (nth 2 expr)))
2037 nil (eq (car (nth 2 expr)) '-))
2038 var low high))
2039 ((and (eq (car expr) '/)
2040 (not (math-primp (nth 1 expr)))
2041 (setq t1 (math-sum-const-factors (nth 1 expr) var)))
2042 (math-mul (car t1)
2043 (math-sum-rec (math-div (cdr t1) (nth 2 expr))
2044 var low high)))
2045 ((and (eq (car expr) '/)
2046 (setq t1 (math-sum-const-factors (nth 2 expr) var)))
2047 (math-div (math-sum-rec (math-div (nth 1 expr) (cdr t1))
2048 var low high)
2049 (car t1)))
2050 ((eq (car expr) 'neg)
2051 (math-neg (math-sum-rec (nth 1 expr) var low high)))
2052 ((and (eq (car expr) '^)
2053 (not (math-expr-contains (nth 1 expr) var))
2054 (setq t1 (math-is-polynomial (nth 2 expr) var 1)))
2055 (let ((x (math-pow (nth 1 expr) (nth 1 t1))))
2056 (math-div (math-mul (math-sub (math-pow x (math-add 1 high))
2057 (math-pow x low))
2058 (math-pow (nth 1 expr) (car t1)))
2059 (math-sub x 1))))
2060 ((and (setq t1 (math-to-exponentials expr))
2061 (setq t1 (math-sum-rec t1 var low high))
2062 (not (math-expr-calls t1 '(calcFunc-sum))))
2063 (math-to-exps t1))
2064 ((memq (car expr) '(calcFunc-ln calcFunc-log10))
2065 (list (car expr) (calcFunc-prod (nth 1 expr) var low high)))
2066 ((and (eq (car expr) 'calcFunc-log)
2067 (= (length expr) 3)
2068 (not (math-expr-contains (nth 2 expr) var)))
2069 (list 'calcFunc-log
2070 (calcFunc-prod (nth 1 expr) var low high)
2071 (nth 2 expr)))))
2072 (if (equal val '(var nan var-nan)) (setq val nil))
2073 (or val
2074 (let* ((math-tabulate-initial 0)
2075 (math-tabulate-function 'calcFunc-sum))
2076 (calcFunc-table expr var low high)))))
2078 (defun calcFunc-asum (expr var low &optional high step no-mul-flag)
2079 (or high (setq high low low 1))
2080 (if (and step (not (math-equal-int step 1)))
2081 (if (math-negp step)
2082 (math-mul (math-pow -1 low)
2083 (calcFunc-asum expr var high low (math-neg step) t))
2084 (let ((lo (math-simplify (math-div low step))))
2085 (if (math-num-integerp lo)
2086 (calcFunc-asum (math-normalize
2087 (math-expr-subst expr var
2088 (math-mul step var)))
2089 var lo (math-simplify (math-div high step)))
2090 (calcFunc-asum (math-normalize
2091 (math-expr-subst expr var
2092 (math-add (math-mul step var)
2093 low)))
2094 var 0
2095 (math-simplify (math-div (math-sub high low)
2096 step))))))
2097 (math-mul (if no-mul-flag 1 (math-pow -1 low))
2098 (calcFunc-sum (math-mul (math-pow -1 var) expr) var low high))))
2100 (defun math-sum-const-factors (expr var)
2101 (let ((const nil)
2102 (not-const nil)
2103 (p expr))
2104 (while (eq (car-safe p) '*)
2105 (if (math-expr-contains (nth 1 p) var)
2106 (setq not-const (cons (nth 1 p) not-const))
2107 (setq const (cons (nth 1 p) const)))
2108 (setq p (nth 2 p)))
2109 (if (math-expr-contains p var)
2110 (setq not-const (cons p not-const))
2111 (setq const (cons p const)))
2112 (and const
2113 (cons (let ((temp (car const)))
2114 (while (setq const (cdr const))
2115 (setq temp (list '* (car const) temp)))
2116 temp)
2117 (let ((temp (or (car not-const) 1)))
2118 (while (setq not-const (cdr not-const))
2119 (setq temp (list '* (car not-const) temp)))
2120 temp)))))
2122 (defvar math-sum-int-pow-cache (list '(0 1)))
2123 ;; Following is from CRC Math Tables, 27th ed, pp. 52-53.
2124 (defun math-sum-integer-power (pow)
2125 (let ((calc-prefer-frac t)
2126 (n (length math-sum-int-pow-cache)))
2127 (while (<= n pow)
2128 (let* ((new (list 0 0))
2129 (lin new)
2130 (pp (cdr (nth (1- n) math-sum-int-pow-cache)))
2131 (p 2)
2132 (sum 0)
2134 (while pp
2135 (setq q (math-div (car pp) p)
2136 new (cons (math-mul q n) new)
2137 sum (math-add sum q)
2138 p (1+ p)
2139 pp (cdr pp)))
2140 (setcar lin (math-sub 1 (math-mul n sum)))
2141 (setq math-sum-int-pow-cache
2142 (nconc math-sum-int-pow-cache (list (nreverse new)))
2143 n (1+ n))))
2144 (nth pow math-sum-int-pow-cache)))
2146 (defun math-to-exponentials (expr)
2147 (and (consp expr)
2148 (= (length expr) 2)
2149 (let ((x (nth 1 expr))
2150 (pi (if calc-symbolic-mode '(var pi var-pi) (math-pi)))
2151 (i (if calc-symbolic-mode '(var i var-i) '(cplx 0 1))))
2152 (cond ((eq (car expr) 'calcFunc-exp)
2153 (list '^ '(var e var-e) x))
2154 ((eq (car expr) 'calcFunc-sin)
2155 (or (eq calc-angle-mode 'rad)
2156 (setq x (list '/ (list '* x pi) 180)))
2157 (list '/ (list '-
2158 (list '^ '(var e var-e) (list '* x i))
2159 (list '^ '(var e var-e)
2160 (list 'neg (list '* x i))))
2161 (list '* 2 i)))
2162 ((eq (car expr) 'calcFunc-cos)
2163 (or (eq calc-angle-mode 'rad)
2164 (setq x (list '/ (list '* x pi) 180)))
2165 (list '/ (list '+
2166 (list '^ '(var e var-e)
2167 (list '* x i))
2168 (list '^ '(var e var-e)
2169 (list 'neg (list '* x i))))
2171 ((eq (car expr) 'calcFunc-sinh)
2172 (list '/ (list '-
2173 (list '^ '(var e var-e) x)
2174 (list '^ '(var e var-e) (list 'neg x)))
2176 ((eq (car expr) 'calcFunc-cosh)
2177 (list '/ (list '+
2178 (list '^ '(var e var-e) x)
2179 (list '^ '(var e var-e) (list 'neg x)))
2181 (t nil)))))
2183 (defun math-to-exps (expr)
2184 (cond (calc-symbolic-mode expr)
2185 ((Math-primp expr)
2186 (if (equal expr '(var e var-e)) (math-e) expr))
2187 ((and (eq (car expr) '^)
2188 (equal (nth 1 expr) '(var e var-e)))
2189 (list 'calcFunc-exp (nth 2 expr)))
2191 (cons (car expr) (mapcar 'math-to-exps (cdr expr))))))
2194 (defvar math-disable-prods nil)
2195 (defun calcFunc-prod (expr var &optional low high step)
2196 (if math-disable-prods (math-reject-arg))
2197 (let* ((res (let* ((calc-internal-prec (+ calc-internal-prec 2)))
2198 (math-prod-rec expr var low high step)))
2199 (math-disable-prods t))
2200 (math-normalize res)))
2202 (defun math-prod-rec (expr var &optional low high step)
2203 (or low (setq low '(neg (var inf var-inf)) high '(var inf var-inf)))
2204 (and low (not high) (setq high '(var inf var-inf)))
2205 (let (t1 t2 t3 val)
2206 (setq val
2207 (cond
2208 ((not (math-expr-contains expr var))
2209 (math-pow expr (math-add (math-div (math-sub high low) (or step 1))
2210 1)))
2211 ((and step (not (math-equal-int step 1)))
2212 (if (math-negp step)
2213 (math-prod-rec expr var high low (math-neg step))
2214 (let ((lo (math-simplify (math-div low step))))
2215 (if (math-known-num-integerp lo)
2216 (math-prod-rec (math-normalize
2217 (math-expr-subst expr var
2218 (math-mul step var)))
2219 var lo (math-simplify (math-div high step)))
2220 (math-prod-rec (math-normalize
2221 (math-expr-subst expr var
2222 (math-add (math-mul step
2223 var)
2224 low)))
2225 var 0
2226 (math-simplify (math-div (math-sub high low)
2227 step)))))))
2228 ((and (memq (car expr) '(* /))
2229 (setq t1 (math-prod-rec (nth 1 expr) var low high)
2230 t2 (math-prod-rec (nth 2 expr) var low high))
2231 (not (and (math-expr-calls t1 '(calcFunc-prod))
2232 (math-expr-calls t2 '(calcFunc-prod)))))
2233 (list (car expr) t1 t2))
2234 ((and (eq (car expr) '^)
2235 (not (math-expr-contains (nth 2 expr) var)))
2236 (math-pow (math-prod-rec (nth 1 expr) var low high)
2237 (nth 2 expr)))
2238 ((and (eq (car expr) '^)
2239 (not (math-expr-contains (nth 1 expr) var)))
2240 (math-pow (nth 1 expr)
2241 (calcFunc-sum (nth 2 expr) var low high)))
2242 ((eq (car expr) 'sqrt)
2243 (math-normalize (list 'calcFunc-sqrt
2244 (list 'calcFunc-prod (nth 1 expr)
2245 var low high))))
2246 ((eq (car expr) 'neg)
2247 (math-mul (math-pow -1 (math-add (math-sub high low) 1))
2248 (math-prod-rec (nth 1 expr) var low high)))
2249 ((eq (car expr) 'calcFunc-exp)
2250 (list 'calcFunc-exp (calcFunc-sum (nth 1 expr) var low high)))
2251 ((and (setq t1 (math-is-polynomial expr var 1))
2252 (setq t2
2253 (cond
2254 ((or (and (math-equal-int (nth 1 t1) 1)
2255 (setq low (math-simplify
2256 (math-add low (car t1)))
2257 high (math-simplify
2258 (math-add high (car t1)))))
2259 (and (math-equal-int (nth 1 t1) -1)
2260 (setq t2 low
2261 low (math-simplify
2262 (math-sub (car t1) high))
2263 high (math-simplify
2264 (math-sub (car t1) t2)))))
2265 (if (or (math-zerop low) (math-zerop high))
2267 (if (and (or (math-negp low) (math-negp high))
2268 (or (math-num-integerp low)
2269 (math-num-integerp high)))
2270 (if (math-posp high)
2272 (math-mul (math-pow -1
2273 (math-add
2274 (math-add low high) 1))
2275 (list '/
2276 (list 'calcFunc-fact
2277 (math-neg low))
2278 (list 'calcFunc-fact
2279 (math-sub -1 high)))))
2280 (list '/
2281 (list 'calcFunc-fact high)
2282 (list 'calcFunc-fact (math-sub low 1))))))
2283 ((and (or (and (math-equal-int (nth 1 t1) 2)
2284 (setq t2 (math-simplify
2285 (math-add (math-mul low 2)
2286 (car t1)))
2287 t3 (math-simplify
2288 (math-add (math-mul high 2)
2289 (car t1)))))
2290 (and (math-equal-int (nth 1 t1) -2)
2291 (setq t2 (math-simplify
2292 (math-sub (car t1)
2293 (math-mul high 2)))
2294 t3 (math-simplify
2295 (math-sub (car t1)
2296 (math-mul low
2297 2))))))
2298 (or (math-integerp t2)
2299 (and (math-messy-integerp t2)
2300 (setq t2 (math-trunc t2)))
2301 (math-integerp t3)
2302 (and (math-messy-integerp t3)
2303 (setq t3 (math-trunc t3)))))
2304 (if (or (math-zerop t2) (math-zerop t3))
2306 (if (or (math-evenp t2) (math-evenp t3))
2307 (if (or (math-negp t2) (math-negp t3))
2308 (if (math-posp high)
2310 (list '/
2311 (list 'calcFunc-dfact
2312 (math-neg t2))
2313 (list 'calcFunc-dfact
2314 (math-sub -2 t3))))
2315 (list '/
2316 (list 'calcFunc-dfact t3)
2317 (list 'calcFunc-dfact
2318 (math-sub t2 2))))
2319 (if (math-negp t3)
2320 (list '*
2321 (list '^ -1
2322 (list '/ (list '- (list '- t2 t3)
2325 (list '/
2326 (list 'calcFunc-dfact
2327 (math-neg t2))
2328 (list 'calcFunc-dfact
2329 (math-sub -2 t3))))
2330 (if (math-posp t2)
2331 (list '/
2332 (list 'calcFunc-dfact t3)
2333 (list 'calcFunc-dfact
2334 (math-sub t2 2)))
2335 nil))))))))
2336 t2)))
2337 (if (equal val '(var nan var-nan)) (setq val nil))
2338 (or val
2339 (let* ((math-tabulate-initial 1)
2340 (math-tabulate-function 'calcFunc-prod))
2341 (calcFunc-table expr var low high)))))
2346 (defvar math-solve-ranges nil)
2347 (defvar math-solve-sign)
2348 ;;; Attempt to reduce math-solve-lhs = math-solve-rhs to
2349 ;;; math-solve-var = math-solve-rhs', where math-solve-var appears
2350 ;;; in math-solve-lhs but not in math-solve-rhs or math-solve-rhs';
2351 ;;; return math-solve-rhs'.
2352 ;;; Uses global values: math-solve-var, math-solve-full.
2353 (defvar math-solve-var)
2354 (defvar math-solve-full)
2356 ;; The variables math-solve-lhs, math-solve-rhs and math-try-solve-sign
2357 ;; are local to math-try-solve-for, but are used by math-try-solve-prod.
2358 ;; (math-solve-lhs and math-solve-rhs are is also local to
2359 ;; math-decompose-poly, but used by math-solve-poly-funny-powers.)
2360 (defvar math-solve-lhs)
2361 (defvar math-solve-rhs)
2362 (defvar math-try-solve-sign)
2364 (defun math-try-solve-for
2365 (math-solve-lhs math-solve-rhs &optional math-try-solve-sign no-poly)
2366 (let (math-t1 math-t2 math-t3)
2367 (cond ((equal math-solve-lhs math-solve-var)
2368 (setq math-solve-sign math-try-solve-sign)
2369 (if (eq math-solve-full 'all)
2370 (let ((vec (list 'vec (math-evaluate-expr math-solve-rhs)))
2371 newvec var p)
2372 (while math-solve-ranges
2373 (setq p (car math-solve-ranges)
2374 var (car p)
2375 newvec (list 'vec))
2376 (while (setq p (cdr p))
2377 (setq newvec (nconc newvec
2378 (cdr (math-expr-subst
2379 vec var (car p))))))
2380 (setq vec newvec
2381 math-solve-ranges (cdr math-solve-ranges)))
2382 (math-normalize vec))
2383 math-solve-rhs))
2384 ((Math-primp math-solve-lhs)
2385 nil)
2386 ((and (eq (car math-solve-lhs) '-)
2387 (eq (car-safe (nth 1 math-solve-lhs)) (car-safe (nth 2 math-solve-lhs)))
2388 (Math-zerop math-solve-rhs)
2389 (= (length (nth 1 math-solve-lhs)) 2)
2390 (= (length (nth 2 math-solve-lhs)) 2)
2391 (setq math-t1 (get (car (nth 1 math-solve-lhs)) 'math-inverse))
2392 (setq math-t2 (funcall math-t1 '(var SOLVEDUM SOLVEDUM)))
2393 (eq (math-expr-contains-count math-t2 '(var SOLVEDUM SOLVEDUM)) 1)
2394 (setq math-t3 (math-solve-above-dummy math-t2))
2395 (setq math-t1 (math-try-solve-for
2396 (math-sub (nth 1 (nth 1 math-solve-lhs))
2397 (math-expr-subst
2398 math-t2 math-t3
2399 (nth 1 (nth 2 math-solve-lhs))))
2400 0)))
2401 math-t1)
2402 ((eq (car math-solve-lhs) 'neg)
2403 (math-try-solve-for (nth 1 math-solve-lhs) (math-neg math-solve-rhs)
2404 (and math-try-solve-sign (- math-try-solve-sign))))
2405 ((and (not (eq math-solve-full 't)) (math-try-solve-prod)))
2406 ((and (not no-poly)
2407 (setq math-t2
2408 (math-decompose-poly math-solve-lhs
2409 math-solve-var 15 math-solve-rhs)))
2410 (setq math-t1 (cdr (nth 1 math-t2))
2411 math-t1 (let ((math-solve-ranges math-solve-ranges))
2412 (cond ((= (length math-t1) 5)
2413 (apply 'math-solve-quartic (car math-t2) math-t1))
2414 ((= (length math-t1) 4)
2415 (apply 'math-solve-cubic (car math-t2) math-t1))
2416 ((= (length math-t1) 3)
2417 (apply 'math-solve-quadratic (car math-t2) math-t1))
2418 ((= (length math-t1) 2)
2419 (apply 'math-solve-linear
2420 (car math-t2) math-try-solve-sign math-t1))
2421 (math-solve-full
2422 (math-poly-all-roots (car math-t2) math-t1))
2423 (calc-symbolic-mode nil)
2425 (math-try-solve-for
2426 (car math-t2)
2427 (math-poly-any-root (reverse math-t1) 0 t)
2428 nil t)))))
2429 (if math-t1
2430 (if (eq (nth 2 math-t2) 1)
2431 math-t1
2432 (math-solve-prod math-t1 (math-try-solve-for (nth 2 math-t2) 0 nil t)))
2433 (calc-record-why "*Unable to find a symbolic solution")
2434 nil))
2435 ((and (math-solve-find-root-term math-solve-lhs nil)
2436 (eq (math-expr-contains-count math-solve-lhs math-t1) 1)) ; just in case
2437 (math-try-solve-for (math-simplify
2438 (math-sub (if (or math-t3 (math-evenp math-t2))
2439 (math-pow math-t1 math-t2)
2440 (math-neg (math-pow math-t1 math-t2)))
2441 (math-expand-power
2442 (math-sub (math-normalize
2443 (math-expr-subst
2444 math-solve-lhs math-t1 0))
2445 math-solve-rhs)
2446 math-t2 math-solve-var)))
2448 ((eq (car math-solve-lhs) '+)
2449 (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var))
2450 (math-try-solve-for (nth 2 math-solve-lhs)
2451 (math-sub math-solve-rhs (nth 1 math-solve-lhs))
2452 math-try-solve-sign))
2453 ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var))
2454 (math-try-solve-for (nth 1 math-solve-lhs)
2455 (math-sub math-solve-rhs (nth 2 math-solve-lhs))
2456 math-try-solve-sign))))
2457 ((eq (car math-solve-lhs) 'calcFunc-eq)
2458 (math-try-solve-for (math-sub (nth 1 math-solve-lhs) (nth 2 math-solve-lhs))
2459 math-solve-rhs math-try-solve-sign no-poly))
2460 ((eq (car math-solve-lhs) '-)
2461 (cond ((or (and (eq (car-safe (nth 1 math-solve-lhs)) 'calcFunc-sin)
2462 (eq (car-safe (nth 2 math-solve-lhs)) 'calcFunc-cos))
2463 (and (eq (car-safe (nth 1 math-solve-lhs)) 'calcFunc-cos)
2464 (eq (car-safe (nth 2 math-solve-lhs)) 'calcFunc-sin)))
2465 (math-try-solve-for (math-sub (nth 1 math-solve-lhs)
2466 (list (car (nth 1 math-solve-lhs))
2467 (math-sub
2468 (math-quarter-circle t)
2469 (nth 1 (nth 2 math-solve-lhs)))))
2470 math-solve-rhs))
2471 ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var))
2472 (math-try-solve-for (nth 2 math-solve-lhs)
2473 (math-sub (nth 1 math-solve-lhs) math-solve-rhs)
2474 (and math-try-solve-sign
2475 (- math-try-solve-sign))))
2476 ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var))
2477 (math-try-solve-for (nth 1 math-solve-lhs)
2478 (math-add math-solve-rhs (nth 2 math-solve-lhs))
2479 math-try-solve-sign))))
2480 ((and (eq math-solve-full 't) (math-try-solve-prod)))
2481 ((and (eq (car math-solve-lhs) '%)
2482 (not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var)))
2483 (math-try-solve-for (nth 1 math-solve-lhs) (math-add math-solve-rhs
2484 (math-solve-get-int
2485 (nth 2 math-solve-lhs)))))
2486 ((eq (car math-solve-lhs) 'calcFunc-log)
2487 (cond ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var))
2488 (math-try-solve-for (nth 1 math-solve-lhs)
2489 (math-pow (nth 2 math-solve-lhs) math-solve-rhs)))
2490 ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var))
2491 (math-try-solve-for (nth 2 math-solve-lhs) (math-pow
2492 (nth 1 math-solve-lhs)
2493 (math-div 1 math-solve-rhs))))))
2494 ((and (= (length math-solve-lhs) 2)
2495 (symbolp (car math-solve-lhs))
2496 (setq math-t1 (get (car math-solve-lhs) 'math-inverse))
2497 (setq math-t2 (funcall math-t1 math-solve-rhs)))
2498 (setq math-t1 (get (car math-solve-lhs) 'math-inverse-sign))
2499 (math-try-solve-for (nth 1 math-solve-lhs) (math-normalize math-t2)
2500 (and math-try-solve-sign math-t1
2501 (if (integerp math-t1)
2502 (* math-t1 math-try-solve-sign)
2503 (funcall math-t1 math-solve-lhs
2504 math-try-solve-sign)))))
2505 ((and (symbolp (car math-solve-lhs))
2506 (setq math-t1 (get (car math-solve-lhs) 'math-inverse-n))
2507 (setq math-t2 (funcall math-t1 math-solve-lhs math-solve-rhs)))
2508 math-t2)
2509 ((setq math-t1 (math-expand-formula math-solve-lhs))
2510 (math-try-solve-for math-t1 math-solve-rhs math-try-solve-sign))
2512 (calc-record-why "*No inverse known" math-solve-lhs)
2513 nil))))
2516 (defun math-try-solve-prod ()
2517 (cond ((eq (car math-solve-lhs) '*)
2518 (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var))
2519 (math-try-solve-for (nth 2 math-solve-lhs)
2520 (math-div math-solve-rhs (nth 1 math-solve-lhs))
2521 (math-solve-sign math-try-solve-sign
2522 (nth 1 math-solve-lhs))))
2523 ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var))
2524 (math-try-solve-for (nth 1 math-solve-lhs)
2525 (math-div math-solve-rhs (nth 2 math-solve-lhs))
2526 (math-solve-sign math-try-solve-sign
2527 (nth 2 math-solve-lhs))))
2528 ((Math-zerop math-solve-rhs)
2529 (math-solve-prod (let ((math-solve-ranges math-solve-ranges))
2530 (math-try-solve-for (nth 2 math-solve-lhs) 0))
2531 (math-try-solve-for (nth 1 math-solve-lhs) 0)))))
2532 ((eq (car math-solve-lhs) '/)
2533 (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var))
2534 (math-try-solve-for (nth 2 math-solve-lhs)
2535 (math-div (nth 1 math-solve-lhs) math-solve-rhs)
2536 (math-solve-sign math-try-solve-sign
2537 (nth 1 math-solve-lhs))))
2538 ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var))
2539 (math-try-solve-for (nth 1 math-solve-lhs)
2540 (math-mul math-solve-rhs (nth 2 math-solve-lhs))
2541 (math-solve-sign math-try-solve-sign
2542 (nth 2 math-solve-lhs))))
2543 ((setq math-t1 (math-try-solve-for (math-sub (nth 1 math-solve-lhs)
2544 (math-mul (nth 2 math-solve-lhs)
2545 math-solve-rhs))
2547 math-t1)))
2548 ((eq (car math-solve-lhs) '^)
2549 (cond ((not (math-expr-contains (nth 1 math-solve-lhs) math-solve-var))
2550 (math-try-solve-for
2551 (nth 2 math-solve-lhs)
2552 (math-add (math-normalize
2553 (list 'calcFunc-log math-solve-rhs (nth 1 math-solve-lhs)))
2554 (math-div
2555 (math-mul 2
2556 (math-mul '(var pi var-pi)
2557 (math-solve-get-int
2558 '(var i var-i))))
2559 (math-normalize
2560 (list 'calcFunc-ln (nth 1 math-solve-lhs)))))))
2561 ((not (math-expr-contains (nth 2 math-solve-lhs) math-solve-var))
2562 (cond ((and (integerp (nth 2 math-solve-lhs))
2563 (>= (nth 2 math-solve-lhs) 2)
2564 (setq math-t1 (math-integer-log2 (nth 2 math-solve-lhs))))
2565 (setq math-t2 math-solve-rhs)
2566 (if (and (eq math-solve-full t)
2567 (math-known-realp (nth 1 math-solve-lhs)))
2568 (progn
2569 (while (>= (setq math-t1 (1- math-t1)) 0)
2570 (setq math-t2 (list 'calcFunc-sqrt math-t2)))
2571 (setq math-t2 (math-solve-get-sign math-t2)))
2572 (while (>= (setq math-t1 (1- math-t1)) 0)
2573 (setq math-t2 (math-solve-get-sign
2574 (math-normalize
2575 (list 'calcFunc-sqrt math-t2))))))
2576 (math-try-solve-for
2577 (nth 1 math-solve-lhs)
2578 (math-normalize math-t2)))
2579 ((math-looks-negp (nth 2 math-solve-lhs))
2580 (math-try-solve-for
2581 (list '^ (nth 1 math-solve-lhs)
2582 (math-neg (nth 2 math-solve-lhs)))
2583 (math-div 1 math-solve-rhs)))
2584 ((and (eq math-solve-full t)
2585 (Math-integerp (nth 2 math-solve-lhs))
2586 (math-known-realp (nth 1 math-solve-lhs)))
2587 (setq math-t1 (math-normalize
2588 (list 'calcFunc-nroot math-solve-rhs
2589 (nth 2 math-solve-lhs))))
2590 (if (math-evenp (nth 2 math-solve-lhs))
2591 (setq math-t1 (math-solve-get-sign math-t1)))
2592 (math-try-solve-for
2593 (nth 1 math-solve-lhs) math-t1
2594 (and math-try-solve-sign
2595 (math-oddp (nth 2 math-solve-lhs))
2596 (math-solve-sign math-try-solve-sign
2597 (nth 2 math-solve-lhs)))))
2598 (t (math-try-solve-for
2599 (nth 1 math-solve-lhs)
2600 (math-mul
2601 (math-normalize
2602 (list 'calcFunc-exp
2603 (if (Math-realp (nth 2 math-solve-lhs))
2604 (math-div (math-mul
2605 '(var pi var-pi)
2606 (math-solve-get-int
2607 '(var i var-i)
2608 (and (integerp (nth 2 math-solve-lhs))
2609 (math-abs
2610 (nth 2 math-solve-lhs)))))
2611 (math-div (nth 2 math-solve-lhs) 2))
2612 (math-div (math-mul
2614 (math-mul
2615 '(var pi var-pi)
2616 (math-solve-get-int
2617 '(var i var-i)
2618 (and (integerp (nth 2 math-solve-lhs))
2619 (math-abs
2620 (nth 2 math-solve-lhs))))))
2621 (nth 2 math-solve-lhs)))))
2622 (math-normalize
2623 (list 'calcFunc-nroot
2624 math-solve-rhs
2625 (nth 2 math-solve-lhs))))
2626 (and math-try-solve-sign
2627 (math-oddp (nth 2 math-solve-lhs))
2628 (math-solve-sign math-try-solve-sign
2629 (nth 2 math-solve-lhs)))))))))
2630 (t nil)))
2632 (defun math-solve-prod (lsoln rsoln)
2633 (cond ((null lsoln)
2634 rsoln)
2635 ((null rsoln)
2636 lsoln)
2637 ((eq math-solve-full 'all)
2638 (cons 'vec (append (cdr lsoln) (cdr rsoln))))
2639 (math-solve-full
2640 (list 'calcFunc-if
2641 (list 'calcFunc-gt (math-solve-get-sign 1) 0)
2642 lsoln
2643 rsoln))
2644 (t lsoln)))
2646 ;;; This deals with negative, fractional, and symbolic powers of "x".
2647 ;; The variable math-solve-b is local to math-decompose-poly,
2648 ;; but is used by math-solve-poly-funny-powers.
2649 (defvar math-solve-b)
2651 (defun math-solve-poly-funny-powers (sub-rhs) ; uses "t1", "t2"
2652 (setq math-t1 math-solve-lhs)
2653 (let ((pp math-poly-neg-powers)
2654 fac)
2655 (while pp
2656 (setq fac (math-pow (car pp) (or math-poly-mult-powers 1))
2657 math-t1 (math-mul math-t1 fac)
2658 math-solve-rhs (math-mul math-solve-rhs fac)
2659 pp (cdr pp))))
2660 (if sub-rhs (setq math-t1 (math-sub math-t1 math-solve-rhs)))
2661 (let ((math-poly-neg-powers nil))
2662 (setq math-t2 (math-mul (or math-poly-mult-powers 1)
2663 (let ((calc-prefer-frac t))
2664 (math-div 1 math-poly-frac-powers)))
2665 math-t1 (math-is-polynomial
2666 (math-simplify (calcFunc-expand math-t1)) math-solve-b 50))))
2668 ;;; This converts "a x^8 + b x^5 + c x^2" to "(a (x^3)^2 + b (x^3) + c) * x^2".
2669 (defun math-solve-crunch-poly (max-degree) ; uses "t1", "t3"
2670 (let ((count 0))
2671 (while (and math-t1 (Math-zerop (car math-t1)))
2672 (setq math-t1 (cdr math-t1)
2673 count (1+ count)))
2674 (and math-t1
2675 (let* ((degree (1- (length math-t1)))
2676 (scale degree))
2677 (while (and (> scale 1) (= (car math-t3) 1))
2678 (and (= (% degree scale) 0)
2679 (let ((p math-t1)
2680 (n 0)
2681 (new-t1 nil)
2682 (okay t))
2683 (while (and p okay)
2684 (if (= (% n scale) 0)
2685 (setq new-t1 (nconc new-t1 (list (car p))))
2686 (or (Math-zerop (car p))
2687 (setq okay nil)))
2688 (setq p (cdr p)
2689 n (1+ n)))
2690 (if okay
2691 (setq math-t3 (cons scale (cdr math-t3))
2692 math-t1 new-t1))))
2693 (setq scale (1- scale)))
2694 (setq math-t3 (list (math-mul (car math-t3) math-t2)
2695 (math-mul count math-t2)))
2696 (<= (1- (length math-t1)) max-degree)))))
2698 (defun calcFunc-poly (expr var &optional degree)
2699 (if degree
2700 (or (natnump degree) (math-reject-arg degree 'fixnatnump))
2701 (setq degree 50))
2702 (let ((p (math-is-polynomial expr var degree 'gen)))
2703 (if p
2704 (if (equal p '(0))
2705 (list 'vec)
2706 (cons 'vec p))
2707 (math-reject-arg expr "Expected a polynomial"))))
2709 (defun calcFunc-gpoly (expr var &optional degree)
2710 (if degree
2711 (or (natnump degree) (math-reject-arg degree 'fixnatnump))
2712 (setq degree 50))
2713 (let* ((math-poly-base-variable var)
2714 (d (math-decompose-poly expr var degree nil)))
2715 (if d
2716 (cons 'vec d)
2717 (math-reject-arg expr "Expected a polynomial"))))
2719 (defun math-decompose-poly (math-solve-lhs math-solve-var degree sub-rhs)
2720 (let ((math-solve-rhs (or sub-rhs 1))
2721 math-t1 math-t2 math-t3)
2722 (setq math-t2 (math-polynomial-base
2723 math-solve-lhs
2724 (function
2725 (lambda (math-solve-b)
2726 (let ((math-poly-neg-powers '(1))
2727 (math-poly-mult-powers nil)
2728 (math-poly-frac-powers 1)
2729 (math-poly-exp-base t))
2730 (and (not (equal math-solve-b math-solve-lhs))
2731 (or (not (memq (car-safe math-solve-b) '(+ -))) sub-rhs)
2732 (setq math-t3 '(1 0) math-t2 1
2733 math-t1 (math-is-polynomial math-solve-lhs
2734 math-solve-b 50))
2735 (if (and (equal math-poly-neg-powers '(1))
2736 (memq math-poly-mult-powers '(nil 1))
2737 (eq math-poly-frac-powers 1)
2738 sub-rhs)
2739 (setq math-t1 (cons (math-sub (car math-t1) math-solve-rhs)
2740 (cdr math-t1)))
2741 (math-solve-poly-funny-powers sub-rhs))
2742 (math-solve-crunch-poly degree)
2743 (or (math-expr-contains math-solve-b math-solve-var)
2744 (math-expr-contains (car math-t3) math-solve-var))))))))
2745 (if math-t2
2746 (list (math-pow math-t2 (car math-t3))
2747 (cons 'vec math-t1)
2748 (if sub-rhs
2749 (math-pow math-t2 (nth 1 math-t3))
2750 (math-div (math-pow math-t2 (nth 1 math-t3)) math-solve-rhs))))))
2752 (defun math-solve-linear (var sign b a)
2753 (math-try-solve-for var
2754 (math-div (math-neg b) a)
2755 (math-solve-sign sign a)
2758 (defun math-solve-quadratic (var c b a)
2759 (math-try-solve-for
2761 (if (math-looks-evenp b)
2762 (let ((halfb (math-div b 2)))
2763 (math-div
2764 (math-add
2765 (math-neg halfb)
2766 (math-solve-get-sign
2767 (math-normalize
2768 (list 'calcFunc-sqrt
2769 (math-add (math-sqr halfb)
2770 (math-mul (math-neg c) a))))))
2772 (math-div
2773 (math-add
2774 (math-neg b)
2775 (math-solve-get-sign
2776 (math-normalize
2777 (list 'calcFunc-sqrt
2778 (math-add (math-sqr b)
2779 (math-mul 4 (math-mul (math-neg c) a)))))))
2780 (math-mul 2 a)))
2781 nil t))
2783 (defun math-solve-cubic (var d c b a)
2784 (let* ((p (math-div b a))
2785 (q (math-div c a))
2786 (r (math-div d a))
2787 (psqr (math-sqr p))
2788 (aa (math-sub q (math-div psqr 3)))
2789 (bb (math-add r
2790 (math-div (math-sub (math-mul 2 (math-mul psqr p))
2791 (math-mul 9 (math-mul p q)))
2792 27)))
2794 (if (Math-zerop aa)
2795 (math-try-solve-for (math-pow (math-add var (math-div p 3)) 3)
2796 (math-neg bb) nil t)
2797 (if (Math-zerop bb)
2798 (math-try-solve-for
2799 (math-mul (math-add var (math-div p 3))
2800 (math-add (math-sqr (math-add var (math-div p 3)))
2801 aa))
2802 0 nil t)
2803 (setq m (math-mul 2 (list 'calcFunc-sqrt (math-div aa -3))))
2804 (math-try-solve-for
2806 (math-sub
2807 (math-normalize
2808 (math-mul
2810 (list 'calcFunc-cos
2811 (math-div
2812 (math-sub (list 'calcFunc-arccos
2813 (math-div (math-mul 3 bb)
2814 (math-mul aa m)))
2815 (math-mul 2
2816 (math-mul
2817 (math-add 1 (math-solve-get-int
2818 1 3))
2819 (math-half-circle
2820 calc-symbolic-mode))))
2821 3))))
2822 (math-div p 3))
2823 nil t)))))
2825 (defun math-solve-quartic (var d c b a aa)
2826 (setq a (math-div a aa))
2827 (setq b (math-div b aa))
2828 (setq c (math-div c aa))
2829 (setq d (math-div d aa))
2830 (math-try-solve-for
2832 (let* ((asqr (math-sqr a))
2833 (asqr4 (math-div asqr 4))
2834 (y (let ((math-solve-full nil)
2835 calc-next-why)
2836 (math-solve-cubic math-solve-var
2837 (math-sub (math-sub
2838 (math-mul 4 (math-mul b d))
2839 (math-mul asqr d))
2840 (math-sqr c))
2841 (math-sub (math-mul a c)
2842 (math-mul 4 d))
2843 (math-neg b)
2844 1)))
2845 (rsqr (math-add (math-sub asqr4 b) y))
2846 (r (list 'calcFunc-sqrt rsqr))
2847 (sign1 (math-solve-get-sign 1))
2848 (de (list 'calcFunc-sqrt
2849 (math-add
2850 (math-sub (math-mul 3 asqr4)
2851 (math-mul 2 b))
2852 (if (Math-zerop rsqr)
2853 (math-mul
2855 (math-mul sign1
2856 (list 'calcFunc-sqrt
2857 (math-sub (math-sqr y)
2858 (math-mul 4 d)))))
2859 (math-sub
2860 (math-mul sign1
2861 (math-div
2862 (math-sub (math-sub
2863 (math-mul 4 (math-mul a b))
2864 (math-mul 8 c))
2865 (math-mul asqr a))
2866 (math-mul 4 r)))
2867 rsqr))))))
2868 (math-normalize
2869 (math-sub (math-add (math-mul sign1 (math-div r 2))
2870 (math-solve-get-sign (math-div de 2)))
2871 (math-div a 4))))
2872 nil t))
2874 (defvar math-symbolic-solve nil)
2875 (defvar math-int-coefs nil)
2877 ;; The variable math-int-threshold is local to math-poly-all-roots,
2878 ;; but is used by math-poly-newton-root.
2879 (defvar math-int-threshold)
2880 ;; The variables math-int-scale, math-int-factors and math-double-roots
2881 ;; are local to math-poly-all-roots, but are used by math-poly-integer-root.
2882 (defvar math-int-scale)
2883 (defvar math-int-factors)
2884 (defvar math-double-roots)
2886 (defun math-poly-all-roots (var p &optional math-factoring)
2887 (catch 'ouch
2888 (let* ((math-symbolic-solve calc-symbolic-mode)
2889 (roots nil)
2890 (deg (1- (length p)))
2891 (orig-p (reverse p))
2892 (math-int-coefs nil)
2893 (math-int-scale nil)
2894 (math-double-roots nil)
2895 (math-int-factors nil)
2896 (math-int-threshold nil)
2897 (pp p))
2898 ;; If rational coefficients, look for exact rational factors.
2899 (while (and pp (Math-ratp (car pp)))
2900 (setq pp (cdr pp)))
2901 (if pp
2902 (if (or math-factoring math-symbolic-solve)
2903 (throw 'ouch nil))
2904 (let ((lead (car orig-p))
2905 (calc-prefer-frac t)
2906 (scale (apply 'math-lcm-denoms p)))
2907 (setq math-int-scale (math-abs (math-mul scale lead))
2908 math-int-threshold (math-div '(float 5 -2) math-int-scale)
2909 math-int-coefs (cdr (math-div (cons 'vec orig-p) lead)))))
2910 (if (> deg 4)
2911 (let ((calc-prefer-frac nil)
2912 (calc-symbolic-mode nil)
2913 (pp p)
2914 (def-p (copy-sequence orig-p)))
2915 (while pp
2916 (if (Math-numberp (car pp))
2917 (setq pp (cdr pp))
2918 (throw 'ouch nil)))
2919 (while (> deg (if math-symbolic-solve 2 4))
2920 (let* ((x (math-poly-any-root def-p '(float 0 0) nil))
2921 b c pp)
2922 (if (and (eq (car-safe x) 'cplx)
2923 (math-nearly-zerop (nth 2 x) (nth 1 x)))
2924 (setq x (calcFunc-re x)))
2925 (or math-factoring
2926 (setq roots (cons x roots)))
2927 (or (math-numberp x)
2928 (setq x (math-evaluate-expr x)))
2929 (setq pp def-p
2930 b (car def-p))
2931 (while (setq pp (cdr pp))
2932 (setq c (car pp))
2933 (setcar pp b)
2934 (setq b (math-add (math-mul x b) c)))
2935 (setq def-p (cdr def-p)
2936 deg (1- deg))))
2937 (setq p (reverse def-p))))
2938 (if (> deg 1)
2939 (let ((math-solve-var '(var DUMMY var-DUMMY))
2940 (math-solve-sign nil)
2941 (math-solve-ranges nil)
2942 (math-solve-full 'all))
2943 (if (= (length p) (length math-int-coefs))
2944 (setq p (reverse math-int-coefs)))
2945 (setq roots (append (cdr (apply (cond ((= deg 2)
2946 'math-solve-quadratic)
2947 ((= deg 3)
2948 'math-solve-cubic)
2950 'math-solve-quartic))
2951 math-solve-var p))
2952 roots)))
2953 (if (> deg 0)
2954 (setq roots (cons (math-div (math-neg (car p)) (nth 1 p))
2955 roots))))
2956 (if math-factoring
2957 (progn
2958 (while roots
2959 (math-poly-integer-root (car roots))
2960 (setq roots (cdr roots)))
2961 (list math-int-factors (nreverse math-int-coefs) math-int-scale))
2962 (let ((vec nil) res)
2963 (while roots
2964 (let ((root (car roots))
2965 (math-solve-full (and math-solve-full 'all)))
2966 (if (math-floatp root)
2967 (setq root (math-poly-any-root orig-p root t)))
2968 (setq vec (append vec
2969 (cdr (or (math-try-solve-for var root nil t)
2970 (throw 'ouch nil))))))
2971 (setq roots (cdr roots)))
2972 (setq vec (cons 'vec (nreverse vec)))
2973 (if math-symbolic-solve
2974 (setq vec (math-normalize vec)))
2975 (if (eq math-solve-full t)
2976 (list 'calcFunc-subscr
2978 (math-solve-get-int 1 (1- (length orig-p)) 1))
2979 vec))))))
2981 (defun math-lcm-denoms (&rest fracs)
2982 (let ((den 1))
2983 (while fracs
2984 (if (eq (car-safe (car fracs)) 'frac)
2985 (setq den (calcFunc-lcm den (nth 2 (car fracs)))))
2986 (setq fracs (cdr fracs)))
2987 den))
2989 (defun math-poly-any-root (p x polish) ; p is a reverse poly coeff list
2990 (let* ((newt (if (math-zerop x)
2991 (math-poly-newton-root
2992 p '(cplx (float 123 -6) (float 1 -4)) 4)
2993 (math-poly-newton-root p x 4)))
2994 (res (if (math-zerop (cdr newt))
2995 (car newt)
2996 (if (and (math-lessp (cdr newt) '(float 1 -3)) (not polish))
2997 (setq newt (math-poly-newton-root p (car newt) 30)))
2998 (if (math-zerop (cdr newt))
2999 (car newt)
3000 (math-poly-laguerre-root p x polish)))))
3001 (and math-symbolic-solve (math-floatp res)
3002 (throw 'ouch nil))
3003 res))
3005 (defun math-poly-newton-root (p x iters)
3006 (let* ((calc-prefer-frac nil)
3007 (calc-symbolic-mode nil)
3008 (try-integer math-int-coefs)
3009 (dx x) b d)
3010 (while (and (> (setq iters (1- iters)) 0)
3011 (let ((pp p))
3012 (math-working "newton" x)
3013 (setq b (car p)
3014 d 0)
3015 (while (setq pp (cdr pp))
3016 (setq d (math-add (math-mul x d) b)
3017 b (math-add (math-mul x b) (car pp))))
3018 (not (math-zerop d)))
3019 (progn
3020 (setq dx (math-div b d)
3021 x (math-sub x dx))
3022 (if try-integer
3023 (let ((adx (math-abs-approx dx)))
3024 (and (math-lessp adx math-int-threshold)
3025 (let ((iroot (math-poly-integer-root x)))
3026 (if iroot
3027 (setq x iroot dx 0)
3028 (setq try-integer nil))))))
3029 (or (not (or (eq dx 0)
3030 (math-nearly-zerop dx (math-abs-approx x))))
3031 (progn (setq dx 0) nil)))))
3032 (cons x (if (math-zerop x)
3033 1 (math-div (math-abs-approx dx) (math-abs-approx x))))))
3035 (defun math-poly-integer-root (x)
3036 (and (math-lessp (calcFunc-xpon (math-abs-approx x)) calc-internal-prec)
3037 math-int-coefs
3038 (let* ((calc-prefer-frac t)
3039 (xre (calcFunc-re x))
3040 (xim (calcFunc-im x))
3041 (xresq (math-sqr xre))
3042 (ximsq (math-sqr xim)))
3043 (if (math-lessp ximsq (calcFunc-scf xresq -1))
3044 ;; Look for linear factor
3045 (let* ((rnd (math-div (math-round (math-mul xre math-int-scale))
3046 math-int-scale))
3047 (icp math-int-coefs)
3048 (rem (car icp))
3049 (newcoef nil))
3050 (while (setq icp (cdr icp))
3051 (setq newcoef (cons rem newcoef)
3052 rem (math-add (car icp)
3053 (math-mul rem rnd))))
3054 (and (math-zerop rem)
3055 (progn
3056 (setq math-int-coefs (nreverse newcoef)
3057 math-int-factors (cons (list (math-neg rnd))
3058 math-int-factors))
3059 rnd)))
3060 ;; Look for irreducible quadratic factor
3061 (let* ((rnd1 (math-div (math-round
3062 (math-mul xre (math-mul -2 math-int-scale)))
3063 math-int-scale))
3064 (sqscale (math-sqr math-int-scale))
3065 (rnd0 (math-div (math-round (math-mul (math-add xresq ximsq)
3066 sqscale))
3067 sqscale))
3068 (rem1 (car math-int-coefs))
3069 (icp (cdr math-int-coefs))
3070 (rem0 (car icp))
3071 (newcoef nil)
3072 (found (assoc (list rnd0 rnd1 (math-posp xim))
3073 math-double-roots))
3074 this)
3075 (if found
3076 (setq math-double-roots (delq found math-double-roots)
3077 rem0 0 rem1 0)
3078 (while (setq icp (cdr icp))
3079 (setq this rem1
3080 newcoef (cons rem1 newcoef)
3081 rem1 (math-sub rem0 (math-mul this rnd1))
3082 rem0 (math-sub (car icp) (math-mul this rnd0)))))
3083 (and (math-zerop rem0)
3084 (math-zerop rem1)
3085 (let ((aa (math-div rnd1 -2)))
3086 (or found (setq math-int-coefs (reverse newcoef)
3087 math-double-roots (cons (list
3088 (list
3089 rnd0 rnd1
3090 (math-negp xim)))
3091 math-double-roots)
3092 math-int-factors (cons (cons rnd0 rnd1)
3093 math-int-factors)))
3094 (math-add aa
3095 (let ((calc-symbolic-mode math-symbolic-solve))
3096 (math-mul (math-sqrt (math-sub (math-sqr aa)
3097 rnd0))
3098 (if (math-negp xim) -1 1)))))))))))
3100 ;;; The following routine is from Numerical Recipes, section 9.5.
3101 (defun math-poly-laguerre-root (p x polish)
3102 (let* ((calc-prefer-frac nil)
3103 (calc-symbolic-mode nil)
3104 (iters 0)
3105 (m (1- (length p)))
3106 (try-newt (not polish))
3107 (tried-newt nil)
3108 b d f x1 dx dxold)
3109 (while
3110 (and (or (< (setq iters (1+ iters)) 50)
3111 (math-reject-arg x "*Laguerre's method failed to converge"))
3112 (let ((err (math-abs-approx (car p)))
3113 (abx (math-abs-approx x))
3114 (pp p))
3115 (setq b (car p)
3116 d 0 f 0)
3117 (while (setq pp (cdr pp))
3118 (setq f (math-add (math-mul x f) d)
3119 d (math-add (math-mul x d) b)
3120 b (math-add (math-mul x b) (car pp))
3121 err (math-add (math-abs-approx b) (math-mul abx err))))
3122 (math-lessp (calcFunc-scf err (- -2 calc-internal-prec))
3123 (math-abs-approx b)))
3124 (or (not (math-zerop d))
3125 (not (math-zerop f))
3126 (progn
3127 (setq x (math-pow (math-neg b) (list 'frac 1 m)))
3128 nil))
3129 (let* ((g (math-div d b))
3130 (g2 (math-sqr g))
3131 (h (math-sub g2 (math-mul 2 (math-div f b))))
3132 (sq (math-sqrt
3133 (math-mul (1- m) (math-sub (math-mul m h) g2))))
3134 (gp (math-add g sq))
3135 (gm (math-sub g sq)))
3136 (if (math-lessp (calcFunc-abssqr gp) (calcFunc-abssqr gm))
3137 (setq gp gm))
3138 (setq dx (math-div m gp)
3139 x1 (math-sub x dx))
3140 (if (and try-newt
3141 (math-lessp (math-abs-approx dx)
3142 (calcFunc-scf (math-abs-approx x) -3)))
3143 (let ((newt (math-poly-newton-root p x1 7)))
3144 (setq tried-newt t
3145 try-newt nil)
3146 (if (math-zerop (cdr newt))
3147 (setq x (car newt) x1 x)
3148 (if (math-lessp (cdr newt) '(float 1 -6))
3149 (let ((newt2 (math-poly-newton-root
3150 p (car newt) 20)))
3151 (if (math-zerop (cdr newt2))
3152 (setq x (car newt2) x1 x)
3153 (setq x (car newt))))))))
3154 (not (or (eq x x1)
3155 (math-nearly-equal x x1))))
3156 (let ((cdx (math-abs-approx dx)))
3157 (setq x x1
3158 tried-newt nil)
3159 (prog1
3160 (or (<= iters 6)
3161 (math-lessp cdx dxold)
3162 (progn
3163 (if polish
3164 (let ((digs (calcFunc-xpon
3165 (math-div (math-abs-approx x) cdx))))
3166 (calc-record-why
3167 "*Could not attain full precision")
3168 (if (natnump digs)
3169 (let ((calc-internal-prec (max 3 digs)))
3170 (setq x (math-normalize x))))))
3171 nil))
3172 (setq dxold cdx)))
3173 (or polish
3174 (math-lessp (calcFunc-scf (math-abs-approx x)
3175 (- calc-internal-prec))
3176 dxold))))
3177 (or (and (math-floatp x)
3178 (math-poly-integer-root x))
3179 x)))
3181 (defun math-solve-above-dummy (x)
3182 (and (not (Math-primp x))
3183 (if (and (equal (nth 1 x) '(var SOLVEDUM SOLVEDUM))
3184 (= (length x) 2))
3186 (let ((res nil))
3187 (while (and (setq x (cdr x))
3188 (not (setq res (math-solve-above-dummy (car x))))))
3189 res))))
3191 (defun math-solve-find-root-term (x neg) ; sets "t2", "t3"
3192 (if (math-solve-find-root-in-prod x)
3193 (setq math-t3 neg
3194 math-t1 x)
3195 (and (memq (car-safe x) '(+ -))
3196 (or (math-solve-find-root-term (nth 1 x) neg)
3197 (math-solve-find-root-term (nth 2 x)
3198 (if (eq (car x) '-) (not neg) neg))))))
3200 (defun math-solve-find-root-in-prod (x)
3201 (and (consp x)
3202 (math-expr-contains x math-solve-var)
3203 (or (and (eq (car x) 'calcFunc-sqrt)
3204 (setq math-t2 2))
3205 (and (eq (car x) '^)
3206 (or (and (memq (math-quarter-integer (nth 2 x)) '(1 2 3))
3207 (setq math-t2 2))
3208 (and (eq (car-safe (nth 2 x)) 'frac)
3209 (eq (nth 2 (nth 2 x)) 3)
3210 (setq math-t2 3))))
3211 (and (memq (car x) '(* /))
3212 (or (and (not (math-expr-contains (nth 1 x) math-solve-var))
3213 (math-solve-find-root-in-prod (nth 2 x)))
3214 (and (not (math-expr-contains (nth 2 x) math-solve-var))
3215 (math-solve-find-root-in-prod (nth 1 x))))))))
3217 ;; The variable math-solve-vars is local to math-solve-system,
3218 ;; but is used by math-solve-system-rec.
3219 (defvar math-solve-vars)
3221 ;; The variable math-solve-simplifying is local to math-solve-system
3222 ;; and math-solve-system-rec, but is used by math-solve-system-subst.
3223 (defvar math-solve-simplifying)
3225 (defun math-solve-system (exprs math-solve-vars math-solve-full)
3226 (setq exprs (mapcar 'list (if (Math-vectorp exprs)
3227 (cdr exprs)
3228 (list exprs)))
3229 math-solve-vars (if (Math-vectorp math-solve-vars)
3230 (cdr math-solve-vars)
3231 (list math-solve-vars)))
3232 (or (let ((math-solve-simplifying nil))
3233 (math-solve-system-rec exprs math-solve-vars nil))
3234 (let ((math-solve-simplifying t))
3235 (math-solve-system-rec exprs math-solve-vars nil))))
3237 ;;; The following backtracking solver works by choosing a variable
3238 ;;; and equation, and trying to solve the equation for the variable.
3239 ;;; If it succeeds it calls itself recursively with that variable and
3240 ;;; equation removed from their respective lists, and with the solution
3241 ;;; added to solns as well as being substituted into all existing
3242 ;;; equations. The algorithm terminates when any solution path
3243 ;;; manages to remove all the variables from var-list.
3245 ;;; To support calcFunc-roots, entries in eqn-list and solns are
3246 ;;; actually lists of equations.
3248 ;; The variables math-solve-system-res and math-solve-system-vv are
3249 ;; local to math-solve-system-rec, but are used by math-solve-system-subst.
3250 (defvar math-solve-system-vv)
3251 (defvar math-solve-system-res)
3254 (defun math-solve-system-rec (eqn-list var-list solns)
3255 (if var-list
3256 (let ((v var-list)
3257 (math-solve-system-res nil))
3259 ;; Try each variable in turn.
3260 (while
3261 (and
3263 (let* ((math-solve-system-vv (car v))
3264 (e eqn-list)
3265 (elim (eq (car-safe math-solve-system-vv) 'calcFunc-elim)))
3266 (if elim
3267 (setq math-solve-system-vv (nth 1 math-solve-system-vv)))
3269 ;; Try each equation in turn.
3270 (while
3271 (and
3273 (let ((e2 (car e))
3274 (eprev nil)
3275 res2)
3276 (setq math-solve-system-res nil)
3278 ;; Try to solve for math-solve-system-vv the list of equations e2.
3279 (while (and e2
3280 (setq res2 (or (and (eq (car e2) eprev)
3281 res2)
3282 (math-solve-for (car e2) 0
3283 math-solve-system-vv
3284 math-solve-full))))
3285 (setq eprev (car e2)
3286 math-solve-system-res (cons (if (eq math-solve-full 'all)
3287 (cdr res2)
3288 (list res2))
3289 math-solve-system-res)
3290 e2 (cdr e2)))
3291 (if e2
3292 (setq math-solve-system-res nil)
3294 ;; Found a solution. Now try other variables.
3295 (setq math-solve-system-res (nreverse math-solve-system-res)
3296 math-solve-system-res (math-solve-system-rec
3297 (mapcar
3298 'math-solve-system-subst
3299 (delq (car e)
3300 (copy-sequence eqn-list)))
3301 (delq (car v) (copy-sequence var-list))
3302 (let ((math-solve-simplifying nil)
3303 (s (mapcar
3304 (function
3305 (lambda (x)
3306 (cons
3307 (car x)
3308 (math-solve-system-subst
3309 (cdr x)))))
3310 solns)))
3311 (if elim
3313 (cons (cons
3314 math-solve-system-vv
3315 (apply 'append math-solve-system-res))
3316 s)))))
3317 (not math-solve-system-res))))
3318 (setq e (cdr e)))
3319 (not math-solve-system-res)))
3320 (setq v (cdr v)))
3321 math-solve-system-res)
3323 ;; Eliminated all variables, so now put solution into the proper format.
3324 (setq solns (sort solns
3325 (function
3326 (lambda (x y)
3327 (not (memq (car x) (memq (car y) math-solve-vars)))))))
3328 (if (eq math-solve-full 'all)
3329 (math-transpose
3330 (math-normalize
3331 (cons 'vec
3332 (if solns
3333 (mapcar (function (lambda (x) (cons 'vec (cdr x)))) solns)
3334 (mapcar (function (lambda (x) (cons 'vec x))) eqn-list)))))
3335 (math-normalize
3336 (cons 'vec
3337 (if solns
3338 (mapcar (function (lambda (x) (cons 'calcFunc-eq x))) solns)
3339 (mapcar 'car eqn-list)))))))
3341 (defun math-solve-system-subst (x) ; uses "res" and "v"
3342 (let ((accum nil)
3343 (res2 math-solve-system-res))
3344 (while x
3345 (setq accum (nconc accum
3346 (mapcar (function
3347 (lambda (r)
3348 (if math-solve-simplifying
3349 (math-simplify
3350 (math-expr-subst
3351 (car x) math-solve-system-vv r))
3352 (math-expr-subst
3353 (car x) math-solve-system-vv r))))
3354 (car res2)))
3355 x (cdr x)
3356 res2 (cdr res2)))
3357 accum))
3360 ;; calc-command-flags is declared in calc.el
3361 (defvar calc-command-flags)
3363 (defun math-get-from-counter (name)
3364 (let ((ctr (assq name calc-command-flags)))
3365 (if ctr
3366 (setcdr ctr (1+ (cdr ctr)))
3367 (setq ctr (cons name 1)
3368 calc-command-flags (cons ctr calc-command-flags)))
3369 (cdr ctr)))
3371 (defvar var-GenCount)
3373 (defun math-solve-get-sign (val)
3374 (setq val (math-simplify val))
3375 (if (and (eq (car-safe val) '*)
3376 (Math-numberp (nth 1 val)))
3377 (list '* (nth 1 val) (math-solve-get-sign (nth 2 val)))
3378 (and (eq (car-safe val) 'calcFunc-sqrt)
3379 (eq (car-safe (nth 1 val)) '^)
3380 (setq val (math-normalize (list '^
3381 (nth 1 (nth 1 val))
3382 (math-div (nth 2 (nth 1 val)) 2)))))
3383 (if math-solve-full
3384 (if (and (calc-var-value 'var-GenCount)
3385 (Math-natnump var-GenCount)
3386 (not (eq math-solve-full 'all)))
3387 (prog1
3388 (math-mul (list 'calcFunc-as var-GenCount) val)
3389 (setq var-GenCount (math-add var-GenCount 1))
3390 (calc-refresh-evaltos 'var-GenCount))
3391 (let* ((var (concat "s" (int-to-string (math-get-from-counter 'solve-sign))))
3392 (var2 (list 'var (intern var) (intern (concat "var-" var)))))
3393 (if (eq math-solve-full 'all)
3394 (setq math-solve-ranges (cons (list var2 1 -1)
3395 math-solve-ranges)))
3396 (math-mul var2 val)))
3397 (calc-record-why "*Choosing positive solution")
3398 val)))
3400 (defun math-solve-get-int (val &optional range first)
3401 (if math-solve-full
3402 (if (and (calc-var-value 'var-GenCount)
3403 (Math-natnump var-GenCount)
3404 (not (eq math-solve-full 'all)))
3405 (prog1
3406 (math-mul val (list 'calcFunc-an var-GenCount))
3407 (setq var-GenCount (math-add var-GenCount 1))
3408 (calc-refresh-evaltos 'var-GenCount))
3409 (let* ((var (concat "n" (int-to-string
3410 (math-get-from-counter 'solve-int))))
3411 (var2 (list 'var (intern var) (intern (concat "var-" var)))))
3412 (if (and range (eq math-solve-full 'all))
3413 (setq math-solve-ranges (cons (cons var2
3414 (cdr (calcFunc-index
3415 range (or first 0))))
3416 math-solve-ranges)))
3417 (math-mul val var2)))
3418 (calc-record-why "*Choosing 0 for arbitrary integer in solution")
3421 (defun math-solve-sign (sign expr)
3422 (and sign
3423 (let ((s1 (math-possible-signs expr)))
3424 (cond ((memq s1 '(4 6))
3425 sign)
3426 ((memq s1 '(1 3))
3427 (- sign))))))
3429 (defun math-looks-evenp (expr)
3430 (if (Math-integerp expr)
3431 (math-evenp expr)
3432 (if (memq (car expr) '(* /))
3433 (math-looks-evenp (nth 1 expr)))))
3435 (defun math-solve-for (lhs rhs math-solve-var math-solve-full &optional sign)
3436 (if (math-expr-contains rhs math-solve-var)
3437 (math-solve-for (math-sub lhs rhs) 0 math-solve-var math-solve-full)
3438 (and (math-expr-contains lhs math-solve-var)
3439 (math-with-extra-prec 1
3440 (let* ((math-poly-base-variable math-solve-var)
3441 (res (math-try-solve-for lhs rhs sign)))
3442 (if (and (eq math-solve-full 'all)
3443 (math-known-realp math-solve-var))
3444 (let ((old-len (length res))
3445 new-len)
3446 (setq res (delq nil
3447 (mapcar (function
3448 (lambda (x)
3449 (and (not (memq (car-safe x)
3450 '(cplx polar)))
3451 x)))
3452 res))
3453 new-len (length res))
3454 (if (< new-len old-len)
3455 (calc-record-why (if (= new-len 1)
3456 "*All solutions were complex"
3457 (format
3458 "*Omitted %d complex solutions"
3459 (- old-len new-len)))))))
3460 res)))))
3462 (defun math-solve-eqn (expr var full)
3463 (if (memq (car-safe expr) '(calcFunc-neq calcFunc-lt calcFunc-gt
3464 calcFunc-leq calcFunc-geq))
3465 (let ((res (math-solve-for (cons '- (cdr expr))
3466 0 var full
3467 (if (eq (car expr) 'calcFunc-neq) nil 1))))
3468 (and res
3469 (if (eq math-solve-sign 1)
3470 (list (car expr) var res)
3471 (if (eq math-solve-sign -1)
3472 (list (car expr) res var)
3473 (or (eq (car expr) 'calcFunc-neq)
3474 (calc-record-why
3475 "*Can't determine direction of inequality"))
3476 (and (memq (car expr) '(calcFunc-neq calcFunc-lt calcFunc-gt))
3477 (list 'calcFunc-neq var res))))))
3478 (let ((res (math-solve-for expr 0 var full)))
3479 (and res
3480 (list 'calcFunc-eq var res)))))
3482 (defun math-reject-solution (expr var func)
3483 (if (math-expr-contains expr var)
3484 (or (equal (car calc-next-why) '(* "Unable to find a symbolic solution"))
3485 (calc-record-why "*Unable to find a solution")))
3486 (list func expr var))
3488 (defun calcFunc-solve (expr var)
3489 (or (if (or (Math-vectorp expr) (Math-vectorp var))
3490 (math-solve-system expr var nil)
3491 (math-solve-eqn expr var nil))
3492 (math-reject-solution expr var 'calcFunc-solve)))
3494 (defun calcFunc-fsolve (expr var)
3495 (or (if (or (Math-vectorp expr) (Math-vectorp var))
3496 (math-solve-system expr var t)
3497 (math-solve-eqn expr var t))
3498 (math-reject-solution expr var 'calcFunc-fsolve)))
3500 (defun calcFunc-roots (expr var)
3501 (let ((math-solve-ranges nil))
3502 (or (if (or (Math-vectorp expr) (Math-vectorp var))
3503 (math-solve-system expr var 'all)
3504 (math-solve-for expr 0 var 'all))
3505 (math-reject-solution expr var 'calcFunc-roots))))
3507 (defun calcFunc-finv (expr var)
3508 (let ((res (math-solve-for expr math-integ-var var nil)))
3509 (if res
3510 (math-normalize (math-expr-subst res math-integ-var var))
3511 (math-reject-solution expr var 'calcFunc-finv))))
3513 (defun calcFunc-ffinv (expr var)
3514 (let ((res (math-solve-for expr math-integ-var var t)))
3515 (if res
3516 (math-normalize (math-expr-subst res math-integ-var var))
3517 (math-reject-solution expr var 'calcFunc-finv))))
3520 (put 'calcFunc-inv 'math-inverse
3521 (function (lambda (x) (math-div 1 x))))
3522 (put 'calcFunc-inv 'math-inverse-sign -1)
3524 (put 'calcFunc-sqrt 'math-inverse
3525 (function (lambda (x) (math-sqr x))))
3527 (put 'calcFunc-conj 'math-inverse
3528 (function (lambda (x) (list 'calcFunc-conj x))))
3530 (put 'calcFunc-abs 'math-inverse
3531 (function (lambda (x) (math-solve-get-sign x))))
3533 (put 'calcFunc-deg 'math-inverse
3534 (function (lambda (x) (list 'calcFunc-rad x))))
3535 (put 'calcFunc-deg 'math-inverse-sign 1)
3537 (put 'calcFunc-rad 'math-inverse
3538 (function (lambda (x) (list 'calcFunc-deg x))))
3539 (put 'calcFunc-rad 'math-inverse-sign 1)
3541 (put 'calcFunc-ln 'math-inverse
3542 (function (lambda (x) (list 'calcFunc-exp x))))
3543 (put 'calcFunc-ln 'math-inverse-sign 1)
3545 (put 'calcFunc-log10 'math-inverse
3546 (function (lambda (x) (list 'calcFunc-exp10 x))))
3547 (put 'calcFunc-log10 'math-inverse-sign 1)
3549 (put 'calcFunc-lnp1 'math-inverse
3550 (function (lambda (x) (list 'calcFunc-expm1 x))))
3551 (put 'calcFunc-lnp1 'math-inverse-sign 1)
3553 (put 'calcFunc-exp 'math-inverse
3554 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-ln x))
3555 (math-mul 2
3556 (math-mul '(var pi var-pi)
3557 (math-solve-get-int
3558 '(var i var-i))))))))
3559 (put 'calcFunc-exp 'math-inverse-sign 1)
3561 (put 'calcFunc-expm1 'math-inverse
3562 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-lnp1 x))
3563 (math-mul 2
3564 (math-mul '(var pi var-pi)
3565 (math-solve-get-int
3566 '(var i var-i))))))))
3567 (put 'calcFunc-expm1 'math-inverse-sign 1)
3569 (put 'calcFunc-sin 'math-inverse
3570 (function (lambda (x) (let ((n (math-solve-get-int 1)))
3571 (math-add (math-mul (math-normalize
3572 (list 'calcFunc-arcsin x))
3573 (math-pow -1 n))
3574 (math-mul (math-half-circle t)
3575 n))))))
3577 (put 'calcFunc-cos 'math-inverse
3578 (function (lambda (x) (math-add (math-solve-get-sign
3579 (math-normalize
3580 (list 'calcFunc-arccos x)))
3581 (math-solve-get-int
3582 (math-full-circle t))))))
3584 (put 'calcFunc-tan 'math-inverse
3585 (function (lambda (x) (math-add (math-normalize (list 'calcFunc-arctan x))
3586 (math-solve-get-int
3587 (math-half-circle t))))))
3589 (put 'calcFunc-arcsin 'math-inverse
3590 (function (lambda (x) (math-normalize (list 'calcFunc-sin x)))))
3592 (put 'calcFunc-arccos 'math-inverse
3593 (function (lambda (x) (math-normalize (list 'calcFunc-cos x)))))
3595 (put 'calcFunc-arctan 'math-inverse
3596 (function (lambda (x) (math-normalize (list 'calcFunc-tan x)))))
3598 (put 'calcFunc-sinh 'math-inverse
3599 (function (lambda (x) (let ((n (math-solve-get-int 1)))
3600 (math-add (math-mul (math-normalize
3601 (list 'calcFunc-arcsinh x))
3602 (math-pow -1 n))
3603 (math-mul (math-half-circle t)
3604 (math-mul
3605 '(var i var-i)
3606 n)))))))
3607 (put 'calcFunc-sinh 'math-inverse-sign 1)
3609 (put 'calcFunc-cosh 'math-inverse
3610 (function (lambda (x) (math-add (math-solve-get-sign
3611 (math-normalize
3612 (list 'calcFunc-arccosh x)))
3613 (math-mul (math-full-circle t)
3614 (math-solve-get-int
3615 '(var i var-i)))))))
3617 (put 'calcFunc-tanh 'math-inverse
3618 (function (lambda (x) (math-add (math-normalize
3619 (list 'calcFunc-arctanh x))
3620 (math-mul (math-half-circle t)
3621 (math-solve-get-int
3622 '(var i var-i)))))))
3623 (put 'calcFunc-tanh 'math-inverse-sign 1)
3625 (put 'calcFunc-arcsinh 'math-inverse
3626 (function (lambda (x) (math-normalize (list 'calcFunc-sinh x)))))
3627 (put 'calcFunc-arcsinh 'math-inverse-sign 1)
3629 (put 'calcFunc-arccosh 'math-inverse
3630 (function (lambda (x) (math-normalize (list 'calcFunc-cosh x)))))
3632 (put 'calcFunc-arctanh 'math-inverse
3633 (function (lambda (x) (math-normalize (list 'calcFunc-tanh x)))))
3634 (put 'calcFunc-arctanh 'math-inverse-sign 1)
3638 (defun calcFunc-taylor (expr var num)
3639 (let ((x0 0) (v var))
3640 (if (memq (car-safe var) '(+ - calcFunc-eq))
3641 (setq x0 (if (eq (car var) '+) (math-neg (nth 2 var)) (nth 2 var))
3642 v (nth 1 var)))
3643 (or (and (eq (car-safe v) 'var)
3644 (math-expr-contains expr v)
3645 (natnump num)
3646 (let ((accum (math-expr-subst expr v x0))
3647 (var2 (if (eq (car var) 'calcFunc-eq)
3648 (cons '- (cdr var))
3649 var))
3650 (n 0)
3651 (nfac 1)
3652 (fprime expr))
3653 (while (and (<= (setq n (1+ n)) num)
3654 (setq fprime (calcFunc-deriv fprime v nil t)))
3655 (setq fprime (math-simplify fprime)
3656 nfac (math-mul nfac n)
3657 accum (math-add accum
3658 (math-div (math-mul (math-pow var2 n)
3659 (math-expr-subst
3660 fprime v x0))
3661 nfac))))
3662 (and fprime
3663 (math-normalize accum))))
3664 (list 'calcFunc-taylor expr var num))))
3666 (provide 'calcalg2)
3668 ;;; calcalg2.el ends here