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