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