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