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