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