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