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