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