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