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