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