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