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