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