Commit | Line | Data |
---|---|---|
3132f345 CW |
1 | ;;; calc-math.el --- mathematical functions for Calc |
2 | ||
58ba2f8f | 3 | ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004, |
ae940284 | 4 | ;; 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc. |
3132f345 CW |
5 | |
6 | ;; Author: David Gillespie <daveg@synaptics.com> | |
e8fff8ed | 7 | ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com> |
136211a9 EZ |
8 | |
9 | ;; This file is part of GNU Emacs. | |
10 | ||
662c9c64 | 11 | ;; GNU Emacs is free software: you can redistribute it and/or modify |
7c671b23 | 12 | ;; it under the terms of the GNU General Public License as published by |
662c9c64 GM |
13 | ;; the Free Software Foundation, either version 3 of the License, or |
14 | ;; (at your option) any later version. | |
7c671b23 | 15 | |
136211a9 | 16 | ;; GNU Emacs is distributed in the hope that it will be useful, |
7c671b23 GM |
17 | ;; but WITHOUT ANY WARRANTY; without even the implied warranty of |
18 | ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
19 | ;; GNU General Public License for more details. | |
20 | ||
21 | ;; You should have received a copy of the GNU General Public License | |
662c9c64 | 22 | ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>. |
136211a9 | 23 | |
3132f345 | 24 | ;;; Commentary: |
136211a9 | 25 | |
3132f345 | 26 | ;;; Code: |
136211a9 EZ |
27 | |
28 | ;; This file is autoloaded from calc-ext.el. | |
136211a9 | 29 | |
95995a85 | 30 | (require 'calc-ext) |
136211a9 EZ |
31 | (require 'calc-macs) |
32 | ||
178b8baf JB |
33 | |
34 | ;;; Find out how many 9s in 9.9999... will give distinct Emacs floats, | |
35 | ;;; then back off by one. | |
36 | ||
37 | (defvar math-emacs-precision | |
38 | (let* ((n 1) | |
39 | (x 9) | |
40 | (xx (+ x (* 9 (expt 10 (- n)))))) | |
41 | (while (/= x xx) | |
42 | (progn | |
43 | (setq n (1+ n)) | |
44 | (setq x xx) | |
45 | (setq xx (+ x (* 9 (expt 10 (- n))))))) | |
46 | (1- n)) | |
47 | "The number of digits in an Emacs float.") | |
48 | ||
49 | ;;; Find the largest power of 10 which is an Emacs float, | |
50 | ;;; then back off by one so that any float d.dddd...eN | |
51 | ;;; is an Emacs float, for acceptable d.dddd.... | |
52 | ||
53 | (defvar math-largest-emacs-expt | |
86e405cf JB |
54 | (let ((x 1) |
55 | (pow 1e2)) | |
56 | ;; The following loop is for efficiency; it should stop when | |
57 | ;; 10^(2x) is too large. This could be indicated by a range | |
16635351 | 58 | ;; error when computing 10^(2x) or an infinite value for 10^(2x). |
86e405cf JB |
59 | (while (and |
60 | pow | |
16635351 | 61 | (< pow 1.0e+INF)) |
86e405cf JB |
62 | (setq x (* 2 x)) |
63 | (setq pow (condition-case nil | |
64 | (expt 10.0 (* 2 x)) | |
65 | (error nil)))) | |
66 | ;; The following loop should stop when 10^(x+1) is too large. | |
67 | (setq pow (condition-case nil | |
68 | (expt 10.0 (1+ x)) | |
69 | (error nil))) | |
70 | (while (and | |
71 | pow | |
16635351 | 72 | (< pow 1.0e+INF)) |
86e405cf JB |
73 | (setq x (1+ x)) |
74 | (setq pow (condition-case nil | |
75 | (expt 10.0 (1+ x)) | |
76 | (error nil)))) | |
77 | (1- x)) | |
178b8baf JB |
78 | "The largest exponent which Calc will convert to an Emacs float.") |
79 | ||
80 | (defvar math-smallest-emacs-expt | |
81 | (let ((x -1)) | |
82 | (while (condition-case nil | |
f1640784 | 83 | (> (expt 10.0 x) 0.0) |
178b8baf JB |
84 | (error nil)) |
85 | (setq x (* 2 x))) | |
86 | (setq x (/ x 2)) | |
87 | (while (condition-case nil | |
f1640784 | 88 | (> (expt 10.0 x) 0.0) |
178b8baf JB |
89 | (error nil)) |
90 | (setq x (1- x))) | |
91 | (+ x 2)) | |
92 | "The smallest exponent which Calc will convert to an Emacs float.") | |
93 | ||
94 | (defun math-use-emacs-fn (fn x) | |
95 | "Use the native Emacs function FN to evaluate the Calc number X. | |
96 | If this can't be done, return NIL." | |
97 | (and | |
98 | (<= calc-internal-prec math-emacs-precision) | |
99 | (math-realp x) | |
100 | (let* ((fx (math-float x)) | |
101 | (xpon (+ (nth 2 x) (1- (math-numdigs (nth 1 x)))))) | |
102 | (and (<= math-smallest-emacs-expt xpon) | |
103 | (<= xpon math-largest-emacs-expt) | |
104 | (condition-case nil | |
105 | (math-read-number | |
106 | (number-to-string | |
107 | (funcall fn | |
7b3752ab | 108 | (string-to-number |
2e9b968b JB |
109 | (let |
110 | ((calc-number-radix 10) | |
111 | (calc-float-format (list 'float calc-internal-prec)) | |
112 | (calc-group-digits nil) | |
113 | (calc-point-char ".")) | |
7b3752ab | 114 | (math-format-number (math-float x))))))) |
178b8baf JB |
115 | (error nil)))))) |
116 | ||
136211a9 EZ |
117 | (defun calc-sqrt (arg) |
118 | (interactive "P") | |
119 | (calc-slow-wrapper | |
120 | (if (calc-is-inverse) | |
121 | (calc-unary-op "^2" 'calcFunc-sqr arg) | |
491c3062 | 122 | (calc-unary-op "sqrt" 'calcFunc-sqrt arg)))) |
136211a9 EZ |
123 | |
124 | (defun calc-isqrt (arg) | |
125 | (interactive "P") | |
126 | (calc-slow-wrapper | |
127 | (if (calc-is-inverse) | |
128 | (calc-unary-op "^2" 'calcFunc-sqr arg) | |
491c3062 | 129 | (calc-unary-op "isqt" 'calcFunc-isqrt arg)))) |
136211a9 EZ |
130 | |
131 | ||
132 | (defun calc-hypot (arg) | |
133 | (interactive "P") | |
134 | (calc-slow-wrapper | |
491c3062 | 135 | (calc-binary-op "hypt" 'calcFunc-hypot arg))) |
136211a9 EZ |
136 | |
137 | (defun calc-ln (arg) | |
138 | (interactive "P") | |
139 | (calc-invert-func) | |
491c3062 | 140 | (calc-exp arg)) |
136211a9 EZ |
141 | |
142 | (defun calc-log10 (arg) | |
143 | (interactive "P") | |
144 | (calc-hyperbolic-func) | |
491c3062 | 145 | (calc-ln arg)) |
136211a9 EZ |
146 | |
147 | (defun calc-log (arg) | |
148 | (interactive "P") | |
149 | (calc-slow-wrapper | |
150 | (if (calc-is-inverse) | |
151 | (calc-binary-op "alog" 'calcFunc-alog arg) | |
491c3062 | 152 | (calc-binary-op "log" 'calcFunc-log arg)))) |
136211a9 EZ |
153 | |
154 | (defun calc-ilog (arg) | |
155 | (interactive "P") | |
156 | (calc-slow-wrapper | |
157 | (if (calc-is-inverse) | |
158 | (calc-binary-op "alog" 'calcFunc-alog arg) | |
491c3062 | 159 | (calc-binary-op "ilog" 'calcFunc-ilog arg)))) |
136211a9 EZ |
160 | |
161 | (defun calc-lnp1 (arg) | |
162 | (interactive "P") | |
163 | (calc-invert-func) | |
491c3062 | 164 | (calc-expm1 arg)) |
136211a9 EZ |
165 | |
166 | (defun calc-exp (arg) | |
167 | (interactive "P") | |
168 | (calc-slow-wrapper | |
169 | (if (calc-is-hyperbolic) | |
170 | (if (calc-is-inverse) | |
171 | (calc-unary-op "lg10" 'calcFunc-log10 arg) | |
172 | (calc-unary-op "10^" 'calcFunc-exp10 arg)) | |
173 | (if (calc-is-inverse) | |
174 | (calc-unary-op "ln" 'calcFunc-ln arg) | |
491c3062 | 175 | (calc-unary-op "exp" 'calcFunc-exp arg))))) |
136211a9 EZ |
176 | |
177 | (defun calc-expm1 (arg) | |
178 | (interactive "P") | |
179 | (calc-slow-wrapper | |
180 | (if (calc-is-inverse) | |
181 | (calc-unary-op "ln+1" 'calcFunc-lnp1 arg) | |
491c3062 | 182 | (calc-unary-op "ex-1" 'calcFunc-expm1 arg)))) |
136211a9 EZ |
183 | |
184 | (defun calc-pi () | |
185 | (interactive) | |
186 | (calc-slow-wrapper | |
187 | (if (calc-is-inverse) | |
188 | (if (calc-is-hyperbolic) | |
189 | (if calc-symbolic-mode | |
190 | (calc-pop-push-record 0 "phi" '(var phi var-phi)) | |
191 | (calc-pop-push-record 0 "phi" (math-phi))) | |
192 | (if calc-symbolic-mode | |
193 | (calc-pop-push-record 0 "gmma" '(var gamma var-gamma)) | |
194 | (calc-pop-push-record 0 "gmma" (math-gamma-const)))) | |
195 | (if (calc-is-hyperbolic) | |
196 | (if calc-symbolic-mode | |
197 | (calc-pop-push-record 0 "e" '(var e var-e)) | |
198 | (calc-pop-push-record 0 "e" (math-e))) | |
199 | (if calc-symbolic-mode | |
200 | (calc-pop-push-record 0 "pi" '(var pi var-pi)) | |
491c3062 | 201 | (calc-pop-push-record 0 "pi" (math-pi))))))) |
136211a9 EZ |
202 | |
203 | (defun calc-sin (arg) | |
204 | (interactive "P") | |
205 | (calc-slow-wrapper | |
206 | (if (calc-is-hyperbolic) | |
207 | (if (calc-is-inverse) | |
208 | (calc-unary-op "asnh" 'calcFunc-arcsinh arg) | |
209 | (calc-unary-op "sinh" 'calcFunc-sinh arg)) | |
210 | (if (calc-is-inverse) | |
211 | (calc-unary-op "asin" 'calcFunc-arcsin arg) | |
491c3062 | 212 | (calc-unary-op "sin" 'calcFunc-sin arg))))) |
136211a9 EZ |
213 | |
214 | (defun calc-arcsin (arg) | |
215 | (interactive "P") | |
216 | (calc-invert-func) | |
491c3062 | 217 | (calc-sin arg)) |
136211a9 EZ |
218 | |
219 | (defun calc-sinh (arg) | |
220 | (interactive "P") | |
221 | (calc-hyperbolic-func) | |
491c3062 | 222 | (calc-sin arg)) |
136211a9 EZ |
223 | |
224 | (defun calc-arcsinh (arg) | |
225 | (interactive "P") | |
226 | (calc-invert-func) | |
227 | (calc-hyperbolic-func) | |
491c3062 | 228 | (calc-sin arg)) |
136211a9 | 229 | |
f53e6c20 JB |
230 | (defun calc-sec (arg) |
231 | (interactive "P") | |
232 | (calc-slow-wrapper | |
233 | (if (calc-is-hyperbolic) | |
234 | (calc-unary-op "sech" 'calcFunc-sech arg) | |
235 | (calc-unary-op "sec" 'calcFunc-sec arg)))) | |
236 | ||
237 | (defun calc-sech (arg) | |
238 | (interactive "P") | |
239 | (calc-hyperbolic-func) | |
240 | (calc-sec arg)) | |
241 | ||
136211a9 EZ |
242 | (defun calc-cos (arg) |
243 | (interactive "P") | |
244 | (calc-slow-wrapper | |
245 | (if (calc-is-hyperbolic) | |
246 | (if (calc-is-inverse) | |
247 | (calc-unary-op "acsh" 'calcFunc-arccosh arg) | |
248 | (calc-unary-op "cosh" 'calcFunc-cosh arg)) | |
249 | (if (calc-is-inverse) | |
250 | (calc-unary-op "acos" 'calcFunc-arccos arg) | |
491c3062 | 251 | (calc-unary-op "cos" 'calcFunc-cos arg))))) |
136211a9 EZ |
252 | |
253 | (defun calc-arccos (arg) | |
254 | (interactive "P") | |
255 | (calc-invert-func) | |
491c3062 | 256 | (calc-cos arg)) |
136211a9 EZ |
257 | |
258 | (defun calc-cosh (arg) | |
259 | (interactive "P") | |
260 | (calc-hyperbolic-func) | |
491c3062 | 261 | (calc-cos arg)) |
136211a9 EZ |
262 | |
263 | (defun calc-arccosh (arg) | |
264 | (interactive "P") | |
265 | (calc-invert-func) | |
266 | (calc-hyperbolic-func) | |
491c3062 | 267 | (calc-cos arg)) |
136211a9 | 268 | |
f53e6c20 JB |
269 | (defun calc-csc (arg) |
270 | (interactive "P") | |
271 | (calc-slow-wrapper | |
272 | (if (calc-is-hyperbolic) | |
273 | (calc-unary-op "csch" 'calcFunc-csch arg) | |
274 | (calc-unary-op "csc" 'calcFunc-csc arg)))) | |
275 | ||
276 | (defun calc-csch (arg) | |
277 | (interactive "P") | |
278 | (calc-hyperbolic-func) | |
279 | (calc-csc arg)) | |
280 | ||
136211a9 EZ |
281 | (defun calc-sincos () |
282 | (interactive) | |
283 | (calc-slow-wrapper | |
284 | (if (calc-is-inverse) | |
285 | (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1))) | |
491c3062 | 286 | (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1)))))) |
136211a9 EZ |
287 | |
288 | (defun calc-tan (arg) | |
289 | (interactive "P") | |
290 | (calc-slow-wrapper | |
291 | (if (calc-is-hyperbolic) | |
292 | (if (calc-is-inverse) | |
293 | (calc-unary-op "atnh" 'calcFunc-arctanh arg) | |
294 | (calc-unary-op "tanh" 'calcFunc-tanh arg)) | |
295 | (if (calc-is-inverse) | |
296 | (calc-unary-op "atan" 'calcFunc-arctan arg) | |
491c3062 | 297 | (calc-unary-op "tan" 'calcFunc-tan arg))))) |
136211a9 EZ |
298 | |
299 | (defun calc-arctan (arg) | |
300 | (interactive "P") | |
301 | (calc-invert-func) | |
491c3062 | 302 | (calc-tan arg)) |
136211a9 EZ |
303 | |
304 | (defun calc-tanh (arg) | |
305 | (interactive "P") | |
306 | (calc-hyperbolic-func) | |
491c3062 | 307 | (calc-tan arg)) |
136211a9 EZ |
308 | |
309 | (defun calc-arctanh (arg) | |
310 | (interactive "P") | |
311 | (calc-invert-func) | |
312 | (calc-hyperbolic-func) | |
491c3062 | 313 | (calc-tan arg)) |
136211a9 | 314 | |
f53e6c20 JB |
315 | (defun calc-cot (arg) |
316 | (interactive "P") | |
317 | (calc-slow-wrapper | |
318 | (if (calc-is-hyperbolic) | |
319 | (calc-unary-op "coth" 'calcFunc-coth arg) | |
320 | (calc-unary-op "cot" 'calcFunc-cot arg)))) | |
321 | ||
eb1ef455 | 322 | (defun calc-coth (arg) |
f53e6c20 JB |
323 | (interactive "P") |
324 | (calc-hyperbolic-func) | |
eb1ef455 | 325 | (calc-cot arg)) |
f53e6c20 | 326 | |
136211a9 EZ |
327 | (defun calc-arctan2 () |
328 | (interactive) | |
329 | (calc-slow-wrapper | |
491c3062 | 330 | (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2))))) |
136211a9 EZ |
331 | |
332 | (defun calc-conj (arg) | |
333 | (interactive "P") | |
334 | (calc-wrapper | |
491c3062 | 335 | (calc-unary-op "conj" 'calcFunc-conj arg))) |
136211a9 EZ |
336 | |
337 | (defun calc-imaginary () | |
338 | (interactive) | |
339 | (calc-slow-wrapper | |
491c3062 | 340 | (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1))))) |
136211a9 | 341 | |
136211a9 EZ |
342 | (defun calc-to-degrees (arg) |
343 | (interactive "P") | |
344 | (calc-wrapper | |
491c3062 | 345 | (calc-unary-op ">deg" 'calcFunc-deg arg))) |
136211a9 EZ |
346 | |
347 | (defun calc-to-radians (arg) | |
348 | (interactive "P") | |
349 | (calc-wrapper | |
491c3062 | 350 | (calc-unary-op ">rad" 'calcFunc-rad arg))) |
136211a9 EZ |
351 | |
352 | ||
353 | (defun calc-degrees-mode (arg) | |
354 | (interactive "p") | |
355 | (cond ((= arg 1) | |
356 | (calc-wrapper | |
357 | (calc-change-mode 'calc-angle-mode 'deg) | |
3132f345 | 358 | (message "Angles measured in degrees"))) |
136211a9 EZ |
359 | ((= arg 2) (calc-radians-mode)) |
360 | ((= arg 3) (calc-hms-mode)) | |
491c3062 | 361 | (t (error "Prefix argument out of range")))) |
136211a9 EZ |
362 | |
363 | (defun calc-radians-mode () | |
364 | (interactive) | |
365 | (calc-wrapper | |
366 | (calc-change-mode 'calc-angle-mode 'rad) | |
3132f345 | 367 | (message "Angles measured in radians"))) |
136211a9 EZ |
368 | |
369 | ||
370 | ;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public] | |
371 | ;;; This method takes advantage of the fact that Newton's method starting | |
372 | ;;; with an overestimate always works, even using truncating integer division! | |
373 | (defun math-isqrt (a) | |
374 | (cond ((Math-zerop a) a) | |
375 | ((not (math-natnump a)) | |
376 | (math-reject-arg a 'natnump)) | |
377 | ((integerp a) | |
378 | (math-isqrt-small a)) | |
379 | (t | |
491c3062 | 380 | (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a)))))))) |
136211a9 EZ |
381 | |
382 | (defun calcFunc-isqrt (a) | |
383 | (if (math-realp a) | |
384 | (math-isqrt (math-floor a)) | |
491c3062 | 385 | (math-floor (math-sqrt a)))) |
136211a9 EZ |
386 | |
387 | ||
f0529b5b | 388 | ;;; This returns (flag . result) where the flag is t if A is a perfect square. |
136211a9 EZ |
389 | (defun math-isqrt-bignum (a) ; [P.l L] |
390 | (let ((len (length a))) | |
391 | (if (= (% len 2) 0) | |
392 | (let* ((top (nthcdr (- len 2) a))) | |
393 | (math-isqrt-bignum-iter | |
394 | a | |
98888d77 | 395 | (math-scale-bignum-digit-size |
136211a9 EZ |
396 | (math-bignum-big |
397 | (1+ (math-isqrt-small | |
98888d77 | 398 | (+ (* (nth 1 top) math-bignum-digit-size) (car top))))) |
136211a9 EZ |
399 | (1- (/ len 2))))) |
400 | (let* ((top (nth (1- len) a))) | |
401 | (math-isqrt-bignum-iter | |
402 | a | |
98888d77 | 403 | (math-scale-bignum-digit-size |
136211a9 | 404 | (list (1+ (math-isqrt-small top))) |
491c3062 | 405 | (/ len 2))))))) |
136211a9 EZ |
406 | |
407 | (defun math-isqrt-bignum-iter (a guess) ; [l L l] | |
408 | (math-working "isqrt" (cons 'bigpos guess)) | |
409 | (let* ((q (math-div-bignum a guess)) | |
410 | (s (math-add-bignum (car q) guess)) | |
411 | (g2 (math-div2-bignum s)) | |
412 | (comp (math-compare-bignum g2 guess))) | |
413 | (if (< comp 0) | |
414 | (math-isqrt-bignum-iter a g2) | |
415 | (cons (and (= comp 0) | |
416 | (math-zerop-bignum (cdr q)) | |
417 | (= (% (car s) 2) 0)) | |
491c3062 | 418 | guess)))) |
136211a9 EZ |
419 | |
420 | (defun math-zerop-bignum (a) | |
421 | (and (eq (car a) 0) | |
422 | (progn | |
423 | (while (eq (car (setq a (cdr a))) 0)) | |
491c3062 | 424 | (null a)))) |
136211a9 | 425 | |
98888d77 | 426 | (defun math-scale-bignum-digit-size (a n) ; [L L S] |
136211a9 EZ |
427 | (while (> n 0) |
428 | (setq a (cons 0 a) | |
429 | n (1- n))) | |
491c3062 | 430 | a) |
136211a9 EZ |
431 | |
432 | (defun math-isqrt-small (a) ; A > 0. [S S] | |
98888d77 JB |
433 | (let ((g (cond ((>= a 1000000) 10000) |
434 | ((>= a 10000) 1000) | |
136211a9 EZ |
435 | ((>= a 100) 100) |
436 | (t 10))) | |
437 | g2) | |
438 | (while (< (setq g2 (/ (+ g (/ a g)) 2)) g) | |
439 | (setq g g2)) | |
491c3062 | 440 | g)) |
136211a9 EZ |
441 | |
442 | ||
443 | ||
444 | ||
445 | ;;; Compute the square root of a number. | |
446 | ;;; [T N] if possible, else [F N] if possible, else [C N]. [Public] | |
447 | (defun math-sqrt (a) | |
448 | (or | |
449 | (and (Math-zerop a) a) | |
450 | (and (math-known-nonposp a) | |
451 | (math-imaginary (math-sqrt (math-neg a)))) | |
452 | (and (integerp a) | |
453 | (let ((sqrt (math-isqrt-small a))) | |
454 | (if (= (* sqrt sqrt) a) | |
455 | sqrt | |
456 | (if calc-symbolic-mode | |
457 | (list 'calcFunc-sqrt a) | |
458 | (math-sqrt-float (math-float a) (math-float sqrt)))))) | |
459 | (and (eq (car-safe a) 'bigpos) | |
460 | (let* ((res (math-isqrt-bignum (cdr a))) | |
461 | (sqrt (math-normalize (cons 'bigpos (cdr res))))) | |
462 | (if (car res) | |
463 | sqrt | |
464 | (if calc-symbolic-mode | |
465 | (list 'calcFunc-sqrt a) | |
466 | (math-sqrt-float (math-float a) (math-float sqrt)))))) | |
467 | (and (eq (car-safe a) 'frac) | |
468 | (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a))))) | |
469 | (num-sqrt (math-normalize (cons 'bigpos (cdr num-res)))) | |
470 | (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a))))) | |
471 | (den-sqrt (math-normalize (cons 'bigpos (cdr den-res))))) | |
472 | (if (and (car num-res) (car den-res)) | |
473 | (list 'frac num-sqrt den-sqrt) | |
474 | (if calc-symbolic-mode | |
475 | (if (or (car num-res) (car den-res)) | |
476 | (math-div (if (car num-res) | |
477 | num-sqrt (list 'calcFunc-sqrt (nth 1 a))) | |
478 | (if (car den-res) | |
479 | den-sqrt (list 'calcFunc-sqrt (nth 2 a)))) | |
480 | (list 'calcFunc-sqrt a)) | |
481 | (math-sqrt-float (math-float a) | |
482 | (math-div (math-float num-sqrt) den-sqrt)))))) | |
483 | (and (eq (car-safe a) 'float) | |
484 | (if calc-symbolic-mode | |
485 | (if (= (% (nth 2 a) 2) 0) | |
486 | (let ((res (math-isqrt-bignum | |
487 | (cdr (Math-bignum-test (nth 1 a)))))) | |
488 | (if (car res) | |
489 | (math-make-float (math-normalize | |
490 | (cons 'bigpos (cdr res))) | |
491 | (/ (nth 2 a) 2)) | |
492 | (signal 'inexact-result nil))) | |
493 | (signal 'inexact-result nil)) | |
494 | (math-sqrt-float a))) | |
495 | (and (eq (car-safe a) 'cplx) | |
496 | (math-with-extra-prec 2 | |
497 | (let* ((d (math-abs a)) | |
498 | (imag (math-sqrt (math-mul (math-sub d (nth 1 a)) | |
499 | '(float 5 -1))))) | |
500 | (list 'cplx | |
501 | (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1))) | |
502 | (if (math-negp (nth 2 a)) (math-neg imag) imag))))) | |
503 | (and (eq (car-safe a) 'polar) | |
504 | (list 'polar | |
505 | (math-sqrt (nth 1 a)) | |
506 | (math-mul (nth 2 a) '(float 5 -1)))) | |
507 | (and (eq (car-safe a) 'sdev) | |
508 | (let ((sqrt (math-sqrt (nth 1 a)))) | |
509 | (math-make-sdev sqrt | |
510 | (math-div (nth 2 a) (math-mul sqrt 2))))) | |
511 | (and (eq (car-safe a) 'intv) | |
512 | (not (math-negp (nth 2 a))) | |
513 | (math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a)))) | |
514 | (and (eq (car-safe a) '*) | |
515 | (or (math-known-nonnegp (nth 1 a)) | |
516 | (math-known-nonnegp (nth 2 a))) | |
517 | (math-mul (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a)))) | |
518 | (and (eq (car-safe a) '/) | |
519 | (or (and (math-known-nonnegp (nth 2 a)) | |
520 | (math-div (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a)))) | |
521 | (and (math-known-nonnegp (nth 1 a)) | |
522 | (not (math-equal-int (nth 1 a) 1)) | |
523 | (math-mul (math-sqrt (nth 1 a)) | |
524 | (math-sqrt (math-div 1 (nth 2 a))))))) | |
525 | (and (eq (car-safe a) '^) | |
526 | (math-known-evenp (nth 2 a)) | |
527 | (math-known-realp (nth 1 a)) | |
528 | (math-abs (math-pow (nth 1 a) (math-div (nth 2 a) 2)))) | |
529 | (let ((inf (math-infinitep a))) | |
530 | (and inf | |
531 | (math-mul (math-sqrt (math-infinite-dir a inf)) inf))) | |
532 | (progn | |
533 | (calc-record-why 'numberp a) | |
491c3062 | 534 | (list 'calcFunc-sqrt a)))) |
3132f345 | 535 | (defalias 'calcFunc-sqrt 'math-sqrt) |
136211a9 EZ |
536 | |
537 | (defun math-infinite-dir (a &optional inf) | |
538 | (or inf (setq inf (math-infinitep a))) | |
491c3062 | 539 | (math-normalize (math-expr-subst a inf 1))) |
136211a9 EZ |
540 | |
541 | (defun math-sqrt-float (a &optional guess) ; [F F F] | |
542 | (if calc-symbolic-mode | |
543 | (signal 'inexact-result nil) | |
491c3062 | 544 | (math-with-extra-prec 1 (math-sqrt-raw a guess)))) |
136211a9 EZ |
545 | |
546 | (defun math-sqrt-raw (a &optional guess) ; [F F F] | |
547 | (if (not (Math-posp a)) | |
548 | (math-sqrt a) | |
17cc361e JB |
549 | (cond |
550 | ((math-use-emacs-fn 'sqrt a)) | |
551 | (t | |
552 | (if (null guess) | |
553 | (let ((ldiff (- (math-numdigs (nth 1 a)) 6))) | |
554 | (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff))) | |
555 | (setq guess (math-make-float (math-isqrt-small | |
556 | (math-scale-int (nth 1 a) (- ldiff))) | |
557 | (/ (+ (nth 2 a) ldiff) 2))))) | |
558 | (math-sqrt-float-iter a guess))))) | |
136211a9 EZ |
559 | |
560 | (defun math-sqrt-float-iter (a guess) ; [F F F] | |
561 | (math-working "sqrt" guess) | |
562 | (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess)) | |
563 | '(float 5 -1)))) | |
564 | (if (math-nearly-equal-float g2 guess) | |
565 | g2 | |
491c3062 | 566 | (math-sqrt-float-iter a g2)))) |
136211a9 EZ |
567 | |
568 | ;;; True if A and B differ only in the last digit of precision. [P F F] | |
569 | (defun math-nearly-equal-float (a b) | |
570 | (let ((ediff (- (nth 2 a) (nth 2 b)))) | |
571 | (cond ((= ediff 0) ;; Expanded out for speed | |
572 | (setq ediff (math-add (Math-integer-neg (nth 1 a)) (nth 1 b))) | |
573 | (or (eq ediff 0) | |
574 | (and (not (consp ediff)) | |
575 | (< ediff 10) | |
576 | (> ediff -10) | |
577 | (= (math-numdigs (nth 1 a)) calc-internal-prec)))) | |
578 | ((= ediff 1) | |
579 | (setq ediff (math-add (Math-integer-neg (nth 1 b)) | |
580 | (math-scale-int (nth 1 a) 1))) | |
581 | (and (not (consp ediff)) | |
582 | (< ediff 10) | |
583 | (> ediff -10) | |
584 | (= (math-numdigs (nth 1 b)) calc-internal-prec))) | |
585 | ((= ediff -1) | |
586 | (setq ediff (math-add (Math-integer-neg (nth 1 a)) | |
587 | (math-scale-int (nth 1 b) 1))) | |
588 | (and (not (consp ediff)) | |
589 | (< ediff 10) | |
590 | (> ediff -10) | |
491c3062 | 591 | (= (math-numdigs (nth 1 a)) calc-internal-prec)))))) |
136211a9 EZ |
592 | |
593 | (defun math-nearly-equal (a b) ; [P N N] [Public] | |
594 | (setq a (math-float a)) | |
595 | (setq b (math-float b)) | |
596 | (if (eq (car a) 'polar) (setq a (math-complex a))) | |
597 | (if (eq (car b) 'polar) (setq b (math-complex b))) | |
598 | (if (eq (car a) 'cplx) | |
599 | (if (eq (car b) 'cplx) | |
600 | (and (or (math-nearly-equal-float (nth 1 a) (nth 1 b)) | |
601 | (and (math-nearly-zerop-float (nth 1 a) (nth 2 a)) | |
602 | (math-nearly-zerop-float (nth 1 b) (nth 2 b)))) | |
603 | (or (math-nearly-equal-float (nth 2 a) (nth 2 b)) | |
604 | (and (math-nearly-zerop-float (nth 2 a) (nth 1 a)) | |
605 | (math-nearly-zerop-float (nth 2 b) (nth 1 b))))) | |
606 | (and (math-nearly-equal-float (nth 1 a) b) | |
607 | (math-nearly-zerop-float (nth 2 a) b))) | |
608 | (if (eq (car b) 'cplx) | |
609 | (and (math-nearly-equal-float a (nth 1 b)) | |
610 | (math-nearly-zerop-float a (nth 2 b))) | |
491c3062 | 611 | (math-nearly-equal-float a b)))) |
136211a9 EZ |
612 | |
613 | ;;; True if A is nearly zero compared to B. [P F F] | |
614 | (defun math-nearly-zerop-float (a b) | |
615 | (or (eq (nth 1 a) 0) | |
616 | (<= (+ (math-numdigs (nth 1 a)) (nth 2 a)) | |
491c3062 | 617 | (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec))))) |
136211a9 EZ |
618 | |
619 | (defun math-nearly-zerop (a b) ; [P N R] [Public] | |
620 | (setq a (math-float a)) | |
621 | (setq b (math-float b)) | |
622 | (if (eq (car a) 'cplx) | |
623 | (and (math-nearly-zerop-float (nth 1 a) b) | |
624 | (math-nearly-zerop-float (nth 2 a) b)) | |
625 | (if (eq (car a) 'polar) | |
626 | (math-nearly-zerop-float (nth 1 a) b) | |
491c3062 | 627 | (math-nearly-zerop-float a b)))) |
136211a9 EZ |
628 | |
629 | ;;; This implementation could be improved, accuracy-wise. | |
630 | (defun math-hypot (a b) | |
631 | (cond ((Math-zerop a) (math-abs b)) | |
632 | ((Math-zerop b) (math-abs a)) | |
633 | ((not (Math-scalarp a)) | |
634 | (if (math-infinitep a) | |
635 | (if (math-infinitep b) | |
636 | (if (equal a b) | |
637 | a | |
638 | '(var nan var-nan)) | |
639 | a) | |
640 | (calc-record-why 'scalarp a) | |
641 | (list 'calcFunc-hypot a b))) | |
642 | ((not (Math-scalarp b)) | |
643 | (if (math-infinitep b) | |
644 | b | |
645 | (calc-record-why 'scalarp b) | |
646 | (list 'calcFunc-hypot a b))) | |
647 | ((and (Math-numberp a) (Math-numberp b)) | |
648 | (math-with-extra-prec 1 | |
649 | (math-sqrt (math-add (calcFunc-abssqr a) (calcFunc-abssqr b))))) | |
650 | ((eq (car-safe a) 'hms) | |
651 | (if (eq (car-safe b) 'hms) ; this helps sdev's of hms forms | |
652 | (math-to-hms (math-hypot (math-from-hms a 'deg) | |
653 | (math-from-hms b 'deg))) | |
654 | (math-to-hms (math-hypot (math-from-hms a 'deg) b)))) | |
655 | ((eq (car-safe b) 'hms) | |
656 | (math-to-hms (math-hypot a (math-from-hms b 'deg)))) | |
491c3062 | 657 | (t nil))) |
3132f345 | 658 | (defalias 'calcFunc-hypot 'math-hypot) |
136211a9 EZ |
659 | |
660 | (defun calcFunc-sqr (x) | |
491c3062 | 661 | (math-pow x 2)) |
136211a9 EZ |
662 | |
663 | ||
664 | ||
665 | (defun math-nth-root (a n) | |
666 | (cond ((= n 2) (math-sqrt a)) | |
667 | ((Math-zerop a) a) | |
668 | ((Math-negp a) nil) | |
669 | ((Math-integerp a) | |
670 | (let ((root (math-nth-root-integer a n))) | |
671 | (if (car root) | |
672 | (cdr root) | |
673 | (and (not calc-symbolic-mode) | |
674 | (math-nth-root-float (math-float a) n | |
675 | (math-float (cdr root))))))) | |
676 | ((eq (car-safe a) 'frac) | |
677 | (let* ((num-root (math-nth-root-integer (nth 1 a) n)) | |
678 | (den-root (math-nth-root-integer (nth 2 a) n))) | |
679 | (if (and (car num-root) (car den-root)) | |
680 | (list 'frac (cdr num-root) (cdr den-root)) | |
681 | (and (not calc-symbolic-mode) | |
682 | (math-nth-root-float | |
683 | (math-float a) n | |
684 | (math-div-float (math-float (cdr num-root)) | |
685 | (math-float (cdr den-root)))))))) | |
686 | ((eq (car-safe a) 'float) | |
687 | (and (not calc-symbolic-mode) | |
688 | (math-nth-root-float a n))) | |
689 | ((eq (car-safe a) 'polar) | |
690 | (let ((root (math-nth-root (nth 1 a) n))) | |
691 | (and root (list 'polar root (math-div (nth 2 a) n))))) | |
491c3062 | 692 | (t nil))) |
136211a9 | 693 | |
86498823 JB |
694 | ;; The variables math-nrf-n, math-nrf-nf and math-nrf-nfm1 are local |
695 | ;; to math-nth-root-float, but are used by math-nth-root-float-iter, | |
696 | ;; which is called by math-nth-root-float. | |
697 | (defvar math-nrf-n) | |
698 | (defvar math-nrf-nf) | |
699 | (defvar math-nrf-nfm1) | |
700 | ||
701 | (defun math-nth-root-float (a math-nrf-n &optional guess) | |
136211a9 EZ |
702 | (math-inexact-result) |
703 | (math-with-extra-prec 1 | |
86498823 JB |
704 | (let ((math-nrf-nf (math-float math-nrf-n)) |
705 | (math-nrf-nfm1 (math-float (1- math-nrf-n)))) | |
136211a9 EZ |
706 | (math-nth-root-float-iter a (or guess |
707 | (math-make-float | |
708 | 1 (/ (+ (math-numdigs (nth 1 a)) | |
709 | (nth 2 a) | |
86498823 JB |
710 | (/ math-nrf-n 2)) |
711 | math-nrf-n))))))) | |
136211a9 | 712 | |
86498823 | 713 | (defun math-nth-root-float-iter (a guess) |
136211a9 | 714 | (math-working "root" guess) |
86498823 | 715 | (let ((g2 (math-div-float (math-add-float (math-mul math-nrf-nfm1 guess) |
136211a9 | 716 | (math-div-float |
86498823 JB |
717 | a (math-ipow guess (1- math-nrf-n)))) |
718 | math-nrf-nf))) | |
136211a9 EZ |
719 | (if (math-nearly-equal-float g2 guess) |
720 | g2 | |
491c3062 | 721 | (math-nth-root-float-iter a g2)))) |
136211a9 | 722 | |
86498823 JB |
723 | ;; The variable math-nri-n is local to math-nth-root-integer, but |
724 | ;; is used by math-nth-root-int-iter, which is called by | |
725 | ;; math-nth-root-int. | |
726 | (defvar math-nri-n) | |
727 | ||
728 | (defun math-nth-root-integer (a math-nri-n &optional guess) ; [I I S] | |
136211a9 EZ |
729 | (math-nth-root-int-iter a (or guess |
730 | (math-scale-int 1 (/ (+ (math-numdigs a) | |
86498823 JB |
731 | (1- math-nri-n)) |
732 | math-nri-n))))) | |
136211a9 | 733 | |
86498823 | 734 | (defun math-nth-root-int-iter (a guess) |
136211a9 | 735 | (math-working "root" guess) |
86498823 JB |
736 | (let* ((q (math-idivmod a (math-ipow guess (1- math-nri-n)))) |
737 | (s (math-add (car q) (math-mul (1- math-nri-n) guess))) | |
738 | (g2 (math-idivmod s math-nri-n))) | |
136211a9 EZ |
739 | (if (Math-natnum-lessp (car g2) guess) |
740 | (math-nth-root-int-iter a (car g2)) | |
741 | (cons (and (equal (car g2) guess) | |
742 | (eq (cdr q) 0) | |
743 | (eq (cdr g2) 0)) | |
491c3062 | 744 | guess)))) |
136211a9 EZ |
745 | |
746 | (defun calcFunc-nroot (x n) | |
747 | (calcFunc-pow x (if (integerp n) | |
748 | (math-make-frac 1 n) | |
491c3062 | 749 | (math-div 1 n)))) |
136211a9 EZ |
750 | |
751 | ||
752 | ||
753 | ||
754 | ;;;; Transcendental functions. | |
755 | ||
756 | ;;; All of these functions are defined on the complex plane. | |
757 | ;;; (Branch cuts, etc. follow Steele's Common Lisp book.) | |
758 | ||
759 | ;;; Most functions increase calc-internal-prec by 2 digits, then round | |
760 | ;;; down afterward. "-raw" functions use the current precision, require | |
761 | ;;; their arguments to be in float (or complex float) format, and always | |
762 | ;;; work in radians (where applicable). | |
763 | ||
764 | (defun math-to-radians (a) ; [N N] | |
765 | (cond ((eq (car-safe a) 'hms) | |
766 | (math-from-hms a 'rad)) | |
767 | ((memq calc-angle-mode '(deg hms)) | |
768 | (math-mul a (math-pi-over-180))) | |
491c3062 | 769 | (t a))) |
136211a9 EZ |
770 | |
771 | (defun math-from-radians (a) ; [N N] | |
772 | (cond ((eq calc-angle-mode 'deg) | |
773 | (if (math-constp a) | |
774 | (math-div a (math-pi-over-180)) | |
775 | (list 'calcFunc-deg a))) | |
776 | ((eq calc-angle-mode 'hms) | |
777 | (math-to-hms a 'rad)) | |
491c3062 | 778 | (t a))) |
136211a9 EZ |
779 | |
780 | (defun math-to-radians-2 (a) ; [N N] | |
781 | (cond ((eq (car-safe a) 'hms) | |
782 | (math-from-hms a 'rad)) | |
783 | ((memq calc-angle-mode '(deg hms)) | |
784 | (if calc-symbolic-mode | |
785 | (math-div (math-mul a '(var pi var-pi)) 180) | |
786 | (math-mul a (math-pi-over-180)))) | |
491c3062 | 787 | (t a))) |
136211a9 EZ |
788 | |
789 | (defun math-from-radians-2 (a) ; [N N] | |
790 | (cond ((memq calc-angle-mode '(deg hms)) | |
791 | (if calc-symbolic-mode | |
792 | (math-div (math-mul 180 a) '(var pi var-pi)) | |
793 | (math-div a (math-pi-over-180)))) | |
491c3062 | 794 | (t a))) |
136211a9 EZ |
795 | |
796 | ||
797 | ||
798 | ;;; Sine, cosine, and tangent. | |
799 | ||
800 | (defun calcFunc-sin (x) ; [N N] [Public] | |
801 | (cond ((and (integerp x) | |
802 | (if (eq calc-angle-mode 'deg) | |
803 | (= (% x 90) 0) | |
804 | (= x 0))) | |
805 | (aref [0 1 0 -1] (math-mod (/ x 90) 4))) | |
806 | ((Math-scalarp x) | |
807 | (math-with-extra-prec 2 | |
808 | (math-sin-raw (math-to-radians (math-float x))))) | |
809 | ((eq (car x) 'sdev) | |
810 | (if (math-constp x) | |
811 | (math-with-extra-prec 2 | |
812 | (let* ((xx (math-to-radians (math-float (nth 1 x)))) | |
813 | (xs (math-to-radians (math-float (nth 2 x)))) | |
814 | (sc (math-sin-cos-raw xx))) | |
815 | (math-make-sdev (car sc) (math-mul xs (cdr sc))))) | |
816 | (math-make-sdev (calcFunc-sin (nth 1 x)) | |
817 | (math-mul (nth 2 x) (calcFunc-cos (nth 1 x)))))) | |
818 | ((and (eq (car x) 'intv) (math-intv-constp x)) | |
819 | (calcFunc-cos (math-sub x (math-quarter-circle nil)))) | |
820 | ((equal x '(var nan var-nan)) | |
821 | x) | |
822 | (t (calc-record-why 'scalarp x) | |
491c3062 | 823 | (list 'calcFunc-sin x)))) |
136211a9 EZ |
824 | |
825 | (defun calcFunc-cos (x) ; [N N] [Public] | |
826 | (cond ((and (integerp x) | |
827 | (if (eq calc-angle-mode 'deg) | |
828 | (= (% x 90) 0) | |
829 | (= x 0))) | |
830 | (aref [1 0 -1 0] (math-mod (/ x 90) 4))) | |
831 | ((Math-scalarp x) | |
832 | (math-with-extra-prec 2 | |
833 | (math-cos-raw (math-to-radians (math-float x))))) | |
834 | ((eq (car x) 'sdev) | |
835 | (if (math-constp x) | |
836 | (math-with-extra-prec 2 | |
837 | (let* ((xx (math-to-radians (math-float (nth 1 x)))) | |
838 | (xs (math-to-radians (math-float (nth 2 x)))) | |
839 | (sc (math-sin-cos-raw xx))) | |
840 | (math-make-sdev (cdr sc) (math-mul xs (car sc))))) | |
841 | (math-make-sdev (calcFunc-cos (nth 1 x)) | |
842 | (math-mul (nth 2 x) (calcFunc-sin (nth 1 x)))))) | |
843 | ((and (eq (car x) 'intv) (math-intv-constp x)) | |
844 | (math-with-extra-prec 2 | |
845 | (let* ((xx (math-to-radians (math-float x))) | |
846 | (na (math-floor (math-div (nth 2 xx) (math-pi)))) | |
847 | (nb (math-floor (math-div (nth 3 xx) (math-pi)))) | |
848 | (span (math-sub nb na))) | |
849 | (if (memq span '(0 1)) | |
850 | (let ((int (math-sort-intv (nth 1 x) | |
851 | (math-cos-raw (nth 2 xx)) | |
852 | (math-cos-raw (nth 3 xx))))) | |
853 | (if (eq span 1) | |
854 | (if (math-evenp na) | |
855 | (math-make-intv (logior (nth 1 x) 2) | |
856 | -1 | |
857 | (nth 3 int)) | |
858 | (math-make-intv (logior (nth 1 x) 1) | |
859 | (nth 2 int) | |
860 | 1)) | |
861 | int)) | |
862 | (list 'intv 3 -1 1))))) | |
863 | ((equal x '(var nan var-nan)) | |
864 | x) | |
865 | (t (calc-record-why 'scalarp x) | |
491c3062 | 866 | (list 'calcFunc-cos x)))) |
136211a9 EZ |
867 | |
868 | (defun calcFunc-sincos (x) ; [V N] [Public] | |
869 | (if (Math-scalarp x) | |
870 | (math-with-extra-prec 2 | |
871 | (let ((sc (math-sin-cos-raw (math-to-radians (math-float x))))) | |
872 | (list 'vec (cdr sc) (car sc)))) ; the vector [cos, sin] | |
491c3062 | 873 | (list 'vec (calcFunc-sin x) (calcFunc-cos x)))) |
136211a9 EZ |
874 | |
875 | (defun calcFunc-tan (x) ; [N N] [Public] | |
876 | (cond ((and (integerp x) | |
877 | (if (eq calc-angle-mode 'deg) | |
878 | (= (% x 180) 0) | |
879 | (= x 0))) | |
880 | 0) | |
881 | ((Math-scalarp x) | |
882 | (math-with-extra-prec 2 | |
883 | (math-tan-raw (math-to-radians (math-float x))))) | |
884 | ((eq (car x) 'sdev) | |
885 | (if (math-constp x) | |
886 | (math-with-extra-prec 2 | |
887 | (let* ((xx (math-to-radians (math-float (nth 1 x)))) | |
888 | (xs (math-to-radians (math-float (nth 2 x)))) | |
889 | (sc (math-sin-cos-raw xx))) | |
890 | (if (and (math-zerop (cdr sc)) (not calc-infinite-mode)) | |
891 | (progn | |
892 | (calc-record-why "*Division by zero") | |
893 | (list 'calcFunc-tan x)) | |
894 | (math-make-sdev (math-div-float (car sc) (cdr sc)) | |
895 | (math-div-float xs (math-sqr (cdr sc))))))) | |
896 | (math-make-sdev (calcFunc-tan (nth 1 x)) | |
897 | (math-div (nth 2 x) | |
898 | (math-sqr (calcFunc-cos (nth 1 x))))))) | |
899 | ((and (eq (car x) 'intv) (math-intv-constp x)) | |
900 | (or (math-with-extra-prec 2 | |
901 | (let* ((xx (math-to-radians (math-float x))) | |
902 | (na (math-floor (math-div (math-sub (nth 2 xx) | |
903 | (math-pi-over-2)) | |
904 | (math-pi)))) | |
905 | (nb (math-floor (math-div (math-sub (nth 3 xx) | |
906 | (math-pi-over-2)) | |
907 | (math-pi))))) | |
908 | (and (equal na nb) | |
909 | (math-sort-intv (nth 1 x) | |
910 | (math-tan-raw (nth 2 xx)) | |
911 | (math-tan-raw (nth 3 xx)))))) | |
912 | '(intv 3 (neg (var inf var-inf)) (var inf var-inf)))) | |
913 | ((equal x '(var nan var-nan)) | |
914 | x) | |
915 | (t (calc-record-why 'scalarp x) | |
491c3062 | 916 | (list 'calcFunc-tan x)))) |
136211a9 | 917 | |
f53e6c20 JB |
918 | (defun calcFunc-sec (x) |
919 | (cond ((and (integerp x) | |
920 | (eq calc-angle-mode 'deg) | |
921 | (= (% x 180) 0)) | |
922 | (if (= (% x 360) 0) | |
923 | 1 | |
924 | -1)) | |
925 | ((and (integerp x) | |
926 | (eq calc-angle-mode 'rad) | |
927 | (= x 0)) | |
928 | 1) | |
929 | ((Math-scalarp x) | |
930 | (math-with-extra-prec 2 | |
931 | (math-sec-raw (math-to-radians (math-float x))))) | |
932 | ((eq (car x) 'sdev) | |
933 | (if (math-constp x) | |
934 | (math-with-extra-prec 2 | |
935 | (let* ((xx (math-to-radians (math-float (nth 1 x)))) | |
936 | (xs (math-to-radians (math-float (nth 2 x)))) | |
937 | (sc (math-sin-cos-raw xx))) | |
938 | (if (and (math-zerop (cdr sc)) | |
939 | (not calc-infinite-mode)) | |
940 | (progn | |
941 | (calc-record-why "*Division by zero") | |
942 | (list 'calcFunc-sec x)) | |
943 | (math-make-sdev (math-div-float '(float 1 0) (cdr sc)) | |
944 | (math-div-float | |
945 | (math-mul xs (car sc)) | |
946 | (math-sqr (cdr sc))))))) | |
947 | (math-make-sdev (calcFunc-sec (nth 1 x)) | |
948 | (math-div | |
949 | (math-mul (nth 2 x) | |
950 | (calcFunc-sin (nth 1 x))) | |
951 | (math-sqr (calcFunc-cos (nth 1 x))))))) | |
952 | ((and (eq (car x) 'intv) | |
953 | (math-intv-constp x)) | |
954 | (math-with-extra-prec 2 | |
955 | (let* ((xx (math-to-radians (math-float x))) | |
956 | (na (math-floor (math-div (math-sub (nth 2 xx) | |
957 | (math-pi-over-2)) | |
958 | (math-pi)))) | |
959 | (nb (math-floor (math-div (math-sub (nth 3 xx) | |
960 | (math-pi-over-2)) | |
961 | (math-pi)))) | |
962 | (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2)))) | |
963 | (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2)))) | |
964 | (span (math-sub nbb naa))) | |
965 | (if (not (equal na nb)) | |
966 | '(intv 3 (neg (var inf var-inf)) (var inf var-inf)) | |
967 | (let ((int (math-sort-intv (nth 1 x) | |
968 | (math-sec-raw (nth 2 xx)) | |
969 | (math-sec-raw (nth 3 xx))))) | |
970 | (if (eq span 1) | |
971 | (if (math-evenp (math-div (math-add naa 1) 2)) | |
972 | (math-make-intv (logior (nth 1 int) 2) | |
973 | 1 | |
974 | (nth 3 int)) | |
975 | (math-make-intv (logior (nth 1 int) 1) | |
976 | (nth 2 int) | |
977 | -1)) | |
978 | int)))))) | |
979 | ((equal x '(var nan var-nan)) | |
980 | x) | |
981 | (t (calc-record-why 'scalarp x) | |
982 | (list 'calcFunc-sec x)))) | |
983 | ||
984 | (defun calcFunc-csc (x) | |
985 | (cond ((and (integerp x) | |
986 | (eq calc-angle-mode 'deg) | |
987 | (= (% (- x 90) 180) 0)) | |
988 | (if (= (% (- x 90) 360) 0) | |
989 | 1 | |
990 | -1)) | |
991 | ((Math-scalarp x) | |
992 | (math-with-extra-prec 2 | |
993 | (math-csc-raw (math-to-radians (math-float x))))) | |
994 | ((eq (car x) 'sdev) | |
995 | (if (math-constp x) | |
996 | (math-with-extra-prec 2 | |
997 | (let* ((xx (math-to-radians (math-float (nth 1 x)))) | |
998 | (xs (math-to-radians (math-float (nth 2 x)))) | |
999 | (sc (math-sin-cos-raw xx))) | |
1000 | (if (and (math-zerop (car sc)) | |
1001 | (not calc-infinite-mode)) | |
1002 | (progn | |
1003 | (calc-record-why "*Division by zero") | |
1004 | (list 'calcFunc-csc x)) | |
1005 | (math-make-sdev (math-div-float '(float 1 0) (car sc)) | |
1006 | (math-div-float | |
1007 | (math-mul xs (cdr sc)) | |
1008 | (math-sqr (car sc))))))) | |
1009 | (math-make-sdev (calcFunc-csc (nth 1 x)) | |
1010 | (math-div | |
1011 | (math-mul (nth 2 x) | |
1012 | (calcFunc-cos (nth 1 x))) | |
1013 | (math-sqr (calcFunc-sin (nth 1 x))))))) | |
1014 | ((and (eq (car x) 'intv) | |
1015 | (math-intv-constp x)) | |
1016 | (math-with-extra-prec 2 | |
1017 | (let* ((xx (math-to-radians (math-float x))) | |
1018 | (na (math-floor (math-div (nth 2 xx) (math-pi)))) | |
1019 | (nb (math-floor (math-div (nth 3 xx) (math-pi)))) | |
1020 | (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2)))) | |
1021 | (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2)))) | |
1022 | (span (math-sub nbb naa))) | |
1023 | (if (not (equal na nb)) | |
1024 | '(intv 3 (neg (var inf var-inf)) (var inf var-inf)) | |
1025 | (let ((int (math-sort-intv (nth 1 x) | |
1026 | (math-csc-raw (nth 2 xx)) | |
1027 | (math-csc-raw (nth 3 xx))))) | |
1028 | (if (eq span 1) | |
1029 | (if (math-evenp (math-div naa 2)) | |
1030 | (math-make-intv (logior (nth 1 int) 2) | |
1031 | 1 | |
1032 | (nth 3 int)) | |
1033 | (math-make-intv (logior (nth 1 int) 1) | |
1034 | (nth 2 int) | |
1035 | -1)) | |
1036 | int)))))) | |
1037 | ((equal x '(var nan var-nan)) | |
1038 | x) | |
1039 | (t (calc-record-why 'scalarp x) | |
1040 | (list 'calcFunc-csc x)))) | |
1041 | ||
1042 | (defun calcFunc-cot (x) ; [N N] [Public] | |
1043 | (cond ((and (integerp x) | |
1044 | (if (eq calc-angle-mode 'deg) | |
1045 | (= (% (- x 90) 180) 0) | |
1046 | (= x 0))) | |
1047 | 0) | |
1048 | ((Math-scalarp x) | |
1049 | (math-with-extra-prec 2 | |
1050 | (math-cot-raw (math-to-radians (math-float x))))) | |
1051 | ((eq (car x) 'sdev) | |
1052 | (if (math-constp x) | |
1053 | (math-with-extra-prec 2 | |
1054 | (let* ((xx (math-to-radians (math-float (nth 1 x)))) | |
1055 | (xs (math-to-radians (math-float (nth 2 x)))) | |
1056 | (sc (math-sin-cos-raw xx))) | |
1057 | (if (and (math-zerop (car sc)) (not calc-infinite-mode)) | |
1058 | (progn | |
1059 | (calc-record-why "*Division by zero") | |
1060 | (list 'calcFunc-cot x)) | |
1061 | (math-make-sdev (math-div-float (cdr sc) (car sc)) | |
1062 | (math-div-float xs (math-sqr (car sc))))))) | |
1063 | (math-make-sdev (calcFunc-cot (nth 1 x)) | |
1064 | (math-div (nth 2 x) | |
1065 | (math-sqr (calcFunc-sin (nth 1 x))))))) | |
1066 | ((and (eq (car x) 'intv) (math-intv-constp x)) | |
1067 | (or (math-with-extra-prec 2 | |
1068 | (let* ((xx (math-to-radians (math-float x))) | |
1069 | (na (math-floor (math-div (nth 2 xx) (math-pi)))) | |
eb1ef455 | 1070 | (nb (math-floor (math-div (nth 3 xx) (math-pi))))) |
f53e6c20 JB |
1071 | (and (equal na nb) |
1072 | (math-sort-intv (nth 1 x) | |
1073 | (math-cot-raw (nth 2 xx)) | |
eb1ef455 | 1074 | (math-cot-raw (nth 3 xx)))))) |
f53e6c20 JB |
1075 | '(intv 3 (neg (var inf var-inf)) (var inf var-inf)))) |
1076 | ((equal x '(var nan var-nan)) | |
1077 | x) | |
1078 | (t (calc-record-why 'scalarp x) | |
1079 | (list 'calcFunc-cot x)))) | |
1080 | ||
b0d9db86 | 1081 | (defun math-sin-raw (x &optional orgx) ; [N N] |
136211a9 EZ |
1082 | (cond ((eq (car x) 'cplx) |
1083 | (let* ((expx (math-exp-raw (nth 2 x))) | |
1084 | (expmx (math-div-float '(float 1 0) expx)) | |
1085 | (sc (math-sin-cos-raw (nth 1 x)))) | |
1086 | (list 'cplx | |
1087 | (math-mul-float (car sc) | |
1088 | (math-mul-float (math-add-float expx expmx) | |
1089 | '(float 5 -1))) | |
1090 | (math-mul-float (cdr sc) | |
1091 | (math-mul-float (math-sub-float expx expmx) | |
1092 | '(float 5 -1)))))) | |
1093 | ((eq (car x) 'polar) | |
1094 | (math-polar (math-sin-raw (math-complex x)))) | |
1095 | ((Math-integer-negp (nth 1 x)) | |
b0d9db86 | 1096 | (math-neg-float (math-sin-raw (math-neg-float x) (if orgx orgx x)))) |
136211a9 | 1097 | ((math-lessp-float '(float 7 0) x) ; avoid inf loops due to roundoff |
b0d9db86 JB |
1098 | (math-sin-raw (math-mod x (math-two-pi)) (if orgx orgx x))) |
1099 | (t (math-sin-raw-2 x (if orgx orgx x))))) | |
136211a9 EZ |
1100 | |
1101 | (defun math-cos-raw (x) ; [N N] | |
1102 | (if (eq (car-safe x) 'polar) | |
1103 | (math-polar (math-cos-raw (math-complex x))) | |
b0d9db86 | 1104 | (math-sin-raw (math-sub (math-pi-over-2) x) x))) |
136211a9 | 1105 | |
f53e6c20 JB |
1106 | (defun math-sec-raw (x) ; [N N] |
1107 | (cond ((eq (car x) 'cplx) | |
1108 | (let* ((x (math-mul x '(float 1 0))) | |
1109 | (expx (math-exp-raw (nth 2 x))) | |
1110 | (expmx (math-div-float '(float 1 0) expx)) | |
1111 | (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1))) | |
1112 | (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1))) | |
1113 | (sc (math-sin-cos-raw (nth 1 x))) | |
1114 | (d (math-add-float | |
1115 | (math-mul-float (math-sqr (car sc)) | |
1116 | (math-sqr sh)) | |
1117 | (math-mul-float (math-sqr (cdr sc)) | |
1118 | (math-sqr ch))))) | |
1119 | (and (not (eq (nth 1 d) 0)) | |
1120 | (list 'cplx | |
1121 | (math-div-float (math-mul-float (cdr sc) ch) d) | |
1122 | (math-div-float (math-mul-float (car sc) sh) d))))) | |
1123 | ((eq (car x) 'polar) | |
1124 | (math-polar (math-sec-raw (math-complex x)))) | |
1125 | (t | |
1126 | (let ((cs (math-cos-raw x))) | |
1127 | (if (eq cs 0) | |
1128 | (math-div 1 0) | |
1129 | (math-div-float '(float 1 0) cs)))))) | |
1130 | ||
1131 | (defun math-csc-raw (x) ; [N N] | |
1132 | (cond ((eq (car x) 'cplx) | |
1133 | (let* ((x (math-mul x '(float 1 0))) | |
1134 | (expx (math-exp-raw (nth 2 x))) | |
1135 | (expmx (math-div-float '(float 1 0) expx)) | |
1136 | (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1))) | |
1137 | (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1))) | |
1138 | (sc (math-sin-cos-raw (nth 1 x))) | |
1139 | (d (math-add-float | |
1140 | (math-mul-float (math-sqr (car sc)) | |
1141 | (math-sqr ch)) | |
1142 | (math-mul-float (math-sqr (cdr sc)) | |
1143 | (math-sqr sh))))) | |
1144 | (and (not (eq (nth 1 d) 0)) | |
1145 | (list 'cplx | |
1146 | (math-div-float (math-mul-float (car sc) ch) d) | |
1147 | (math-div-float (math-mul-float (cdr sc) sh) d))))) | |
1148 | ((eq (car x) 'polar) | |
6ec30302 | 1149 | (math-polar (math-csc-raw (math-complex x)))) |
f53e6c20 JB |
1150 | (t |
1151 | (let ((sn (math-sin-raw x))) | |
1152 | (if (eq sn 0) | |
1153 | (math-div 1 0) | |
1154 | (math-div-float '(float 1 0) sn)))))) | |
1155 | ||
1156 | (defun math-cot-raw (x) ; [N N] | |
1157 | (cond ((eq (car x) 'cplx) | |
1158 | (let* ((x (math-mul x '(float 1 0))) | |
1159 | (expx (math-exp-raw (nth 2 x))) | |
1160 | (expmx (math-div-float '(float 1 0) expx)) | |
1161 | (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1))) | |
1162 | (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1))) | |
1163 | (sc (math-sin-cos-raw (nth 1 x))) | |
1164 | (d (math-add-float | |
1165 | (math-sqr (car sc)) | |
1166 | (math-sqr sh)))) | |
1167 | (and (not (eq (nth 1 d) 0)) | |
1168 | (list 'cplx | |
1169 | (math-div-float | |
1170 | (math-mul-float (car sc) (cdr sc)) | |
1171 | d) | |
1172 | (math-neg | |
1173 | (math-div-float | |
1174 | (math-mul-float sh ch) | |
1175 | d)))))) | |
1176 | ((eq (car x) 'polar) | |
1177 | (math-polar (math-cot-raw (math-complex x)))) | |
1178 | (t | |
1179 | (let ((sc (math-sin-cos-raw x))) | |
1180 | (if (eq (nth 1 (car sc)) 0) | |
1181 | (math-div (cdr sc) 0) | |
1182 | (math-div-float (cdr sc) (car sc))))))) | |
1183 | ||
1184 | ||
136211a9 EZ |
1185 | ;;; This could use a smarter method: Reduce x as in math-sin-raw, then |
1186 | ;;; compute either sin(x) or cos(x), whichever is smaller, and compute | |
1187 | ;;; the other using the identity sin(x)^2 + cos(x)^2 = 1. | |
1188 | (defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x)) | |
491c3062 | 1189 | (cons (math-sin-raw x) (math-cos-raw x))) |
136211a9 EZ |
1190 | |
1191 | (defun math-tan-raw (x) ; [N N] | |
1192 | (cond ((eq (car x) 'cplx) | |
1193 | (let* ((x (math-mul x '(float 2 0))) | |
1194 | (expx (math-exp-raw (nth 2 x))) | |
1195 | (expmx (math-div-float '(float 1 0) expx)) | |
1196 | (sc (math-sin-cos-raw (nth 1 x))) | |
1197 | (d (math-add-float (cdr sc) | |
1198 | (math-mul-float (math-add-float expx expmx) | |
1199 | '(float 5 -1))))) | |
1200 | (and (not (eq (nth 1 d) 0)) | |
1201 | (list 'cplx | |
1202 | (math-div-float (car sc) d) | |
1203 | (math-div-float (math-mul-float (math-sub-float expx | |
1204 | expmx) | |
1205 | '(float 5 -1)) d))))) | |
1206 | ((eq (car x) 'polar) | |
1207 | (math-polar (math-tan-raw (math-complex x)))) | |
1208 | (t | |
1209 | (let ((sc (math-sin-cos-raw x))) | |
1210 | (if (eq (nth 1 (cdr sc)) 0) | |
1211 | (math-div (car sc) 0) | |
491c3062 | 1212 | (math-div-float (car sc) (cdr sc))))))) |
136211a9 EZ |
1213 | |
1214 | (defun math-sin-raw-2 (x orgx) ; This avoids poss of inf recursion. [F F] | |
1215 | (let ((xmpo2 (math-sub-float (math-pi-over-2) x))) | |
1216 | (cond ((Math-integer-negp (nth 1 xmpo2)) | |
1217 | (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi)) | |
1218 | orgx))) | |
1219 | ((math-lessp-float (math-pi-over-4) x) | |
1220 | (math-cos-raw-2 xmpo2 orgx)) | |
1221 | ((math-lessp-float x (math-neg (math-pi-over-4))) | |
1222 | (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx))) | |
b0d9db86 JB |
1223 | ((math-with-extra-prec -1 (math-nearly-zerop-float x orgx)) |
1224 | '(float 0 0)) | |
17cc361e | 1225 | ((math-use-emacs-fn 'sin x)) |
136211a9 | 1226 | (calc-symbolic-mode (signal 'inexact-result nil)) |
491c3062 | 1227 | (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x))))))) |
136211a9 EZ |
1228 | |
1229 | (defun math-cos-raw-2 (x orgx) ; [F F] | |
b0d9db86 JB |
1230 | (cond ((math-with-extra-prec -1 (math-nearly-zerop-float x orgx)) |
1231 | '(float 1 0)) | |
17cc361e | 1232 | ((math-use-emacs-fn 'cos x)) |
136211a9 EZ |
1233 | (calc-symbolic-mode (signal 'inexact-result nil)) |
1234 | (t (let ((xnegsqr (math-neg-float (math-sqr-float x)))) | |
1235 | (math-sin-series | |
1236 | (math-add-float '(float 1 0) | |
1237 | (math-mul-float xnegsqr '(float 5 -1))) | |
491c3062 | 1238 | 24 5 xnegsqr xnegsqr))))) |
136211a9 EZ |
1239 | |
1240 | (defun math-sin-series (sum nfac n x xnegsqr) | |
1241 | (math-working "sin" sum) | |
1242 | (let* ((nextx (math-mul-float x xnegsqr)) | |
1243 | (nextsum (math-add-float sum (math-div-float nextx | |
1244 | (math-float nfac))))) | |
1245 | (if (math-nearly-equal-float sum nextsum) | |
1246 | sum | |
1247 | (math-sin-series nextsum (math-mul nfac (* n (1+ n))) | |
491c3062 | 1248 | (+ n 2) nextx xnegsqr)))) |
136211a9 EZ |
1249 | |
1250 | ||
1251 | ;;; Inverse sine, cosine, tangent. | |
1252 | ||
1253 | (defun calcFunc-arcsin (x) ; [N N] [Public] | |
1254 | (cond ((eq x 0) 0) | |
1255 | ((and (eq x 1) (eq calc-angle-mode 'deg)) 90) | |
1256 | ((and (eq x -1) (eq calc-angle-mode 'deg)) -90) | |
1257 | (calc-symbolic-mode (signal 'inexact-result nil)) | |
1258 | ((Math-numberp x) | |
1259 | (math-with-extra-prec 2 | |
1260 | (math-from-radians (math-arcsin-raw (math-float x))))) | |
1261 | ((eq (car x) 'sdev) | |
1262 | (math-make-sdev (calcFunc-arcsin (nth 1 x)) | |
1263 | (math-from-radians | |
1264 | (math-div (nth 2 x) | |
1265 | (math-sqrt | |
1266 | (math-sub 1 (math-sqr (nth 1 x)))))))) | |
1267 | ((eq (car x) 'intv) | |
1268 | (math-sort-intv (nth 1 x) | |
1269 | (calcFunc-arcsin (nth 2 x)) | |
1270 | (calcFunc-arcsin (nth 3 x)))) | |
1271 | ((equal x '(var nan var-nan)) | |
1272 | x) | |
1273 | (t (calc-record-why 'numberp x) | |
491c3062 | 1274 | (list 'calcFunc-arcsin x)))) |
136211a9 EZ |
1275 | |
1276 | (defun calcFunc-arccos (x) ; [N N] [Public] | |
1277 | (cond ((eq x 1) 0) | |
1278 | ((and (eq x 0) (eq calc-angle-mode 'deg)) 90) | |
1279 | ((and (eq x -1) (eq calc-angle-mode 'deg)) 180) | |
1280 | (calc-symbolic-mode (signal 'inexact-result nil)) | |
1281 | ((Math-numberp x) | |
1282 | (math-with-extra-prec 2 | |
1283 | (math-from-radians (math-arccos-raw (math-float x))))) | |
1284 | ((eq (car x) 'sdev) | |
1285 | (math-make-sdev (calcFunc-arccos (nth 1 x)) | |
1286 | (math-from-radians | |
1287 | (math-div (nth 2 x) | |
1288 | (math-sqrt | |
1289 | (math-sub 1 (math-sqr (nth 1 x)))))))) | |
1290 | ((eq (car x) 'intv) | |
1291 | (math-sort-intv (nth 1 x) | |
1292 | (calcFunc-arccos (nth 2 x)) | |
1293 | (calcFunc-arccos (nth 3 x)))) | |
1294 | ((equal x '(var nan var-nan)) | |
1295 | x) | |
1296 | (t (calc-record-why 'numberp x) | |
491c3062 | 1297 | (list 'calcFunc-arccos x)))) |
136211a9 EZ |
1298 | |
1299 | (defun calcFunc-arctan (x) ; [N N] [Public] | |
1300 | (cond ((eq x 0) 0) | |
1301 | ((and (eq x 1) (eq calc-angle-mode 'deg)) 45) | |
1302 | ((and (eq x -1) (eq calc-angle-mode 'deg)) -45) | |
1303 | ((Math-numberp x) | |
1304 | (math-with-extra-prec 2 | |
1305 | (math-from-radians (math-arctan-raw (math-float x))))) | |
1306 | ((eq (car x) 'sdev) | |
1307 | (math-make-sdev (calcFunc-arctan (nth 1 x)) | |
1308 | (math-from-radians | |
1309 | (math-div (nth 2 x) | |
1310 | (math-add 1 (math-sqr (nth 1 x))))))) | |
1311 | ((eq (car x) 'intv) | |
1312 | (math-sort-intv (nth 1 x) | |
1313 | (calcFunc-arctan (nth 2 x)) | |
1314 | (calcFunc-arctan (nth 3 x)))) | |
1315 | ((equal x '(var inf var-inf)) | |
1316 | (math-quarter-circle t)) | |
1317 | ((equal x '(neg (var inf var-inf))) | |
1318 | (math-neg (math-quarter-circle t))) | |
1319 | ((equal x '(var nan var-nan)) | |
1320 | x) | |
1321 | (t (calc-record-why 'numberp x) | |
491c3062 | 1322 | (list 'calcFunc-arctan x)))) |
136211a9 EZ |
1323 | |
1324 | (defun math-arcsin-raw (x) ; [N N] | |
1325 | (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x))))) | |
1326 | (if (or (memq (car x) '(cplx polar)) | |
1327 | (memq (car a) '(cplx polar))) | |
1328 | (math-with-extra-prec 2 ; use extra precision for difficult case | |
1329 | (math-mul '(cplx 0 -1) | |
1330 | (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a)))) | |
491c3062 | 1331 | (math-arctan2-raw x a)))) |
136211a9 EZ |
1332 | |
1333 | (defun math-arccos-raw (x) ; [N N] | |
491c3062 | 1334 | (math-sub (math-pi-over-2) (math-arcsin-raw x))) |
136211a9 EZ |
1335 | |
1336 | (defun math-arctan-raw (x) ; [N N] | |
1337 | (cond ((memq (car x) '(cplx polar)) | |
1338 | (math-with-extra-prec 2 ; extra-extra | |
1339 | (math-div (math-sub | |
1340 | (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x))) | |
1341 | (math-ln-raw (math-add 1 (math-mul '(cplx 0 -1) x)))) | |
1342 | '(cplx 0 2)))) | |
1343 | ((Math-integer-negp (nth 1 x)) | |
1344 | (math-neg-float (math-arctan-raw (math-neg-float x)))) | |
1345 | ((math-zerop x) x) | |
17cc361e | 1346 | ((math-use-emacs-fn 'atan x)) |
136211a9 EZ |
1347 | (calc-symbolic-mode (signal 'inexact-result nil)) |
1348 | ((math-equal-int x 1) (math-pi-over-4)) | |
1349 | ((math-equal-int x -1) (math-neg (math-pi-over-4))) | |
1350 | ((math-lessp-float '(float 414214 -6) x) ; if x > sqrt(2) - 1, reduce | |
1351 | (if (math-lessp-float '(float 1 0) x) | |
1352 | (math-sub-float (math-mul-float (math-pi) '(float 5 -1)) | |
1353 | (math-arctan-raw (math-div-float '(float 1 0) x))) | |
1354 | (math-sub-float (math-mul-float (math-pi) '(float 25 -2)) | |
1355 | (math-arctan-raw (math-div-float | |
1356 | (math-sub-float '(float 1 0) x) | |
1357 | (math-add-float '(float 1 0) | |
1358 | x)))))) | |
491c3062 | 1359 | (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x)))))) |
136211a9 EZ |
1360 | |
1361 | (defun math-arctan-series (sum n x xnegsqr) | |
1362 | (math-working "arctan" sum) | |
1363 | (let* ((nextx (math-mul-float x xnegsqr)) | |
1364 | (nextsum (math-add-float sum (math-div-float nextx (math-float n))))) | |
1365 | (if (math-nearly-equal-float sum nextsum) | |
1366 | sum | |
491c3062 | 1367 | (math-arctan-series nextsum (+ n 2) nextx xnegsqr)))) |
136211a9 EZ |
1368 | |
1369 | (defun calcFunc-arctan2 (y x) ; [F R R] [Public] | |
1370 | (if (Math-anglep y) | |
1371 | (if (Math-anglep x) | |
1372 | (math-with-extra-prec 2 | |
1373 | (math-from-radians (math-arctan2-raw (math-float y) | |
1374 | (math-float x)))) | |
1375 | (calc-record-why 'anglep x) | |
1376 | (list 'calcFunc-arctan2 y x)) | |
1377 | (if (and (or (math-infinitep x) (math-anglep x)) | |
1378 | (or (math-infinitep y) (math-anglep y))) | |
1379 | (progn | |
1380 | (if (math-posp x) | |
1381 | (setq x 1) | |
1382 | (if (math-negp x) | |
1383 | (setq x -1) | |
1384 | (or (math-zerop x) | |
1385 | (setq x nil)))) | |
1386 | (if (math-posp y) | |
1387 | (setq y 1) | |
1388 | (if (math-negp y) | |
1389 | (setq y -1) | |
1390 | (or (math-zerop y) | |
1391 | (setq y nil)))) | |
1392 | (if (and y x) | |
1393 | (calcFunc-arctan2 y x) | |
1394 | '(var nan var-nan))) | |
1395 | (calc-record-why 'anglep y) | |
491c3062 | 1396 | (list 'calcFunc-arctan2 y x)))) |
136211a9 EZ |
1397 | |
1398 | (defun math-arctan2-raw (y x) ; [F R R] | |
1399 | (cond ((math-zerop y) | |
1400 | (if (math-negp x) (math-pi) | |
1401 | (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0))) | |
1402 | ((math-zerop x) | |
1403 | (if (math-posp y) | |
1404 | (math-pi-over-2) | |
1405 | (math-neg (math-pi-over-2)))) | |
1406 | ((math-posp x) | |
1407 | (math-arctan-raw (math-div-float y x))) | |
1408 | ((math-posp y) | |
1409 | (math-add-float (math-arctan-raw (math-div-float y x)) | |
1410 | (math-pi))) | |
1411 | (t | |
1412 | (math-sub-float (math-arctan-raw (math-div-float y x)) | |
491c3062 | 1413 | (math-pi))))) |
136211a9 EZ |
1414 | |
1415 | (defun calcFunc-arcsincos (x) ; [V N] [Public] | |
1416 | (if (and (Math-vectorp x) | |
1417 | (= (length x) 3)) | |
1418 | (calcFunc-arctan2 (nth 2 x) (nth 1 x)) | |
491c3062 | 1419 | (math-reject-arg x "*Two-element vector expected"))) |
136211a9 EZ |
1420 | |
1421 | ||
1422 | ||
1423 | ;;; Exponential function. | |
1424 | ||
1425 | (defun calcFunc-exp (x) ; [N N] [Public] | |
1426 | (cond ((eq x 0) 1) | |
1427 | ((and (memq x '(1 -1)) calc-symbolic-mode) | |
1428 | (if (eq x 1) '(var e var-e) (math-div 1 '(var e var-e)))) | |
1429 | ((Math-numberp x) | |
1430 | (math-with-extra-prec 2 (math-exp-raw (math-float x)))) | |
1431 | ((eq (car-safe x) 'sdev) | |
1432 | (let ((ex (calcFunc-exp (nth 1 x)))) | |
1433 | (math-make-sdev ex (math-mul (nth 2 x) ex)))) | |
1434 | ((eq (car-safe x) 'intv) | |
1435 | (math-make-intv (nth 1 x) (calcFunc-exp (nth 2 x)) | |
1436 | (calcFunc-exp (nth 3 x)))) | |
1437 | ((equal x '(var inf var-inf)) | |
1438 | x) | |
1439 | ((equal x '(neg (var inf var-inf))) | |
1440 | 0) | |
1441 | ((equal x '(var nan var-nan)) | |
1442 | x) | |
1443 | (t (calc-record-why 'numberp x) | |
491c3062 | 1444 | (list 'calcFunc-exp x)))) |
136211a9 EZ |
1445 | |
1446 | (defun calcFunc-expm1 (x) ; [N N] [Public] | |
1447 | (cond ((eq x 0) 0) | |
1448 | ((math-zerop x) '(float 0 0)) | |
1449 | (calc-symbolic-mode (signal 'inexact-result nil)) | |
1450 | ((Math-numberp x) | |
1451 | (math-with-extra-prec 2 | |
1452 | (let ((x (math-float x))) | |
1453 | (if (and (eq (car x) 'float) | |
1454 | (math-lessp-float x '(float 1 0)) | |
1455 | (math-lessp-float '(float -1 0) x)) | |
1456 | (math-exp-minus-1-raw x) | |
1457 | (math-add (math-exp-raw x) -1))))) | |
1458 | ((eq (car-safe x) 'sdev) | |
1459 | (if (math-constp x) | |
1460 | (let ((ex (calcFunc-expm1 (nth 1 x)))) | |
1461 | (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1)))) | |
1462 | (math-make-sdev (calcFunc-expm1 (nth 1 x)) | |
1463 | (math-mul (nth 2 x) (calcFunc-exp (nth 1 x)))))) | |
1464 | ((eq (car-safe x) 'intv) | |
1465 | (math-make-intv (nth 1 x) | |
1466 | (calcFunc-expm1 (nth 2 x)) | |
1467 | (calcFunc-expm1 (nth 3 x)))) | |
1468 | ((equal x '(var inf var-inf)) | |
1469 | x) | |
1470 | ((equal x '(neg (var inf var-inf))) | |
1471 | -1) | |
1472 | ((equal x '(var nan var-nan)) | |
1473 | x) | |
1474 | (t (calc-record-why 'numberp x) | |
491c3062 | 1475 | (list 'calcFunc-expm1 x)))) |
136211a9 EZ |
1476 | |
1477 | (defun calcFunc-exp10 (x) ; [N N] [Public] | |
1478 | (if (eq x 0) | |
1479 | 1 | |
491c3062 | 1480 | (math-pow '(float 1 1) x))) |
136211a9 EZ |
1481 | |
1482 | (defun math-exp-raw (x) ; [N N] | |
1483 | (cond ((math-zerop x) '(float 1 0)) | |
1484 | (calc-symbolic-mode (signal 'inexact-result nil)) | |
1485 | ((eq (car x) 'cplx) | |
1486 | (let ((expx (math-exp-raw (nth 1 x))) | |
1487 | (sc (math-sin-cos-raw (nth 2 x)))) | |
1488 | (list 'cplx | |
1489 | (math-mul-float expx (cdr sc)) | |
1490 | (math-mul-float expx (car sc))))) | |
1491 | ((eq (car x) 'polar) | |
1492 | (let ((xc (math-complex x))) | |
1493 | (list 'polar | |
1494 | (math-exp-raw (nth 1 xc)) | |
1495 | (math-from-radians (nth 2 xc))))) | |
178b8baf | 1496 | ((math-use-emacs-fn 'exp x)) |
136211a9 EZ |
1497 | ((or (math-lessp-float '(float 5 -1) x) |
1498 | (math-lessp-float x '(float -5 -1))) | |
1499 | (if (math-lessp-float '(float 921035 1) x) | |
1500 | (math-overflow) | |
1501 | (if (math-lessp-float x '(float -921035 1)) | |
1502 | (math-underflow))) | |
1503 | (let* ((two-x (math-mul-float x '(float 2 0))) | |
1504 | (hint (math-scale-int (nth 1 two-x) (nth 2 two-x))) | |
1505 | (hfrac (math-sub-float x (math-mul-float (math-float hint) | |
1506 | '(float 5 -1))))) | |
1507 | (math-mul-float (math-ipow (math-sqrt-e) hint) | |
1508 | (math-add-float '(float 1 0) | |
1509 | (math-exp-minus-1-raw hfrac))))) | |
491c3062 | 1510 | (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x))))) |
136211a9 EZ |
1511 | |
1512 | (defun math-exp-minus-1-raw (x) ; [F F] | |
491c3062 | 1513 | (math-exp-series x 2 3 x x)) |
136211a9 EZ |
1514 | |
1515 | (defun math-exp-series (sum nfac n xpow x) | |
1516 | (math-working "exp" sum) | |
1517 | (let* ((nextx (math-mul-float xpow x)) | |
1518 | (nextsum (math-add-float sum (math-div-float nextx | |
1519 | (math-float nfac))))) | |
1520 | (if (math-nearly-equal-float sum nextsum) | |
1521 | sum | |
491c3062 | 1522 | (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x)))) |
136211a9 EZ |
1523 | |
1524 | ||
1525 | ||
1526 | ;;; Logarithms. | |
1527 | ||
1528 | (defun calcFunc-ln (x) ; [N N] [Public] | |
1529 | (cond ((math-zerop x) | |
1530 | (if calc-infinite-mode | |
1531 | '(neg (var inf var-inf)) | |
1532 | (math-reject-arg x "*Logarithm of zero"))) | |
1533 | ((eq x 1) 0) | |
1534 | ((Math-numberp x) | |
1535 | (math-with-extra-prec 2 (math-ln-raw (math-float x)))) | |
1536 | ((eq (car-safe x) 'sdev) | |
1537 | (math-make-sdev (calcFunc-ln (nth 1 x)) | |
1538 | (math-div (nth 2 x) (nth 1 x)))) | |
1539 | ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) | |
1540 | (Math-zerop (nth 2 x)) | |
1541 | (not (math-intv-constp x)))) | |
1542 | (let ((calc-infinite-mode t)) | |
1543 | (math-make-intv (nth 1 x) (calcFunc-ln (nth 2 x)) | |
1544 | (calcFunc-ln (nth 3 x))))) | |
1545 | ((equal x '(var e var-e)) | |
1546 | 1) | |
1547 | ((and (eq (car-safe x) '^) | |
1548 | (equal (nth 1 x) '(var e var-e)) | |
1549 | (math-known-realp (nth 2 x))) | |
1550 | (nth 2 x)) | |
1551 | ((math-infinitep x) | |
1552 | (if (equal x '(var nan var-nan)) | |
1553 | x | |
1554 | '(var inf var-inf))) | |
1555 | (t (calc-record-why 'numberp x) | |
491c3062 | 1556 | (list 'calcFunc-ln x)))) |
136211a9 EZ |
1557 | |
1558 | (defun calcFunc-log10 (x) ; [N N] [Public] | |
1559 | (cond ((math-equal-int x 1) | |
1560 | (if (math-floatp x) '(float 0 0) 0)) | |
1561 | ((and (Math-integerp x) | |
1562 | (math-posp x) | |
1563 | (let ((res (math-integer-log x 10))) | |
1564 | (and (car res) | |
1565 | (setq x (cdr res))))) | |
1566 | x) | |
1567 | ((and (eq (car-safe x) 'frac) | |
1568 | (eq (nth 1 x) 1) | |
1569 | (let ((res (math-integer-log (nth 2 x) 10))) | |
1570 | (and (car res) | |
1571 | (setq x (- (cdr res)))))) | |
1572 | x) | |
1573 | ((math-zerop x) | |
1574 | (if calc-infinite-mode | |
1575 | '(neg (var inf var-inf)) | |
1576 | (math-reject-arg x "*Logarithm of zero"))) | |
1577 | (calc-symbolic-mode (signal 'inexact-result nil)) | |
1578 | ((Math-numberp x) | |
1579 | (math-with-extra-prec 2 | |
1580 | (let ((xf (math-float x))) | |
1581 | (if (eq (nth 1 xf) 0) | |
1582 | (math-reject-arg x "*Logarithm of zero")) | |
1583 | (if (Math-integer-posp (nth 1 xf)) | |
1584 | (if (eq (nth 1 xf) 1) ; log10(1*10^n) = n | |
1585 | (math-float (nth 2 xf)) | |
1586 | (let ((xdigs (1- (math-numdigs (nth 1 xf))))) | |
1587 | (math-add-float | |
1588 | (math-div-float (math-ln-raw-2 | |
1589 | (list 'float (nth 1 xf) (- xdigs))) | |
1590 | (math-ln-10)) | |
1591 | (math-float (+ (nth 2 xf) xdigs))))) | |
1592 | (math-div (calcFunc-ln xf) (math-ln-10)))))) | |
1593 | ((eq (car-safe x) 'sdev) | |
1594 | (math-make-sdev (calcFunc-log10 (nth 1 x)) | |
1595 | (math-div (nth 2 x) | |
1596 | (math-mul (nth 1 x) (math-ln-10))))) | |
1597 | ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) | |
1598 | (not (math-intv-constp x)))) | |
1599 | (math-make-intv (nth 1 x) | |
1600 | (calcFunc-log10 (nth 2 x)) | |
1601 | (calcFunc-log10 (nth 3 x)))) | |
1602 | ((math-infinitep x) | |
1603 | (if (equal x '(var nan var-nan)) | |
1604 | x | |
1605 | '(var inf var-inf))) | |
1606 | (t (calc-record-why 'numberp x) | |
491c3062 | 1607 | (list 'calcFunc-log10 x)))) |
136211a9 EZ |
1608 | |
1609 | (defun calcFunc-log (x &optional b) ; [N N N] [Public] | |
1610 | (cond ((or (null b) (equal b '(var e var-e))) | |
1611 | (math-normalize (list 'calcFunc-ln x))) | |
1612 | ((or (eq b 10) (equal b '(float 1 1))) | |
1613 | (math-normalize (list 'calcFunc-log10 x))) | |
1614 | ((math-zerop x) | |
1615 | (if calc-infinite-mode | |
1616 | (math-div (calcFunc-ln x) (calcFunc-ln b)) | |
1617 | (math-reject-arg x "*Logarithm of zero"))) | |
1618 | ((math-zerop b) | |
1619 | (if calc-infinite-mode | |
1620 | (math-div (calcFunc-ln x) (calcFunc-ln b)) | |
1621 | (math-reject-arg b "*Logarithm of zero"))) | |
1622 | ((math-equal-int b 1) | |
1623 | (if calc-infinite-mode | |
1624 | (math-div (calcFunc-ln x) 0) | |
1625 | (math-reject-arg b "*Logarithm base one"))) | |
1626 | ((math-equal-int x 1) | |
86498823 | 1627 | (if (math-floatp b) '(float 0 0) 0)) |
136211a9 EZ |
1628 | ((and (Math-ratp x) (Math-ratp b) |
1629 | (math-posp x) (math-posp b) | |
1630 | (let* ((sign 1) (inv nil) | |
1631 | (xx (if (Math-lessp 1 x) | |
1632 | x | |
1633 | (setq sign -1) | |
1634 | (math-div 1 x))) | |
1635 | (bb (if (Math-lessp 1 b) | |
1636 | b | |
1637 | (setq sign (- sign)) | |
1638 | (math-div 1 b))) | |
1639 | (res (if (Math-lessp xx bb) | |
1640 | (setq inv (math-integer-log bb xx)) | |
1641 | (math-integer-log xx bb)))) | |
1642 | (and (car res) | |
1643 | (setq x (if inv | |
1644 | (math-div 1 (* sign (cdr res))) | |
1645 | (* sign (cdr res))))))) | |
1646 | x) | |
1647 | (calc-symbolic-mode (signal 'inexact-result nil)) | |
1648 | ((and (Math-numberp x) (Math-numberp b)) | |
1649 | (math-with-extra-prec 2 | |
1650 | (math-div (math-ln-raw (math-float x)) | |
1651 | (math-log-base-raw b)))) | |
1652 | ((and (eq (car-safe x) 'sdev) | |
1653 | (Math-numberp b)) | |
1654 | (math-make-sdev (calcFunc-log (nth 1 x) b) | |
1655 | (math-div (nth 2 x) | |
1656 | (math-mul (nth 1 x) | |
1657 | (math-log-base-raw b))))) | |
1658 | ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) | |
1659 | (not (math-intv-constp x))) | |
1660 | (math-realp b)) | |
1661 | (math-make-intv (nth 1 x) | |
1662 | (calcFunc-log (nth 2 x) b) | |
1663 | (calcFunc-log (nth 3 x) b))) | |
1664 | ((or (eq (car-safe x) 'intv) (eq (car-safe b) 'intv)) | |
1665 | (math-div (calcFunc-ln x) (calcFunc-ln b))) | |
1666 | ((or (math-infinitep x) | |
1667 | (math-infinitep b)) | |
1668 | (math-div (calcFunc-ln x) (calcFunc-ln b))) | |
1669 | (t (if (Math-numberp b) | |
1670 | (calc-record-why 'numberp x) | |
1671 | (calc-record-why 'numberp b)) | |
491c3062 | 1672 | (list 'calcFunc-log x b)))) |
136211a9 EZ |
1673 | |
1674 | (defun calcFunc-alog (x &optional b) | |
1675 | (cond ((or (null b) (equal b '(var e var-e))) | |
1676 | (math-normalize (list 'calcFunc-exp x))) | |
491c3062 | 1677 | (t (math-pow b x)))) |
136211a9 EZ |
1678 | |
1679 | (defun calcFunc-ilog (x b) | |
1680 | (if (and (math-natnump x) (not (eq x 0)) | |
1681 | (math-natnump b) (not (eq b 0))) | |
1682 | (if (eq b 1) | |
1683 | (math-reject-arg x "*Logarithm base one") | |
1684 | (if (Math-natnum-lessp x b) | |
1685 | 0 | |
1686 | (cdr (math-integer-log x b)))) | |
491c3062 | 1687 | (math-floor (calcFunc-log x b)))) |
136211a9 EZ |
1688 | |
1689 | (defun math-integer-log (x b) | |
1690 | (let ((pows (list b)) | |
1691 | (pow (math-sqr b)) | |
1692 | next | |
1693 | sum n) | |
1694 | (while (not (Math-lessp x pow)) | |
1695 | (setq pows (cons pow pows) | |
1696 | pow (math-sqr pow))) | |
1697 | (setq n (lsh 1 (1- (length pows))) | |
1698 | sum n | |
1699 | pow (car pows)) | |
1700 | (while (and (setq pows (cdr pows)) | |
1701 | (Math-lessp pow x)) | |
1702 | (setq n (/ n 2) | |
1703 | next (math-mul pow (car pows))) | |
1704 | (or (Math-lessp x next) | |
1705 | (setq pow next | |
1706 | sum (+ sum n)))) | |
491c3062 | 1707 | (cons (equal pow x) sum))) |
136211a9 EZ |
1708 | |
1709 | ||
3132f345 | 1710 | (defvar math-log-base-cache nil) |
136211a9 EZ |
1711 | (defun math-log-base-raw (b) ; [N N] |
1712 | (if (not (and (equal (car math-log-base-cache) b) | |
1713 | (eq (nth 1 math-log-base-cache) calc-internal-prec))) | |
1714 | (setq math-log-base-cache (list b calc-internal-prec | |
1715 | (math-ln-raw (math-float b))))) | |
491c3062 | 1716 | (nth 2 math-log-base-cache)) |
136211a9 EZ |
1717 | |
1718 | (defun calcFunc-lnp1 (x) ; [N N] [Public] | |
1719 | (cond ((Math-equal-int x -1) | |
1720 | (if calc-infinite-mode | |
1721 | '(neg (var inf var-inf)) | |
1722 | (math-reject-arg x "*Logarithm of zero"))) | |
1723 | ((eq x 0) 0) | |
1724 | ((math-zerop x) '(float 0 0)) | |
1725 | (calc-symbolic-mode (signal 'inexact-result nil)) | |
1726 | ((Math-numberp x) | |
1727 | (math-with-extra-prec 2 | |
1728 | (let ((x (math-float x))) | |
1729 | (if (and (eq (car x) 'float) | |
1730 | (math-lessp-float x '(float 5 -1)) | |
1731 | (math-lessp-float '(float -5 -1) x)) | |
1732 | (math-ln-plus-1-raw x) | |
1733 | (math-ln-raw (math-add-float x '(float 1 0))))))) | |
1734 | ((eq (car-safe x) 'sdev) | |
1735 | (math-make-sdev (calcFunc-lnp1 (nth 1 x)) | |
1736 | (math-div (nth 2 x) (math-add (nth 1 x) 1)))) | |
1737 | ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) | |
1738 | (not (math-intv-constp x)))) | |
1739 | (math-make-intv (nth 1 x) | |
1740 | (calcFunc-lnp1 (nth 2 x)) | |
1741 | (calcFunc-lnp1 (nth 3 x)))) | |
1742 | ((math-infinitep x) | |
1743 | (if (equal x '(var nan var-nan)) | |
1744 | x | |
1745 | '(var inf var-inf))) | |
1746 | (t (calc-record-why 'numberp x) | |
491c3062 | 1747 | (list 'calcFunc-lnp1 x)))) |
136211a9 EZ |
1748 | |
1749 | (defun math-ln-raw (x) ; [N N] --- must be float format! | |
1750 | (cond ((eq (car-safe x) 'cplx) | |
1751 | (list 'cplx | |
1752 | (math-mul-float (math-ln-raw | |
1753 | (math-add-float (math-sqr-float (nth 1 x)) | |
1754 | (math-sqr-float (nth 2 x)))) | |
1755 | '(float 5 -1)) | |
1756 | (math-arctan2-raw (nth 2 x) (nth 1 x)))) | |
1757 | ((eq (car x) 'polar) | |
1758 | (math-polar (list 'cplx | |
1759 | (math-ln-raw (nth 1 x)) | |
1760 | (math-to-radians (nth 2 x))))) | |
1761 | ((Math-equal-int x 1) | |
1762 | '(float 0 0)) | |
1763 | (calc-symbolic-mode (signal 'inexact-result nil)) | |
1764 | ((math-posp (nth 1 x)) ; positive and real | |
17cc361e JB |
1765 | (cond |
1766 | ((math-use-emacs-fn 'log x)) | |
1767 | (t | |
1768 | (let ((xdigs (1- (math-numdigs (nth 1 x))))) | |
1769 | (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs))) | |
1770 | (math-mul-float (math-float (+ (nth 2 x) xdigs)) | |
1771 | (math-ln-10))))))) | |
136211a9 EZ |
1772 | ((math-zerop x) |
1773 | (math-reject-arg x "*Logarithm of zero")) | |
1774 | ((eq calc-complex-mode 'polar) ; negative and real | |
1775 | (math-polar | |
1776 | (list 'cplx ; negative and real | |
1777 | (math-ln-raw (math-neg-float x)) | |
1778 | (math-pi)))) | |
1779 | (t (list 'cplx ; negative and real | |
1780 | (math-ln-raw (math-neg-float x)) | |
491c3062 | 1781 | (math-pi))))) |
136211a9 EZ |
1782 | |
1783 | (defun math-ln-raw-2 (x) ; [F F] | |
1784 | (cond ((math-lessp-float '(float 14 -1) x) | |
1785 | (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1))) | |
1786 | (math-ln-2))) | |
1787 | (t ; now .7 < x <= 1.4 | |
1788 | (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0)) | |
491c3062 | 1789 | (math-add-float x '(float 1 0))))))) |
136211a9 EZ |
1790 | |
1791 | (defun math-ln-raw-3 (x) ; [F F] | |
1792 | (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x)) | |
491c3062 | 1793 | '(float 2 0))) |
136211a9 EZ |
1794 | |
1795 | ;;; Compute ln((1+x)/(1-x)) | |
1796 | (defun math-ln-raw-series (sum n x xsqr) | |
1797 | (math-working "log" sum) | |
1798 | (let* ((nextx (math-mul-float x xsqr)) | |
1799 | (nextsum (math-add-float sum (math-div-float nextx (math-float n))))) | |
1800 | (if (math-nearly-equal-float sum nextsum) | |
1801 | sum | |
491c3062 | 1802 | (math-ln-raw-series nextsum (+ n 2) nextx xsqr)))) |
136211a9 EZ |
1803 | |
1804 | (defun math-ln-plus-1-raw (x) | |
491c3062 | 1805 | (math-lnp1-series x 2 x (math-neg x))) |
136211a9 EZ |
1806 | |
1807 | (defun math-lnp1-series (sum n xpow x) | |
1808 | (math-working "lnp1" sum) | |
1809 | (let* ((nextx (math-mul-float xpow x)) | |
1810 | (nextsum (math-add-float sum (math-div-float nextx (math-float n))))) | |
1811 | (if (math-nearly-equal-float sum nextsum) | |
1812 | sum | |
491c3062 | 1813 | (math-lnp1-series nextsum (1+ n) nextx x)))) |
136211a9 | 1814 | |
4603f755 | 1815 | (defconst math-approx-ln-10 |
82472f45 | 1816 | (math-read-number-simple "2.302585092994045684018") |
94c95a35 | 1817 | "An approximation for ln(10).") |
4603f755 JB |
1818 | |
1819 | (math-defcache math-ln-10 math-approx-ln-10 | |
136211a9 EZ |
1820 | (math-ln-raw-2 '(float 1 1))) |
1821 | ||
4603f755 | 1822 | (defconst math-approx-ln-2 |
82472f45 | 1823 | (math-read-number-simple "0.693147180559945309417") |
94c95a35 | 1824 | "An approximation for ln(2).") |
4603f755 JB |
1825 | |
1826 | (math-defcache math-ln-2 math-approx-ln-2 | |
136211a9 EZ |
1827 | (math-ln-raw-3 (math-float '(frac 1 3)))) |
1828 | ||
1829 | ||
1830 | ||
1831 | ;;; Hyperbolic functions. | |
1832 | ||
1833 | (defun calcFunc-sinh (x) ; [N N] [Public] | |
1834 | (cond ((eq x 0) 0) | |
1835 | (math-expand-formulas | |
1836 | (math-normalize | |
1837 | (list '/ (list '- (list 'calcFunc-exp x) | |
1838 | (list 'calcFunc-exp (list 'neg x))) 2))) | |
1839 | ((Math-numberp x) | |
1840 | (if calc-symbolic-mode (signal 'inexact-result nil)) | |
1841 | (math-with-extra-prec 2 | |
1842 | (let ((expx (math-exp-raw (math-float x)))) | |
1843 | (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1))))) | |
1844 | ((eq (car-safe x) 'sdev) | |
1845 | (math-make-sdev (calcFunc-sinh (nth 1 x)) | |
1846 | (math-mul (nth 2 x) (calcFunc-cosh (nth 1 x))))) | |
1847 | ((eq (car x) 'intv) | |
1848 | (math-sort-intv (nth 1 x) | |
1849 | (calcFunc-sinh (nth 2 x)) | |
1850 | (calcFunc-sinh (nth 3 x)))) | |
1851 | ((or (equal x '(var inf var-inf)) | |
1852 | (equal x '(neg (var inf var-inf))) | |
1853 | (equal x '(var nan var-nan))) | |
1854 | x) | |
1855 | (t (calc-record-why 'numberp x) | |
491c3062 | 1856 | (list 'calcFunc-sinh x)))) |
136211a9 EZ |
1857 | (put 'calcFunc-sinh 'math-expandable t) |
1858 | ||
1859 | (defun calcFunc-cosh (x) ; [N N] [Public] | |
1860 | (cond ((eq x 0) 1) | |
1861 | (math-expand-formulas | |
1862 | (math-normalize | |
1863 | (list '/ (list '+ (list 'calcFunc-exp x) | |
1864 | (list 'calcFunc-exp (list 'neg x))) 2))) | |
1865 | ((Math-numberp x) | |
1866 | (if calc-symbolic-mode (signal 'inexact-result nil)) | |
1867 | (math-with-extra-prec 2 | |
1868 | (let ((expx (math-exp-raw (math-float x)))) | |
1869 | (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1))))) | |
1870 | ((eq (car-safe x) 'sdev) | |
1871 | (math-make-sdev (calcFunc-cosh (nth 1 x)) | |
1872 | (math-mul (nth 2 x) | |
1873 | (calcFunc-sinh (nth 1 x))))) | |
1874 | ((and (eq (car x) 'intv) (math-intv-constp x)) | |
1875 | (setq x (math-abs x)) | |
1876 | (math-sort-intv (nth 1 x) | |
1877 | (calcFunc-cosh (nth 2 x)) | |
1878 | (calcFunc-cosh (nth 3 x)))) | |
1879 | ((or (equal x '(var inf var-inf)) | |
1880 | (equal x '(neg (var inf var-inf))) | |
1881 | (equal x '(var nan var-nan))) | |
1882 | (math-abs x)) | |
1883 | (t (calc-record-why 'numberp x) | |
491c3062 | 1884 | (list 'calcFunc-cosh x)))) |
136211a9 EZ |
1885 | (put 'calcFunc-cosh 'math-expandable t) |
1886 | ||
1887 | (defun calcFunc-tanh (x) ; [N N] [Public] | |
1888 | (cond ((eq x 0) 0) | |
1889 | (math-expand-formulas | |
1890 | (math-normalize | |
1891 | (let ((expx (list 'calcFunc-exp x)) | |
1892 | (expmx (list 'calcFunc-exp (list 'neg x)))) | |
1893 | (math-normalize | |
1894 | (list '/ (list '- expx expmx) (list '+ expx expmx)))))) | |
1895 | ((Math-numberp x) | |
1896 | (if calc-symbolic-mode (signal 'inexact-result nil)) | |
1897 | (math-with-extra-prec 2 | |
1898 | (let* ((expx (calcFunc-exp (math-float x))) | |
1899 | (expmx (math-div 1 expx))) | |
1900 | (math-div (math-sub expx expmx) | |
1901 | (math-add expx expmx))))) | |
1902 | ((eq (car-safe x) 'sdev) | |
1903 | (math-make-sdev (calcFunc-tanh (nth 1 x)) | |
1904 | (math-div (nth 2 x) | |
1905 | (math-sqr (calcFunc-cosh (nth 1 x)))))) | |
1906 | ((eq (car x) 'intv) | |
1907 | (math-sort-intv (nth 1 x) | |
1908 | (calcFunc-tanh (nth 2 x)) | |
1909 | (calcFunc-tanh (nth 3 x)))) | |
1910 | ((equal x '(var inf var-inf)) | |
1911 | 1) | |
1912 | ((equal x '(neg (var inf var-inf))) | |
1913 | -1) | |
1914 | ((equal x '(var nan var-nan)) | |
1915 | x) | |
1916 | (t (calc-record-why 'numberp x) | |
491c3062 | 1917 | (list 'calcFunc-tanh x)))) |
136211a9 EZ |
1918 | (put 'calcFunc-tanh 'math-expandable t) |
1919 | ||
f53e6c20 JB |
1920 | (defun calcFunc-sech (x) ; [N N] [Public] |
1921 | (cond ((eq x 0) 1) | |
1922 | (math-expand-formulas | |
1923 | (math-normalize | |
1924 | (list '/ 2 (list '+ (list 'calcFunc-exp x) | |
1925 | (list 'calcFunc-exp (list 'neg x)))))) | |
1926 | ((Math-numberp x) | |
1927 | (if calc-symbolic-mode (signal 'inexact-result nil)) | |
1928 | (math-with-extra-prec 2 | |
1929 | (let ((expx (math-exp-raw (math-float x)))) | |
1930 | (math-div '(float 2 0) (math-add expx (math-div 1 expx)))))) | |
1931 | ((eq (car-safe x) 'sdev) | |
1932 | (math-make-sdev (calcFunc-sech (nth 1 x)) | |
1933 | (math-mul (nth 2 x) | |
1934 | (math-mul (calcFunc-sech (nth 1 x)) | |
1935 | (calcFunc-tanh (nth 1 x)))))) | |
1936 | ((and (eq (car x) 'intv) (math-intv-constp x)) | |
1937 | (setq x (math-abs x)) | |
1938 | (math-sort-intv (nth 1 x) | |
1939 | (calcFunc-sech (nth 2 x)) | |
1940 | (calcFunc-sech (nth 3 x)))) | |
1941 | ((or (equal x '(var inf var-inf)) | |
1942 | (equal x '(neg (var inf var-inf)))) | |
1943 | 0) | |
1944 | ((equal x '(var nan var-nan)) | |
1945 | x) | |
1946 | (t (calc-record-why 'numberp x) | |
1947 | (list 'calcFunc-sech x)))) | |
1948 | (put 'calcFunc-sech 'math-expandable t) | |
1949 | ||
1950 | (defun calcFunc-csch (x) ; [N N] [Public] | |
1951 | (cond ((eq x 0) (math-div 1 0)) | |
1952 | (math-expand-formulas | |
1953 | (math-normalize | |
1954 | (list '/ 2 (list '- (list 'calcFunc-exp x) | |
1955 | (list 'calcFunc-exp (list 'neg x)))))) | |
1956 | ((Math-numberp x) | |
1957 | (if calc-symbolic-mode (signal 'inexact-result nil)) | |
1958 | (math-with-extra-prec 2 | |
1959 | (let ((expx (math-exp-raw (math-float x)))) | |
1960 | (math-div '(float 2 0) (math-add expx (math-div -1 expx)))))) | |
1961 | ((eq (car-safe x) 'sdev) | |
1962 | (math-make-sdev (calcFunc-csch (nth 1 x)) | |
1963 | (math-mul (nth 2 x) | |
1964 | (math-mul (calcFunc-csch (nth 1 x)) | |
1965 | (calcFunc-coth (nth 1 x)))))) | |
1966 | ((eq (car x) 'intv) | |
1967 | (if (and (Math-negp (nth 2 x)) | |
1968 | (Math-posp (nth 3 x))) | |
1969 | '(intv 3 (neg (var inf var-inf)) (var inf var-inf)) | |
1970 | (math-sort-intv (nth 1 x) | |
1971 | (calcFunc-csch (nth 2 x)) | |
1972 | (calcFunc-csch (nth 3 x))))) | |
1973 | ((or (equal x '(var inf var-inf)) | |
1974 | (equal x '(neg (var inf var-inf)))) | |
1975 | 0) | |
1976 | ((equal x '(var nan var-nan)) | |
1977 | x) | |
1978 | (t (calc-record-why 'numberp x) | |
1979 | (list 'calcFunc-csch x)))) | |
1980 | (put 'calcFunc-csch 'math-expandable t) | |
1981 | ||
1982 | (defun calcFunc-coth (x) ; [N N] [Public] | |
1983 | (cond ((eq x 0) (math-div 1 0)) | |
1984 | (math-expand-formulas | |
1985 | (math-normalize | |
1986 | (let ((expx (list 'calcFunc-exp x)) | |
1987 | (expmx (list 'calcFunc-exp (list 'neg x)))) | |
1988 | (math-normalize | |
1989 | (list '/ (list '+ expx expmx) (list '- expx expmx)))))) | |
1990 | ((Math-numberp x) | |
1991 | (if calc-symbolic-mode (signal 'inexact-result nil)) | |
1992 | (math-with-extra-prec 2 | |
1993 | (let* ((expx (calcFunc-exp (math-float x))) | |
1994 | (expmx (math-div 1 expx))) | |
1995 | (math-div (math-add expx expmx) | |
1996 | (math-sub expx expmx))))) | |
1997 | ((eq (car-safe x) 'sdev) | |
1998 | (math-make-sdev (calcFunc-coth (nth 1 x)) | |
1999 | (math-div (nth 2 x) | |
2000 | (math-sqr (calcFunc-sinh (nth 1 x)))))) | |
2001 | ((eq (car x) 'intv) | |
2002 | (if (and (Math-negp (nth 2 x)) | |
2003 | (Math-posp (nth 3 x))) | |
2004 | '(intv 3 (neg (var inf var-inf)) (var inf var-inf)) | |
2005 | (math-sort-intv (nth 1 x) | |
2006 | (calcFunc-coth (nth 2 x)) | |
2007 | (calcFunc-coth (nth 3 x))))) | |
2008 | ((equal x '(var inf var-inf)) | |
2009 | 1) | |
2010 | ((equal x '(neg (var inf var-inf))) | |
2011 | -1) | |
2012 | ((equal x '(var nan var-nan)) | |
2013 | x) | |
2014 | (t (calc-record-why 'numberp x) | |
2015 | (list 'calcFunc-coth x)))) | |
2016 | (put 'calcFunc-coth 'math-expandable t) | |
2017 | ||
136211a9 EZ |
2018 | (defun calcFunc-arcsinh (x) ; [N N] [Public] |
2019 | (cond ((eq x 0) 0) | |
2020 | (math-expand-formulas | |
2021 | (math-normalize | |
2022 | (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt | |
2023 | (list '+ (list '^ x 2) 1)))))) | |
2024 | ((Math-numberp x) | |
2025 | (if calc-symbolic-mode (signal 'inexact-result nil)) | |
2026 | (math-with-extra-prec 2 | |
2027 | (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x) | |
2028 | '(float 1 0))))))) | |
2029 | ((eq (car-safe x) 'sdev) | |
2030 | (math-make-sdev (calcFunc-arcsinh (nth 1 x)) | |
2031 | (math-div (nth 2 x) | |
2032 | (math-sqrt | |
2033 | (math-add (math-sqr (nth 1 x)) 1))))) | |
2034 | ((eq (car x) 'intv) | |
2035 | (math-sort-intv (nth 1 x) | |
2036 | (calcFunc-arcsinh (nth 2 x)) | |
2037 | (calcFunc-arcsinh (nth 3 x)))) | |
2038 | ((or (equal x '(var inf var-inf)) | |
2039 | (equal x '(neg (var inf var-inf))) | |
2040 | (equal x '(var nan var-nan))) | |
2041 | x) | |
2042 | (t (calc-record-why 'numberp x) | |
491c3062 | 2043 | (list 'calcFunc-arcsinh x)))) |
136211a9 EZ |
2044 | (put 'calcFunc-arcsinh 'math-expandable t) |
2045 | ||
2046 | (defun calcFunc-arccosh (x) ; [N N] [Public] | |
2047 | (cond ((eq x 1) 0) | |
2048 | ((and (eq x -1) calc-symbolic-mode) | |
2049 | '(var pi var-pi)) | |
2050 | ((and (eq x 0) calc-symbolic-mode) | |
2051 | (math-div (math-mul '(var pi var-pi) '(var i var-i)) 2)) | |
2052 | (math-expand-formulas | |
2053 | (math-normalize | |
2054 | (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt | |
2055 | (list '- (list '^ x 2) 1)))))) | |
2056 | ((Math-numberp x) | |
2057 | (if calc-symbolic-mode (signal 'inexact-result nil)) | |
2058 | (if (Math-equal-int x -1) | |
2059 | (math-imaginary (math-pi)) | |
2060 | (math-with-extra-prec 2 | |
2061 | (if (or t ; need to do this even in the real case! | |
2062 | (memq (car-safe x) '(cplx polar))) | |
2063 | (let ((xp1 (math-add 1 x))) ; this gets the branch cuts right | |
2064 | (math-ln-raw | |
2065 | (math-add x (math-mul xp1 | |
2066 | (math-sqrt-raw | |
2067 | (math-div (math-sub | |
2068 | x | |
2069 | '(float 1 0)) | |
2070 | xp1)))))) | |
2071 | (math-ln-raw | |
2072 | (math-add x (math-sqrt-raw (math-add (math-sqr x) | |
2073 | '(float -1 0))))))))) | |
2074 | ((eq (car-safe x) 'sdev) | |
2075 | (math-make-sdev (calcFunc-arccosh (nth 1 x)) | |
2076 | (math-div (nth 2 x) | |
2077 | (math-sqrt | |
2078 | (math-add (math-sqr (nth 1 x)) -1))))) | |
2079 | ((eq (car x) 'intv) | |
2080 | (math-sort-intv (nth 1 x) | |
2081 | (calcFunc-arccosh (nth 2 x)) | |
2082 | (calcFunc-arccosh (nth 3 x)))) | |
2083 | ((or (equal x '(var inf var-inf)) | |
2084 | (equal x '(neg (var inf var-inf))) | |
2085 | (equal x '(var nan var-nan))) | |
2086 | x) | |
2087 | (t (calc-record-why 'numberp x) | |
491c3062 | 2088 | (list 'calcFunc-arccosh x)))) |
136211a9 EZ |
2089 | (put 'calcFunc-arccosh 'math-expandable t) |
2090 | ||
2091 | (defun calcFunc-arctanh (x) ; [N N] [Public] | |
2092 | (cond ((eq x 0) 0) | |
2093 | ((and (Math-equal-int x 1) calc-infinite-mode) | |
2094 | '(var inf var-inf)) | |
2095 | ((and (Math-equal-int x -1) calc-infinite-mode) | |
2096 | '(neg (var inf var-inf))) | |
2097 | (math-expand-formulas | |
2098 | (list '/ (list '- | |
2099 | (list 'calcFunc-ln (list '+ 1 x)) | |
2100 | (list 'calcFunc-ln (list '- 1 x))) 2)) | |
2101 | ((Math-numberp x) | |
2102 | (if calc-symbolic-mode (signal 'inexact-result nil)) | |
2103 | (math-with-extra-prec 2 | |
2104 | (if (or (memq (car-safe x) '(cplx polar)) | |
2105 | (Math-lessp 1 x)) | |
2106 | (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x)) | |
2107 | (math-ln-raw (math-sub '(float 1 0) x))) | |
2108 | '(float 5 -1)) | |
2109 | (if (and (math-equal-int x 1) calc-infinite-mode) | |
2110 | '(var inf var-inf) | |
2111 | (if (and (math-equal-int x -1) calc-infinite-mode) | |
2112 | '(neg (var inf var-inf)) | |
2113 | (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x) | |
2114 | (math-sub 1 x))) | |
2115 | '(float 5 -1))))))) | |
2116 | ((eq (car-safe x) 'sdev) | |
2117 | (math-make-sdev (calcFunc-arctanh (nth 1 x)) | |
2118 | (math-div (nth 2 x) | |
2119 | (math-sub 1 (math-sqr (nth 1 x)))))) | |
2120 | ((eq (car x) 'intv) | |
2121 | (math-sort-intv (nth 1 x) | |
2122 | (calcFunc-arctanh (nth 2 x)) | |
2123 | (calcFunc-arctanh (nth 3 x)))) | |
2124 | ((equal x '(var nan var-nan)) | |
2125 | x) | |
2126 | (t (calc-record-why 'numberp x) | |
491c3062 | 2127 | (list 'calcFunc-arctanh x)))) |
136211a9 EZ |
2128 | (put 'calcFunc-arctanh 'math-expandable t) |
2129 | ||
2130 | ||
2131 | ;;; Convert A from HMS or degrees to radians. | |
2132 | (defun calcFunc-rad (a) ; [R R] [Public] | |
2133 | (cond ((or (Math-numberp a) | |
2134 | (eq (car a) 'intv)) | |
2135 | (math-with-extra-prec 2 | |
2136 | (math-mul a (math-pi-over-180)))) | |
2137 | ((eq (car a) 'hms) | |
2138 | (math-from-hms a 'rad)) | |
2139 | ((eq (car a) 'sdev) | |
2140 | (math-make-sdev (calcFunc-rad (nth 1 a)) | |
2141 | (calcFunc-rad (nth 2 a)))) | |
2142 | (math-expand-formulas | |
2143 | (math-div (math-mul a '(var pi var-pi)) 180)) | |
2144 | ((math-infinitep a) a) | |
491c3062 | 2145 | (t (list 'calcFunc-rad a)))) |
136211a9 EZ |
2146 | (put 'calcFunc-rad 'math-expandable t) |
2147 | ||
2148 | ;;; Convert A from HMS or radians to degrees. | |
2149 | (defun calcFunc-deg (a) ; [R R] [Public] | |
2150 | (cond ((or (Math-numberp a) | |
2151 | (eq (car a) 'intv)) | |
2152 | (math-with-extra-prec 2 | |
2153 | (math-div a (math-pi-over-180)))) | |
2154 | ((eq (car a) 'hms) | |
2155 | (math-from-hms a 'deg)) | |
2156 | ((eq (car a) 'sdev) | |
2157 | (math-make-sdev (calcFunc-deg (nth 1 a)) | |
2158 | (calcFunc-deg (nth 2 a)))) | |
2159 | (math-expand-formulas | |
2160 | (math-div (math-mul 180 a) '(var pi var-pi))) | |
2161 | ((math-infinitep a) a) | |
491c3062 | 2162 | (t (list 'calcFunc-deg a)))) |
136211a9 EZ |
2163 | (put 'calcFunc-deg 'math-expandable t) |
2164 | ||
95995a85 JB |
2165 | (provide 'calc-math) |
2166 | ||
cbee283d | 2167 | ;; arch-tag: c7367e8e-d0b8-4f70-8577-2fb3f31dbb4c |
491c3062 | 2168 | ;;; calc-math.el ends here |