1 ;;; calc-math.el --- mathematical functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2005 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
6 ;; Maintainer: Jay Belanger <belanger@truman.edu>
8 ;; This file is part of GNU Emacs.
10 ;; GNU Emacs is distributed in the hope that it will be useful,
11 ;; but WITHOUT ANY WARRANTY. No author or distributor
12 ;; accepts responsibility to anyone for the consequences of using it
13 ;; or for whether it serves any particular purpose or works at all,
14 ;; unless he says so in writing. Refer to the GNU Emacs General Public
15 ;; License for full details.
17 ;; Everyone is granted permission to copy, modify and redistribute
18 ;; GNU Emacs, but only under the conditions described in the
19 ;; GNU Emacs General Public License. A copy of this license is
20 ;; supposed to have been given to you along with GNU Emacs so you
21 ;; can know your rights and responsibilities. It should be in a
22 ;; file named COPYING. Among other things, the copyright notice
23 ;; and this notice must be preserved on all copies.
29 ;; This file is autoloaded from calc-ext.el.
34 (defun calc-sqrt (arg)
38 (calc-unary-op "^2" 'calcFunc-sqr arg
)
39 (calc-unary-op "sqrt" 'calcFunc-sqrt arg
))))
41 (defun calc-isqrt (arg)
45 (calc-unary-op "^2" 'calcFunc-sqr arg
)
46 (calc-unary-op "isqt" 'calcFunc-isqrt arg
))))
49 (defun calc-hypot (arg)
52 (calc-binary-op "hypt" 'calcFunc-hypot arg
)))
59 (defun calc-log10 (arg)
61 (calc-hyperbolic-func)
68 (calc-binary-op "alog" 'calcFunc-alog arg
)
69 (calc-binary-op "log" 'calcFunc-log arg
))))
71 (defun calc-ilog (arg)
75 (calc-binary-op "alog" 'calcFunc-alog arg
)
76 (calc-binary-op "ilog" 'calcFunc-ilog arg
))))
78 (defun calc-lnp1 (arg)
86 (if (calc-is-hyperbolic)
88 (calc-unary-op "lg10" 'calcFunc-log10 arg
)
89 (calc-unary-op "10^" 'calcFunc-exp10 arg
))
91 (calc-unary-op "ln" 'calcFunc-ln arg
)
92 (calc-unary-op "exp" 'calcFunc-exp arg
)))))
94 (defun calc-expm1 (arg)
98 (calc-unary-op "ln+1" 'calcFunc-lnp1 arg
)
99 (calc-unary-op "ex-1" 'calcFunc-expm1 arg
))))
104 (if (calc-is-inverse)
105 (if (calc-is-hyperbolic)
106 (if calc-symbolic-mode
107 (calc-pop-push-record 0 "phi" '(var phi var-phi
))
108 (calc-pop-push-record 0 "phi" (math-phi)))
109 (if calc-symbolic-mode
110 (calc-pop-push-record 0 "gmma" '(var gamma var-gamma
))
111 (calc-pop-push-record 0 "gmma" (math-gamma-const))))
112 (if (calc-is-hyperbolic)
113 (if calc-symbolic-mode
114 (calc-pop-push-record 0 "e" '(var e var-e
))
115 (calc-pop-push-record 0 "e" (math-e)))
116 (if calc-symbolic-mode
117 (calc-pop-push-record 0 "pi" '(var pi var-pi
))
118 (calc-pop-push-record 0 "pi" (math-pi)))))))
120 (defun calc-sin (arg)
123 (if (calc-is-hyperbolic)
124 (if (calc-is-inverse)
125 (calc-unary-op "asnh" 'calcFunc-arcsinh arg
)
126 (calc-unary-op "sinh" 'calcFunc-sinh arg
))
127 (if (calc-is-inverse)
128 (calc-unary-op "asin" 'calcFunc-arcsin arg
)
129 (calc-unary-op "sin" 'calcFunc-sin arg
)))))
131 (defun calc-arcsin (arg)
136 (defun calc-sinh (arg)
138 (calc-hyperbolic-func)
141 (defun calc-arcsinh (arg)
144 (calc-hyperbolic-func)
147 (defun calc-sec (arg)
150 (if (calc-is-hyperbolic)
151 (calc-unary-op "sech" 'calcFunc-sech arg
)
152 (calc-unary-op "sec" 'calcFunc-sec arg
))))
154 (defun calc-sech (arg)
156 (calc-hyperbolic-func)
159 (defun calc-cos (arg)
162 (if (calc-is-hyperbolic)
163 (if (calc-is-inverse)
164 (calc-unary-op "acsh" 'calcFunc-arccosh arg
)
165 (calc-unary-op "cosh" 'calcFunc-cosh arg
))
166 (if (calc-is-inverse)
167 (calc-unary-op "acos" 'calcFunc-arccos arg
)
168 (calc-unary-op "cos" 'calcFunc-cos arg
)))))
170 (defun calc-arccos (arg)
175 (defun calc-cosh (arg)
177 (calc-hyperbolic-func)
180 (defun calc-arccosh (arg)
183 (calc-hyperbolic-func)
186 (defun calc-csc (arg)
189 (if (calc-is-hyperbolic)
190 (calc-unary-op "csch" 'calcFunc-csch arg
)
191 (calc-unary-op "csc" 'calcFunc-csc arg
))))
193 (defun calc-csch (arg)
195 (calc-hyperbolic-func)
198 (defun calc-sincos ()
201 (if (calc-is-inverse)
202 (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos
(calc-top-n 1)))
203 (calc-enter-result 1 "sncs" (list 'calcFunc-sincos
(calc-top-n 1))))))
205 (defun calc-tan (arg)
208 (if (calc-is-hyperbolic)
209 (if (calc-is-inverse)
210 (calc-unary-op "atnh" 'calcFunc-arctanh arg
)
211 (calc-unary-op "tanh" 'calcFunc-tanh arg
))
212 (if (calc-is-inverse)
213 (calc-unary-op "atan" 'calcFunc-arctan arg
)
214 (calc-unary-op "tan" 'calcFunc-tan arg
)))))
216 (defun calc-arctan (arg)
221 (defun calc-tanh (arg)
223 (calc-hyperbolic-func)
226 (defun calc-arctanh (arg)
229 (calc-hyperbolic-func)
232 (defun calc-cot (arg)
235 (if (calc-is-hyperbolic)
236 (calc-unary-op "coth" 'calcFunc-coth arg
)
237 (calc-unary-op "cot" 'calcFunc-cot arg
))))
239 (defun calc-coth (arg)
241 (calc-hyperbolic-func)
244 (defun calc-arctan2 ()
247 (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2
(calc-top-list-n 2)))))
249 (defun calc-conj (arg)
252 (calc-unary-op "conj" 'calcFunc-conj arg
)))
254 (defun calc-imaginary ()
257 (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))))
259 (defun calc-to-degrees (arg)
262 (calc-unary-op ">deg" 'calcFunc-deg arg
)))
264 (defun calc-to-radians (arg)
267 (calc-unary-op ">rad" 'calcFunc-rad arg
)))
270 (defun calc-degrees-mode (arg)
274 (calc-change-mode 'calc-angle-mode
'deg
)
275 (message "Angles measured in degrees")))
276 ((= arg
2) (calc-radians-mode))
277 ((= arg
3) (calc-hms-mode))
278 (t (error "Prefix argument out of range"))))
280 (defun calc-radians-mode ()
283 (calc-change-mode 'calc-angle-mode
'rad
)
284 (message "Angles measured in radians")))
287 ;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public]
288 ;;; This method takes advantage of the fact that Newton's method starting
289 ;;; with an overestimate always works, even using truncating integer division!
290 (defun math-isqrt (a)
291 (cond ((Math-zerop a
) a
)
292 ((not (math-natnump a
))
293 (math-reject-arg a
'natnump
))
295 (math-isqrt-small a
))
297 (math-normalize (cons 'bigpos
(cdr (math-isqrt-bignum (cdr a
))))))))
299 (defun calcFunc-isqrt (a)
301 (math-isqrt (math-floor a
))
302 (math-floor (math-sqrt a
))))
305 ;;; This returns (flag . result) where the flag is t if A is a perfect square.
306 (defun math-isqrt-bignum (a) ; [P.l L]
307 (let ((len (length a
)))
309 (let* ((top (nthcdr (- len
2) a
)))
310 (math-isqrt-bignum-iter
314 (1+ (math-isqrt-small
315 (+ (* (nth 1 top
) 1000) (car top
)))))
317 (let* ((top (nth (1- len
) a
)))
318 (math-isqrt-bignum-iter
321 (list (1+ (math-isqrt-small top
)))
324 (defun math-isqrt-bignum-iter (a guess
) ; [l L l]
325 (math-working "isqrt" (cons 'bigpos guess
))
326 (let* ((q (math-div-bignum a guess
))
327 (s (math-add-bignum (car q
) guess
))
328 (g2 (math-div2-bignum s
))
329 (comp (math-compare-bignum g2 guess
)))
331 (math-isqrt-bignum-iter a g2
)
332 (cons (and (= comp
0)
333 (math-zerop-bignum (cdr q
))
337 (defun math-zerop-bignum (a)
340 (while (eq (car (setq a
(cdr a
))) 0))
343 (defun math-scale-bignum-3 (a n
) ; [L L S]
349 (defun math-isqrt-small (a) ; A > 0. [S S]
350 (let ((g (cond ((>= a
10000) 1000)
354 (while (< (setq g2
(/ (+ g
(/ a g
)) 2)) g
)
361 ;;; Compute the square root of a number.
362 ;;; [T N] if possible, else [F N] if possible, else [C N]. [Public]
365 (and (Math-zerop a
) a
)
366 (and (math-known-nonposp a
)
367 (math-imaginary (math-sqrt (math-neg a
))))
369 (let ((sqrt (math-isqrt-small a
)))
370 (if (= (* sqrt sqrt
) a
)
372 (if calc-symbolic-mode
373 (list 'calcFunc-sqrt a
)
374 (math-sqrt-float (math-float a
) (math-float sqrt
))))))
375 (and (eq (car-safe a
) 'bigpos
)
376 (let* ((res (math-isqrt-bignum (cdr a
)))
377 (sqrt (math-normalize (cons 'bigpos
(cdr res
)))))
380 (if calc-symbolic-mode
381 (list 'calcFunc-sqrt a
)
382 (math-sqrt-float (math-float a
) (math-float sqrt
))))))
383 (and (eq (car-safe a
) 'frac
)
384 (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a
)))))
385 (num-sqrt (math-normalize (cons 'bigpos
(cdr num-res
))))
386 (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a
)))))
387 (den-sqrt (math-normalize (cons 'bigpos
(cdr den-res
)))))
388 (if (and (car num-res
) (car den-res
))
389 (list 'frac num-sqrt den-sqrt
)
390 (if calc-symbolic-mode
391 (if (or (car num-res
) (car den-res
))
392 (math-div (if (car num-res
)
393 num-sqrt
(list 'calcFunc-sqrt
(nth 1 a
)))
395 den-sqrt
(list 'calcFunc-sqrt
(nth 2 a
))))
396 (list 'calcFunc-sqrt a
))
397 (math-sqrt-float (math-float a
)
398 (math-div (math-float num-sqrt
) den-sqrt
))))))
399 (and (eq (car-safe a
) 'float
)
400 (if calc-symbolic-mode
401 (if (= (%
(nth 2 a
) 2) 0)
402 (let ((res (math-isqrt-bignum
403 (cdr (Math-bignum-test (nth 1 a
))))))
405 (math-make-float (math-normalize
406 (cons 'bigpos
(cdr res
)))
408 (signal 'inexact-result nil
)))
409 (signal 'inexact-result nil
))
410 (math-sqrt-float a
)))
411 (and (eq (car-safe a
) 'cplx
)
412 (math-with-extra-prec 2
413 (let* ((d (math-abs a
))
414 (imag (math-sqrt (math-mul (math-sub d
(nth 1 a
))
417 (math-sqrt (math-mul (math-add d
(nth 1 a
)) '(float 5 -
1)))
418 (if (math-negp (nth 2 a
)) (math-neg imag
) imag
)))))
419 (and (eq (car-safe a
) 'polar
)
421 (math-sqrt (nth 1 a
))
422 (math-mul (nth 2 a
) '(float 5 -
1))))
423 (and (eq (car-safe a
) 'sdev
)
424 (let ((sqrt (math-sqrt (nth 1 a
))))
426 (math-div (nth 2 a
) (math-mul sqrt
2)))))
427 (and (eq (car-safe a
) 'intv
)
428 (not (math-negp (nth 2 a
)))
429 (math-make-intv (nth 1 a
) (math-sqrt (nth 2 a
)) (math-sqrt (nth 3 a
))))
430 (and (eq (car-safe a
) '*)
431 (or (math-known-nonnegp (nth 1 a
))
432 (math-known-nonnegp (nth 2 a
)))
433 (math-mul (math-sqrt (nth 1 a
)) (math-sqrt (nth 2 a
))))
434 (and (eq (car-safe a
) '/)
435 (or (and (math-known-nonnegp (nth 2 a
))
436 (math-div (math-sqrt (nth 1 a
)) (math-sqrt (nth 2 a
))))
437 (and (math-known-nonnegp (nth 1 a
))
438 (not (math-equal-int (nth 1 a
) 1))
439 (math-mul (math-sqrt (nth 1 a
))
440 (math-sqrt (math-div 1 (nth 2 a
)))))))
441 (and (eq (car-safe a
) '^
)
442 (math-known-evenp (nth 2 a
))
443 (math-known-realp (nth 1 a
))
444 (math-abs (math-pow (nth 1 a
) (math-div (nth 2 a
) 2))))
445 (let ((inf (math-infinitep a
)))
447 (math-mul (math-sqrt (math-infinite-dir a inf
)) inf
)))
449 (calc-record-why 'numberp a
)
450 (list 'calcFunc-sqrt a
))))
451 (defalias 'calcFunc-sqrt
'math-sqrt
)
453 (defun math-infinite-dir (a &optional inf
)
454 (or inf
(setq inf
(math-infinitep a
)))
455 (math-normalize (math-expr-subst a inf
1)))
457 (defun math-sqrt-float (a &optional guess
) ; [F F F]
458 (if calc-symbolic-mode
459 (signal 'inexact-result nil
)
460 (math-with-extra-prec 1 (math-sqrt-raw a guess
))))
462 (defun math-sqrt-raw (a &optional guess
) ; [F F F]
463 (if (not (Math-posp a
))
466 (let ((ldiff (- (math-numdigs (nth 1 a
)) 6)))
467 (or (= (%
(+ (nth 2 a
) ldiff
) 2) 0) (setq ldiff
(1+ ldiff
)))
468 (setq guess
(math-make-float (math-isqrt-small
469 (math-scale-int (nth 1 a
) (- ldiff
)))
470 (/ (+ (nth 2 a
) ldiff
) 2)))))
471 (math-sqrt-float-iter a guess
)))
473 (defun math-sqrt-float-iter (a guess
) ; [F F F]
474 (math-working "sqrt" guess
)
475 (let ((g2 (math-mul-float (math-add-float guess
(math-div-float a guess
))
477 (if (math-nearly-equal-float g2 guess
)
479 (math-sqrt-float-iter a g2
))))
481 ;;; True if A and B differ only in the last digit of precision. [P F F]
482 (defun math-nearly-equal-float (a b
)
483 (let ((ediff (- (nth 2 a
) (nth 2 b
))))
484 (cond ((= ediff
0) ;; Expanded out for speed
485 (setq ediff
(math-add (Math-integer-neg (nth 1 a
)) (nth 1 b
)))
487 (and (not (consp ediff
))
490 (= (math-numdigs (nth 1 a
)) calc-internal-prec
))))
492 (setq ediff
(math-add (Math-integer-neg (nth 1 b
))
493 (math-scale-int (nth 1 a
) 1)))
494 (and (not (consp ediff
))
497 (= (math-numdigs (nth 1 b
)) calc-internal-prec
)))
499 (setq ediff
(math-add (Math-integer-neg (nth 1 a
))
500 (math-scale-int (nth 1 b
) 1)))
501 (and (not (consp ediff
))
504 (= (math-numdigs (nth 1 a
)) calc-internal-prec
))))))
506 (defun math-nearly-equal (a b
) ; [P N N] [Public]
507 (setq a
(math-float a
))
508 (setq b
(math-float b
))
509 (if (eq (car a
) 'polar
) (setq a
(math-complex a
)))
510 (if (eq (car b
) 'polar
) (setq b
(math-complex b
)))
511 (if (eq (car a
) 'cplx
)
512 (if (eq (car b
) 'cplx
)
513 (and (or (math-nearly-equal-float (nth 1 a
) (nth 1 b
))
514 (and (math-nearly-zerop-float (nth 1 a
) (nth 2 a
))
515 (math-nearly-zerop-float (nth 1 b
) (nth 2 b
))))
516 (or (math-nearly-equal-float (nth 2 a
) (nth 2 b
))
517 (and (math-nearly-zerop-float (nth 2 a
) (nth 1 a
))
518 (math-nearly-zerop-float (nth 2 b
) (nth 1 b
)))))
519 (and (math-nearly-equal-float (nth 1 a
) b
)
520 (math-nearly-zerop-float (nth 2 a
) b
)))
521 (if (eq (car b
) 'cplx
)
522 (and (math-nearly-equal-float a
(nth 1 b
))
523 (math-nearly-zerop-float a
(nth 2 b
)))
524 (math-nearly-equal-float a b
))))
526 ;;; True if A is nearly zero compared to B. [P F F]
527 (defun math-nearly-zerop-float (a b
)
529 (<= (+ (math-numdigs (nth 1 a
)) (nth 2 a
))
530 (1+ (- (+ (math-numdigs (nth 1 b
)) (nth 2 b
)) calc-internal-prec
)))))
532 (defun math-nearly-zerop (a b
) ; [P N R] [Public]
533 (setq a
(math-float a
))
534 (setq b
(math-float b
))
535 (if (eq (car a
) 'cplx
)
536 (and (math-nearly-zerop-float (nth 1 a
) b
)
537 (math-nearly-zerop-float (nth 2 a
) b
))
538 (if (eq (car a
) 'polar
)
539 (math-nearly-zerop-float (nth 1 a
) b
)
540 (math-nearly-zerop-float a b
))))
542 ;;; This implementation could be improved, accuracy-wise.
543 (defun math-hypot (a b
)
544 (cond ((Math-zerop a
) (math-abs b
))
545 ((Math-zerop b
) (math-abs a
))
546 ((not (Math-scalarp a
))
547 (if (math-infinitep a
)
548 (if (math-infinitep b
)
553 (calc-record-why 'scalarp a
)
554 (list 'calcFunc-hypot a b
)))
555 ((not (Math-scalarp b
))
556 (if (math-infinitep b
)
558 (calc-record-why 'scalarp b
)
559 (list 'calcFunc-hypot a b
)))
560 ((and (Math-numberp a
) (Math-numberp b
))
561 (math-with-extra-prec 1
562 (math-sqrt (math-add (calcFunc-abssqr a
) (calcFunc-abssqr b
)))))
563 ((eq (car-safe a
) 'hms
)
564 (if (eq (car-safe b
) 'hms
) ; this helps sdev's of hms forms
565 (math-to-hms (math-hypot (math-from-hms a
'deg
)
566 (math-from-hms b
'deg
)))
567 (math-to-hms (math-hypot (math-from-hms a
'deg
) b
))))
568 ((eq (car-safe b
) 'hms
)
569 (math-to-hms (math-hypot a
(math-from-hms b
'deg
))))
571 (defalias 'calcFunc-hypot
'math-hypot
)
573 (defun calcFunc-sqr (x)
578 (defun math-nth-root (a n
)
579 (cond ((= n
2) (math-sqrt a
))
583 (let ((root (math-nth-root-integer a n
)))
586 (and (not calc-symbolic-mode
)
587 (math-nth-root-float (math-float a
) n
588 (math-float (cdr root
)))))))
589 ((eq (car-safe a
) 'frac
)
590 (let* ((num-root (math-nth-root-integer (nth 1 a
) n
))
591 (den-root (math-nth-root-integer (nth 2 a
) n
)))
592 (if (and (car num-root
) (car den-root
))
593 (list 'frac
(cdr num-root
) (cdr den-root
))
594 (and (not calc-symbolic-mode
)
597 (math-div-float (math-float (cdr num-root
))
598 (math-float (cdr den-root
))))))))
599 ((eq (car-safe a
) 'float
)
600 (and (not calc-symbolic-mode
)
601 (math-nth-root-float a n
)))
602 ((eq (car-safe a
) 'polar
)
603 (let ((root (math-nth-root (nth 1 a
) n
)))
604 (and root
(list 'polar root
(math-div (nth 2 a
) n
)))))
607 ;; The variables math-nrf-n, math-nrf-nf and math-nrf-nfm1 are local
608 ;; to math-nth-root-float, but are used by math-nth-root-float-iter,
609 ;; which is called by math-nth-root-float.
612 (defvar math-nrf-nfm1
)
614 (defun math-nth-root-float (a math-nrf-n
&optional guess
)
615 (math-inexact-result)
616 (math-with-extra-prec 1
617 (let ((math-nrf-nf (math-float math-nrf-n
))
618 (math-nrf-nfm1 (math-float (1- math-nrf-n
))))
619 (math-nth-root-float-iter a
(or guess
621 1 (/ (+ (math-numdigs (nth 1 a
))
626 (defun math-nth-root-float-iter (a guess
)
627 (math-working "root" guess
)
628 (let ((g2 (math-div-float (math-add-float (math-mul math-nrf-nfm1 guess
)
630 a
(math-ipow guess
(1- math-nrf-n
))))
632 (if (math-nearly-equal-float g2 guess
)
634 (math-nth-root-float-iter a g2
))))
636 ;; The variable math-nri-n is local to math-nth-root-integer, but
637 ;; is used by math-nth-root-int-iter, which is called by
638 ;; math-nth-root-int.
641 (defun math-nth-root-integer (a math-nri-n
&optional guess
) ; [I I S]
642 (math-nth-root-int-iter a
(or guess
643 (math-scale-int 1 (/ (+ (math-numdigs a
)
647 (defun math-nth-root-int-iter (a guess
)
648 (math-working "root" guess
)
649 (let* ((q (math-idivmod a
(math-ipow guess
(1- math-nri-n
))))
650 (s (math-add (car q
) (math-mul (1- math-nri-n
) guess
)))
651 (g2 (math-idivmod s math-nri-n
)))
652 (if (Math-natnum-lessp (car g2
) guess
)
653 (math-nth-root-int-iter a
(car g2
))
654 (cons (and (equal (car g2
) guess
)
659 (defun calcFunc-nroot (x n
)
660 (calcFunc-pow x
(if (integerp n
)
667 ;;;; Transcendental functions.
669 ;;; All of these functions are defined on the complex plane.
670 ;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
672 ;;; Most functions increase calc-internal-prec by 2 digits, then round
673 ;;; down afterward. "-raw" functions use the current precision, require
674 ;;; their arguments to be in float (or complex float) format, and always
675 ;;; work in radians (where applicable).
677 (defun math-to-radians (a) ; [N N]
678 (cond ((eq (car-safe a
) 'hms
)
679 (math-from-hms a
'rad
))
680 ((memq calc-angle-mode
'(deg hms
))
681 (math-mul a
(math-pi-over-180)))
684 (defun math-from-radians (a) ; [N N]
685 (cond ((eq calc-angle-mode
'deg
)
687 (math-div a
(math-pi-over-180))
688 (list 'calcFunc-deg a
)))
689 ((eq calc-angle-mode
'hms
)
690 (math-to-hms a
'rad
))
693 (defun math-to-radians-2 (a) ; [N N]
694 (cond ((eq (car-safe a
) 'hms
)
695 (math-from-hms a
'rad
))
696 ((memq calc-angle-mode
'(deg hms
))
697 (if calc-symbolic-mode
698 (math-div (math-mul a
'(var pi var-pi
)) 180)
699 (math-mul a
(math-pi-over-180))))
702 (defun math-from-radians-2 (a) ; [N N]
703 (cond ((memq calc-angle-mode
'(deg hms
))
704 (if calc-symbolic-mode
705 (math-div (math-mul 180 a
) '(var pi var-pi
))
706 (math-div a
(math-pi-over-180))))
711 ;;; Sine, cosine, and tangent.
713 (defun calcFunc-sin (x) ; [N N] [Public]
714 (cond ((and (integerp x
)
715 (if (eq calc-angle-mode
'deg
)
718 (aref [0 1 0 -
1] (math-mod (/ x
90) 4)))
720 (math-with-extra-prec 2
721 (math-sin-raw (math-to-radians (math-float x
)))))
724 (math-with-extra-prec 2
725 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
726 (xs (math-to-radians (math-float (nth 2 x
))))
727 (sc (math-sin-cos-raw xx
)))
728 (math-make-sdev (car sc
) (math-mul xs
(cdr sc
)))))
729 (math-make-sdev (calcFunc-sin (nth 1 x
))
730 (math-mul (nth 2 x
) (calcFunc-cos (nth 1 x
))))))
731 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
732 (calcFunc-cos (math-sub x
(math-quarter-circle nil
))))
733 ((equal x
'(var nan var-nan
))
735 (t (calc-record-why 'scalarp x
)
736 (list 'calcFunc-sin x
))))
738 (defun calcFunc-cos (x) ; [N N] [Public]
739 (cond ((and (integerp x
)
740 (if (eq calc-angle-mode
'deg
)
743 (aref [1 0 -
1 0] (math-mod (/ x
90) 4)))
745 (math-with-extra-prec 2
746 (math-cos-raw (math-to-radians (math-float x
)))))
749 (math-with-extra-prec 2
750 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
751 (xs (math-to-radians (math-float (nth 2 x
))))
752 (sc (math-sin-cos-raw xx
)))
753 (math-make-sdev (cdr sc
) (math-mul xs
(car sc
)))))
754 (math-make-sdev (calcFunc-cos (nth 1 x
))
755 (math-mul (nth 2 x
) (calcFunc-sin (nth 1 x
))))))
756 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
757 (math-with-extra-prec 2
758 (let* ((xx (math-to-radians (math-float x
)))
759 (na (math-floor (math-div (nth 2 xx
) (math-pi))))
760 (nb (math-floor (math-div (nth 3 xx
) (math-pi))))
761 (span (math-sub nb na
)))
762 (if (memq span
'(0 1))
763 (let ((int (math-sort-intv (nth 1 x
)
764 (math-cos-raw (nth 2 xx
))
765 (math-cos-raw (nth 3 xx
)))))
768 (math-make-intv (logior (nth 1 x
) 2)
771 (math-make-intv (logior (nth 1 x
) 1)
775 (list 'intv
3 -
1 1)))))
776 ((equal x
'(var nan var-nan
))
778 (t (calc-record-why 'scalarp x
)
779 (list 'calcFunc-cos x
))))
781 (defun calcFunc-sincos (x) ; [V N] [Public]
783 (math-with-extra-prec 2
784 (let ((sc (math-sin-cos-raw (math-to-radians (math-float x
)))))
785 (list 'vec
(cdr sc
) (car sc
)))) ; the vector [cos, sin]
786 (list 'vec
(calcFunc-sin x
) (calcFunc-cos x
))))
788 (defun calcFunc-tan (x) ; [N N] [Public]
789 (cond ((and (integerp x
)
790 (if (eq calc-angle-mode
'deg
)
795 (math-with-extra-prec 2
796 (math-tan-raw (math-to-radians (math-float x
)))))
799 (math-with-extra-prec 2
800 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
801 (xs (math-to-radians (math-float (nth 2 x
))))
802 (sc (math-sin-cos-raw xx
)))
803 (if (and (math-zerop (cdr sc
)) (not calc-infinite-mode
))
805 (calc-record-why "*Division by zero")
806 (list 'calcFunc-tan x
))
807 (math-make-sdev (math-div-float (car sc
) (cdr sc
))
808 (math-div-float xs
(math-sqr (cdr sc
)))))))
809 (math-make-sdev (calcFunc-tan (nth 1 x
))
811 (math-sqr (calcFunc-cos (nth 1 x
)))))))
812 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
813 (or (math-with-extra-prec 2
814 (let* ((xx (math-to-radians (math-float x
)))
815 (na (math-floor (math-div (math-sub (nth 2 xx
)
818 (nb (math-floor (math-div (math-sub (nth 3 xx
)
822 (math-sort-intv (nth 1 x
)
823 (math-tan-raw (nth 2 xx
))
824 (math-tan-raw (nth 3 xx
))))))
825 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))))
826 ((equal x
'(var nan var-nan
))
828 (t (calc-record-why 'scalarp x
)
829 (list 'calcFunc-tan x
))))
831 (defun calcFunc-sec (x)
832 (cond ((and (integerp x
)
833 (eq calc-angle-mode
'deg
)
839 (eq calc-angle-mode
'rad
)
843 (math-with-extra-prec 2
844 (math-sec-raw (math-to-radians (math-float x
)))))
847 (math-with-extra-prec 2
848 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
849 (xs (math-to-radians (math-float (nth 2 x
))))
850 (sc (math-sin-cos-raw xx
)))
851 (if (and (math-zerop (cdr sc
))
852 (not calc-infinite-mode
))
854 (calc-record-why "*Division by zero")
855 (list 'calcFunc-sec x
))
856 (math-make-sdev (math-div-float '(float 1 0) (cdr sc
))
858 (math-mul xs
(car sc
))
859 (math-sqr (cdr sc
)))))))
860 (math-make-sdev (calcFunc-sec (nth 1 x
))
863 (calcFunc-sin (nth 1 x
)))
864 (math-sqr (calcFunc-cos (nth 1 x
)))))))
865 ((and (eq (car x
) 'intv
)
866 (math-intv-constp x
))
867 (math-with-extra-prec 2
868 (let* ((xx (math-to-radians (math-float x
)))
869 (na (math-floor (math-div (math-sub (nth 2 xx
)
872 (nb (math-floor (math-div (math-sub (nth 3 xx
)
875 (naa (math-floor (math-div (nth 2 xx
) (math-pi-over-2))))
876 (nbb (math-floor (math-div (nth 3 xx
) (math-pi-over-2))))
877 (span (math-sub nbb naa
)))
878 (if (not (equal na nb
))
879 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
880 (let ((int (math-sort-intv (nth 1 x
)
881 (math-sec-raw (nth 2 xx
))
882 (math-sec-raw (nth 3 xx
)))))
884 (if (math-evenp (math-div (math-add naa
1) 2))
885 (math-make-intv (logior (nth 1 int
) 2)
888 (math-make-intv (logior (nth 1 int
) 1)
892 ((equal x
'(var nan var-nan
))
894 (t (calc-record-why 'scalarp x
)
895 (list 'calcFunc-sec x
))))
897 (defun calcFunc-csc (x)
898 (cond ((and (integerp x
)
899 (eq calc-angle-mode
'deg
)
900 (= (%
(- x
90) 180) 0))
901 (if (= (%
(- x
90) 360) 0)
905 (math-with-extra-prec 2
906 (math-csc-raw (math-to-radians (math-float x
)))))
909 (math-with-extra-prec 2
910 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
911 (xs (math-to-radians (math-float (nth 2 x
))))
912 (sc (math-sin-cos-raw xx
)))
913 (if (and (math-zerop (car sc
))
914 (not calc-infinite-mode
))
916 (calc-record-why "*Division by zero")
917 (list 'calcFunc-csc x
))
918 (math-make-sdev (math-div-float '(float 1 0) (car sc
))
920 (math-mul xs
(cdr sc
))
921 (math-sqr (car sc
)))))))
922 (math-make-sdev (calcFunc-csc (nth 1 x
))
925 (calcFunc-cos (nth 1 x
)))
926 (math-sqr (calcFunc-sin (nth 1 x
)))))))
927 ((and (eq (car x
) 'intv
)
928 (math-intv-constp x
))
929 (math-with-extra-prec 2
930 (let* ((xx (math-to-radians (math-float x
)))
931 (na (math-floor (math-div (nth 2 xx
) (math-pi))))
932 (nb (math-floor (math-div (nth 3 xx
) (math-pi))))
933 (naa (math-floor (math-div (nth 2 xx
) (math-pi-over-2))))
934 (nbb (math-floor (math-div (nth 3 xx
) (math-pi-over-2))))
935 (span (math-sub nbb naa
)))
936 (if (not (equal na nb
))
937 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
938 (let ((int (math-sort-intv (nth 1 x
)
939 (math-csc-raw (nth 2 xx
))
940 (math-csc-raw (nth 3 xx
)))))
942 (if (math-evenp (math-div naa
2))
943 (math-make-intv (logior (nth 1 int
) 2)
946 (math-make-intv (logior (nth 1 int
) 1)
950 ((equal x
'(var nan var-nan
))
952 (t (calc-record-why 'scalarp x
)
953 (list 'calcFunc-csc x
))))
955 (defun calcFunc-cot (x) ; [N N] [Public]
956 (cond ((and (integerp x
)
957 (if (eq calc-angle-mode
'deg
)
958 (= (%
(- x
90) 180) 0)
962 (math-with-extra-prec 2
963 (math-cot-raw (math-to-radians (math-float x
)))))
966 (math-with-extra-prec 2
967 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
968 (xs (math-to-radians (math-float (nth 2 x
))))
969 (sc (math-sin-cos-raw xx
)))
970 (if (and (math-zerop (car sc
)) (not calc-infinite-mode
))
972 (calc-record-why "*Division by zero")
973 (list 'calcFunc-cot x
))
974 (math-make-sdev (math-div-float (cdr sc
) (car sc
))
975 (math-div-float xs
(math-sqr (car sc
)))))))
976 (math-make-sdev (calcFunc-cot (nth 1 x
))
978 (math-sqr (calcFunc-sin (nth 1 x
)))))))
979 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
980 (or (math-with-extra-prec 2
981 (let* ((xx (math-to-radians (math-float x
)))
982 (na (math-floor (math-div (nth 2 xx
) (math-pi))))
983 (nb (math-floor (math-div (nth 3 xx
) (math-pi)))))
985 (math-sort-intv (nth 1 x
)
986 (math-cot-raw (nth 2 xx
))
987 (math-cot-raw (nth 3 xx
))))))
988 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))))
989 ((equal x
'(var nan var-nan
))
991 (t (calc-record-why 'scalarp x
)
992 (list 'calcFunc-cot x
))))
994 (defun math-sin-raw (x) ; [N N]
995 (cond ((eq (car x
) 'cplx
)
996 (let* ((expx (math-exp-raw (nth 2 x
)))
997 (expmx (math-div-float '(float 1 0) expx
))
998 (sc (math-sin-cos-raw (nth 1 x
))))
1000 (math-mul-float (car sc
)
1001 (math-mul-float (math-add-float expx expmx
)
1003 (math-mul-float (cdr sc
)
1004 (math-mul-float (math-sub-float expx expmx
)
1006 ((eq (car x
) 'polar
)
1007 (math-polar (math-sin-raw (math-complex x
))))
1008 ((Math-integer-negp (nth 1 x
))
1009 (math-neg-float (math-sin-raw (math-neg-float x
))))
1010 ((math-lessp-float '(float 7 0) x
) ; avoid inf loops due to roundoff
1011 (math-sin-raw (math-mod x
(math-two-pi))))
1012 (t (math-sin-raw-2 x x
))))
1014 (defun math-cos-raw (x) ; [N N]
1015 (if (eq (car-safe x
) 'polar
)
1016 (math-polar (math-cos-raw (math-complex x
)))
1017 (math-sin-raw (math-sub (math-pi-over-2) x
))))
1019 (defun math-sec-raw (x) ; [N N]
1020 (cond ((eq (car x
) 'cplx
)
1021 (let* ((x (math-mul x
'(float 1 0)))
1022 (expx (math-exp-raw (nth 2 x
)))
1023 (expmx (math-div-float '(float 1 0) expx
))
1024 (sh (math-mul-float (math-sub-float expx expmx
) '(float 5 -
1)))
1025 (ch (math-mul-float (math-add-float expx expmx
) '(float 5 -
1)))
1026 (sc (math-sin-cos-raw (nth 1 x
)))
1028 (math-mul-float (math-sqr (car sc
))
1030 (math-mul-float (math-sqr (cdr sc
))
1032 (and (not (eq (nth 1 d
) 0))
1034 (math-div-float (math-mul-float (cdr sc
) ch
) d
)
1035 (math-div-float (math-mul-float (car sc
) sh
) d
)))))
1036 ((eq (car x
) 'polar
)
1037 (math-polar (math-sec-raw (math-complex x
))))
1039 (let ((cs (math-cos-raw x
)))
1042 (math-div-float '(float 1 0) cs
))))))
1044 (defun math-csc-raw (x) ; [N N]
1045 (cond ((eq (car x
) 'cplx
)
1046 (let* ((x (math-mul x
'(float 1 0)))
1047 (expx (math-exp-raw (nth 2 x
)))
1048 (expmx (math-div-float '(float 1 0) expx
))
1049 (sh (math-mul-float (math-sub-float expx expmx
) '(float 5 -
1)))
1050 (ch (math-mul-float (math-add-float expx expmx
) '(float 5 -
1)))
1051 (sc (math-sin-cos-raw (nth 1 x
)))
1053 (math-mul-float (math-sqr (car sc
))
1055 (math-mul-float (math-sqr (cdr sc
))
1057 (and (not (eq (nth 1 d
) 0))
1059 (math-div-float (math-mul-float (car sc
) ch
) d
)
1060 (math-div-float (math-mul-float (cdr sc
) sh
) d
)))))
1061 ((eq (car x
) 'polar
)
1062 (math-polar (math-csc-raw (math-complex x
))))
1064 (let ((sn (math-sin-raw x
)))
1067 (math-div-float '(float 1 0) sn
))))))
1069 (defun math-cot-raw (x) ; [N N]
1070 (cond ((eq (car x
) 'cplx
)
1071 (let* ((x (math-mul x
'(float 1 0)))
1072 (expx (math-exp-raw (nth 2 x
)))
1073 (expmx (math-div-float '(float 1 0) expx
))
1074 (sh (math-mul-float (math-sub-float expx expmx
) '(float 5 -
1)))
1075 (ch (math-mul-float (math-add-float expx expmx
) '(float 5 -
1)))
1076 (sc (math-sin-cos-raw (nth 1 x
)))
1080 (and (not (eq (nth 1 d
) 0))
1083 (math-mul-float (car sc
) (cdr sc
))
1087 (math-mul-float sh ch
)
1089 ((eq (car x
) 'polar
)
1090 (math-polar (math-cot-raw (math-complex x
))))
1092 (let ((sc (math-sin-cos-raw x
)))
1093 (if (eq (nth 1 (car sc
)) 0)
1094 (math-div (cdr sc
) 0)
1095 (math-div-float (cdr sc
) (car sc
)))))))
1098 ;;; This could use a smarter method: Reduce x as in math-sin-raw, then
1099 ;;; compute either sin(x) or cos(x), whichever is smaller, and compute
1100 ;;; the other using the identity sin(x)^2 + cos(x)^2 = 1.
1101 (defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x))
1102 (cons (math-sin-raw x
) (math-cos-raw x
)))
1104 (defun math-tan-raw (x) ; [N N]
1105 (cond ((eq (car x
) 'cplx
)
1106 (let* ((x (math-mul x
'(float 2 0)))
1107 (expx (math-exp-raw (nth 2 x
)))
1108 (expmx (math-div-float '(float 1 0) expx
))
1109 (sc (math-sin-cos-raw (nth 1 x
)))
1110 (d (math-add-float (cdr sc
)
1111 (math-mul-float (math-add-float expx expmx
)
1113 (and (not (eq (nth 1 d
) 0))
1115 (math-div-float (car sc
) d
)
1116 (math-div-float (math-mul-float (math-sub-float expx
1118 '(float 5 -
1)) d
)))))
1119 ((eq (car x
) 'polar
)
1120 (math-polar (math-tan-raw (math-complex x
))))
1122 (let ((sc (math-sin-cos-raw x
)))
1123 (if (eq (nth 1 (cdr sc
)) 0)
1124 (math-div (car sc
) 0)
1125 (math-div-float (car sc
) (cdr sc
)))))))
1127 (defun math-sin-raw-2 (x orgx
) ; This avoids poss of inf recursion. [F F]
1128 (let ((xmpo2 (math-sub-float (math-pi-over-2) x
)))
1129 (cond ((Math-integer-negp (nth 1 xmpo2
))
1130 (math-neg-float (math-sin-raw-2 (math-sub-float x
(math-pi))
1132 ((math-lessp-float (math-pi-over-4) x
)
1133 (math-cos-raw-2 xmpo2 orgx
))
1134 ((math-lessp-float x
(math-neg (math-pi-over-4)))
1135 (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x
) orgx
)))
1136 ((math-nearly-zerop-float x orgx
) '(float 0 0))
1137 (calc-symbolic-mode (signal 'inexact-result nil
))
1138 (t (math-sin-series x
6 4 x
(math-neg-float (math-sqr-float x
)))))))
1140 (defun math-cos-raw-2 (x orgx
) ; [F F]
1141 (cond ((math-nearly-zerop-float x orgx
) '(float 1 0))
1142 (calc-symbolic-mode (signal 'inexact-result nil
))
1143 (t (let ((xnegsqr (math-neg-float (math-sqr-float x
))))
1145 (math-add-float '(float 1 0)
1146 (math-mul-float xnegsqr
'(float 5 -
1)))
1147 24 5 xnegsqr xnegsqr
)))))
1149 (defun math-sin-series (sum nfac n x xnegsqr
)
1150 (math-working "sin" sum
)
1151 (let* ((nextx (math-mul-float x xnegsqr
))
1152 (nextsum (math-add-float sum
(math-div-float nextx
1153 (math-float nfac
)))))
1154 (if (math-nearly-equal-float sum nextsum
)
1156 (math-sin-series nextsum
(math-mul nfac
(* n
(1+ n
)))
1157 (+ n
2) nextx xnegsqr
))))
1160 ;;; Inverse sine, cosine, tangent.
1162 (defun calcFunc-arcsin (x) ; [N N] [Public]
1164 ((and (eq x
1) (eq calc-angle-mode
'deg
)) 90)
1165 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) -
90)
1166 (calc-symbolic-mode (signal 'inexact-result nil
))
1168 (math-with-extra-prec 2
1169 (math-from-radians (math-arcsin-raw (math-float x
)))))
1171 (math-make-sdev (calcFunc-arcsin (nth 1 x
))
1175 (math-sub 1 (math-sqr (nth 1 x
))))))))
1177 (math-sort-intv (nth 1 x
)
1178 (calcFunc-arcsin (nth 2 x
))
1179 (calcFunc-arcsin (nth 3 x
))))
1180 ((equal x
'(var nan var-nan
))
1182 (t (calc-record-why 'numberp x
)
1183 (list 'calcFunc-arcsin x
))))
1185 (defun calcFunc-arccos (x) ; [N N] [Public]
1187 ((and (eq x
0) (eq calc-angle-mode
'deg
)) 90)
1188 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) 180)
1189 (calc-symbolic-mode (signal 'inexact-result nil
))
1191 (math-with-extra-prec 2
1192 (math-from-radians (math-arccos-raw (math-float x
)))))
1194 (math-make-sdev (calcFunc-arccos (nth 1 x
))
1198 (math-sub 1 (math-sqr (nth 1 x
))))))))
1200 (math-sort-intv (nth 1 x
)
1201 (calcFunc-arccos (nth 2 x
))
1202 (calcFunc-arccos (nth 3 x
))))
1203 ((equal x
'(var nan var-nan
))
1205 (t (calc-record-why 'numberp x
)
1206 (list 'calcFunc-arccos x
))))
1208 (defun calcFunc-arctan (x) ; [N N] [Public]
1210 ((and (eq x
1) (eq calc-angle-mode
'deg
)) 45)
1211 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) -
45)
1213 (math-with-extra-prec 2
1214 (math-from-radians (math-arctan-raw (math-float x
)))))
1216 (math-make-sdev (calcFunc-arctan (nth 1 x
))
1219 (math-add 1 (math-sqr (nth 1 x
)))))))
1221 (math-sort-intv (nth 1 x
)
1222 (calcFunc-arctan (nth 2 x
))
1223 (calcFunc-arctan (nth 3 x
))))
1224 ((equal x
'(var inf var-inf
))
1225 (math-quarter-circle t
))
1226 ((equal x
'(neg (var inf var-inf
)))
1227 (math-neg (math-quarter-circle t
)))
1228 ((equal x
'(var nan var-nan
))
1230 (t (calc-record-why 'numberp x
)
1231 (list 'calcFunc-arctan x
))))
1233 (defun math-arcsin-raw (x) ; [N N]
1234 (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x
)))))
1235 (if (or (memq (car x
) '(cplx polar
))
1236 (memq (car a
) '(cplx polar
)))
1237 (math-with-extra-prec 2 ; use extra precision for difficult case
1238 (math-mul '(cplx 0 -
1)
1239 (math-ln-raw (math-add (math-mul '(cplx 0 1) x
) a
))))
1240 (math-arctan2-raw x a
))))
1242 (defun math-arccos-raw (x) ; [N N]
1243 (math-sub (math-pi-over-2) (math-arcsin-raw x
)))
1245 (defun math-arctan-raw (x) ; [N N]
1246 (cond ((memq (car x
) '(cplx polar
))
1247 (math-with-extra-prec 2 ; extra-extra
1249 (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x
)))
1250 (math-ln-raw (math-add 1 (math-mul '(cplx 0 -
1) x
))))
1252 ((Math-integer-negp (nth 1 x
))
1253 (math-neg-float (math-arctan-raw (math-neg-float x
))))
1255 (calc-symbolic-mode (signal 'inexact-result nil
))
1256 ((math-equal-int x
1) (math-pi-over-4))
1257 ((math-equal-int x -
1) (math-neg (math-pi-over-4)))
1258 ((math-lessp-float '(float 414214 -
6) x
) ; if x > sqrt(2) - 1, reduce
1259 (if (math-lessp-float '(float 1 0) x
)
1260 (math-sub-float (math-mul-float (math-pi) '(float 5 -
1))
1261 (math-arctan-raw (math-div-float '(float 1 0) x
)))
1262 (math-sub-float (math-mul-float (math-pi) '(float 25 -
2))
1263 (math-arctan-raw (math-div-float
1264 (math-sub-float '(float 1 0) x
)
1265 (math-add-float '(float 1 0)
1267 (t (math-arctan-series x
3 x
(math-neg-float (math-sqr-float x
))))))
1269 (defun math-arctan-series (sum n x xnegsqr
)
1270 (math-working "arctan" sum
)
1271 (let* ((nextx (math-mul-float x xnegsqr
))
1272 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1273 (if (math-nearly-equal-float sum nextsum
)
1275 (math-arctan-series nextsum
(+ n
2) nextx xnegsqr
))))
1277 (defun calcFunc-arctan2 (y x
) ; [F R R] [Public]
1280 (math-with-extra-prec 2
1281 (math-from-radians (math-arctan2-raw (math-float y
)
1283 (calc-record-why 'anglep x
)
1284 (list 'calcFunc-arctan2 y x
))
1285 (if (and (or (math-infinitep x
) (math-anglep x
))
1286 (or (math-infinitep y
) (math-anglep y
)))
1301 (calcFunc-arctan2 y x
)
1302 '(var nan var-nan
)))
1303 (calc-record-why 'anglep y
)
1304 (list 'calcFunc-arctan2 y x
))))
1306 (defun math-arctan2-raw (y x
) ; [F R R]
1307 (cond ((math-zerop y
)
1308 (if (math-negp x
) (math-pi)
1309 (if (or (math-floatp x
) (math-floatp y
)) '(float 0 0) 0)))
1313 (math-neg (math-pi-over-2))))
1315 (math-arctan-raw (math-div-float y x
)))
1317 (math-add-float (math-arctan-raw (math-div-float y x
))
1320 (math-sub-float (math-arctan-raw (math-div-float y x
))
1323 (defun calcFunc-arcsincos (x) ; [V N] [Public]
1324 (if (and (Math-vectorp x
)
1326 (calcFunc-arctan2 (nth 2 x
) (nth 1 x
))
1327 (math-reject-arg x
"*Two-element vector expected")))
1331 ;;; Exponential function.
1333 (defun calcFunc-exp (x) ; [N N] [Public]
1335 ((and (memq x
'(1 -
1)) calc-symbolic-mode
)
1336 (if (eq x
1) '(var e var-e
) (math-div 1 '(var e var-e
))))
1338 (math-with-extra-prec 2 (math-exp-raw (math-float x
))))
1339 ((eq (car-safe x
) 'sdev
)
1340 (let ((ex (calcFunc-exp (nth 1 x
))))
1341 (math-make-sdev ex
(math-mul (nth 2 x
) ex
))))
1342 ((eq (car-safe x
) 'intv
)
1343 (math-make-intv (nth 1 x
) (calcFunc-exp (nth 2 x
))
1344 (calcFunc-exp (nth 3 x
))))
1345 ((equal x
'(var inf var-inf
))
1347 ((equal x
'(neg (var inf var-inf
)))
1349 ((equal x
'(var nan var-nan
))
1351 (t (calc-record-why 'numberp x
)
1352 (list 'calcFunc-exp x
))))
1354 (defun calcFunc-expm1 (x) ; [N N] [Public]
1356 ((math-zerop x
) '(float 0 0))
1357 (calc-symbolic-mode (signal 'inexact-result nil
))
1359 (math-with-extra-prec 2
1360 (let ((x (math-float x
)))
1361 (if (and (eq (car x
) 'float
)
1362 (math-lessp-float x
'(float 1 0))
1363 (math-lessp-float '(float -
1 0) x
))
1364 (math-exp-minus-1-raw x
)
1365 (math-add (math-exp-raw x
) -
1)))))
1366 ((eq (car-safe x
) 'sdev
)
1368 (let ((ex (calcFunc-expm1 (nth 1 x
))))
1369 (math-make-sdev ex
(math-mul (nth 2 x
) (math-add ex
1))))
1370 (math-make-sdev (calcFunc-expm1 (nth 1 x
))
1371 (math-mul (nth 2 x
) (calcFunc-exp (nth 1 x
))))))
1372 ((eq (car-safe x
) 'intv
)
1373 (math-make-intv (nth 1 x
)
1374 (calcFunc-expm1 (nth 2 x
))
1375 (calcFunc-expm1 (nth 3 x
))))
1376 ((equal x
'(var inf var-inf
))
1378 ((equal x
'(neg (var inf var-inf
)))
1380 ((equal x
'(var nan var-nan
))
1382 (t (calc-record-why 'numberp x
)
1383 (list 'calcFunc-expm1 x
))))
1385 (defun calcFunc-exp10 (x) ; [N N] [Public]
1388 (math-pow '(float 1 1) x
)))
1390 (defun math-exp-raw (x) ; [N N]
1391 (cond ((math-zerop x
) '(float 1 0))
1392 (calc-symbolic-mode (signal 'inexact-result nil
))
1394 (let ((expx (math-exp-raw (nth 1 x
)))
1395 (sc (math-sin-cos-raw (nth 2 x
))))
1397 (math-mul-float expx
(cdr sc
))
1398 (math-mul-float expx
(car sc
)))))
1399 ((eq (car x
) 'polar
)
1400 (let ((xc (math-complex x
)))
1402 (math-exp-raw (nth 1 xc
))
1403 (math-from-radians (nth 2 xc
)))))
1404 ((or (math-lessp-float '(float 5 -
1) x
)
1405 (math-lessp-float x
'(float -
5 -
1)))
1406 (if (math-lessp-float '(float 921035 1) x
)
1408 (if (math-lessp-float x
'(float -
921035 1))
1410 (let* ((two-x (math-mul-float x
'(float 2 0)))
1411 (hint (math-scale-int (nth 1 two-x
) (nth 2 two-x
)))
1412 (hfrac (math-sub-float x
(math-mul-float (math-float hint
)
1414 (math-mul-float (math-ipow (math-sqrt-e) hint
)
1415 (math-add-float '(float 1 0)
1416 (math-exp-minus-1-raw hfrac
)))))
1417 (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x
)))))
1419 (defun math-exp-minus-1-raw (x) ; [F F]
1420 (math-exp-series x
2 3 x x
))
1422 (defun math-exp-series (sum nfac n xpow x
)
1423 (math-working "exp" sum
)
1424 (let* ((nextx (math-mul-float xpow x
))
1425 (nextsum (math-add-float sum
(math-div-float nextx
1426 (math-float nfac
)))))
1427 (if (math-nearly-equal-float sum nextsum
)
1429 (math-exp-series nextsum
(math-mul nfac n
) (1+ n
) nextx x
))))
1435 (defun calcFunc-ln (x) ; [N N] [Public]
1436 (cond ((math-zerop x
)
1437 (if calc-infinite-mode
1438 '(neg (var inf var-inf
))
1439 (math-reject-arg x
"*Logarithm of zero")))
1442 (math-with-extra-prec 2 (math-ln-raw (math-float x
))))
1443 ((eq (car-safe x
) 'sdev
)
1444 (math-make-sdev (calcFunc-ln (nth 1 x
))
1445 (math-div (nth 2 x
) (nth 1 x
))))
1446 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1447 (Math-zerop (nth 2 x
))
1448 (not (math-intv-constp x
))))
1449 (let ((calc-infinite-mode t
))
1450 (math-make-intv (nth 1 x
) (calcFunc-ln (nth 2 x
))
1451 (calcFunc-ln (nth 3 x
)))))
1452 ((equal x
'(var e var-e
))
1454 ((and (eq (car-safe x
) '^
)
1455 (equal (nth 1 x
) '(var e var-e
))
1456 (math-known-realp (nth 2 x
)))
1459 (if (equal x
'(var nan var-nan
))
1461 '(var inf var-inf
)))
1462 (t (calc-record-why 'numberp x
)
1463 (list 'calcFunc-ln x
))))
1465 (defun calcFunc-log10 (x) ; [N N] [Public]
1466 (cond ((math-equal-int x
1)
1467 (if (math-floatp x
) '(float 0 0) 0))
1468 ((and (Math-integerp x
)
1470 (let ((res (math-integer-log x
10)))
1472 (setq x
(cdr res
)))))
1474 ((and (eq (car-safe x
) 'frac
)
1476 (let ((res (math-integer-log (nth 2 x
) 10)))
1478 (setq x
(- (cdr res
))))))
1481 (if calc-infinite-mode
1482 '(neg (var inf var-inf
))
1483 (math-reject-arg x
"*Logarithm of zero")))
1484 (calc-symbolic-mode (signal 'inexact-result nil
))
1486 (math-with-extra-prec 2
1487 (let ((xf (math-float x
)))
1488 (if (eq (nth 1 xf
) 0)
1489 (math-reject-arg x
"*Logarithm of zero"))
1490 (if (Math-integer-posp (nth 1 xf
))
1491 (if (eq (nth 1 xf
) 1) ; log10(1*10^n) = n
1492 (math-float (nth 2 xf
))
1493 (let ((xdigs (1- (math-numdigs (nth 1 xf
)))))
1495 (math-div-float (math-ln-raw-2
1496 (list 'float
(nth 1 xf
) (- xdigs
)))
1498 (math-float (+ (nth 2 xf
) xdigs
)))))
1499 (math-div (calcFunc-ln xf
) (math-ln-10))))))
1500 ((eq (car-safe x
) 'sdev
)
1501 (math-make-sdev (calcFunc-log10 (nth 1 x
))
1503 (math-mul (nth 1 x
) (math-ln-10)))))
1504 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1505 (not (math-intv-constp x
))))
1506 (math-make-intv (nth 1 x
)
1507 (calcFunc-log10 (nth 2 x
))
1508 (calcFunc-log10 (nth 3 x
))))
1510 (if (equal x
'(var nan var-nan
))
1512 '(var inf var-inf
)))
1513 (t (calc-record-why 'numberp x
)
1514 (list 'calcFunc-log10 x
))))
1516 (defun calcFunc-log (x &optional b
) ; [N N N] [Public]
1517 (cond ((or (null b
) (equal b
'(var e var-e
)))
1518 (math-normalize (list 'calcFunc-ln x
)))
1519 ((or (eq b
10) (equal b
'(float 1 1)))
1520 (math-normalize (list 'calcFunc-log10 x
)))
1522 (if calc-infinite-mode
1523 (math-div (calcFunc-ln x
) (calcFunc-ln b
))
1524 (math-reject-arg x
"*Logarithm of zero")))
1526 (if calc-infinite-mode
1527 (math-div (calcFunc-ln x
) (calcFunc-ln b
))
1528 (math-reject-arg b
"*Logarithm of zero")))
1529 ((math-equal-int b
1)
1530 (if calc-infinite-mode
1531 (math-div (calcFunc-ln x
) 0)
1532 (math-reject-arg b
"*Logarithm base one")))
1533 ((math-equal-int x
1)
1534 (if (math-floatp b
) '(float 0 0) 0))
1535 ((and (Math-ratp x
) (Math-ratp b
)
1536 (math-posp x
) (math-posp b
)
1537 (let* ((sign 1) (inv nil
)
1538 (xx (if (Math-lessp 1 x
)
1542 (bb (if (Math-lessp 1 b
)
1544 (setq sign
(- sign
))
1546 (res (if (Math-lessp xx bb
)
1547 (setq inv
(math-integer-log bb xx
))
1548 (math-integer-log xx bb
))))
1551 (math-div 1 (* sign
(cdr res
)))
1552 (* sign
(cdr res
)))))))
1554 (calc-symbolic-mode (signal 'inexact-result nil
))
1555 ((and (Math-numberp x
) (Math-numberp b
))
1556 (math-with-extra-prec 2
1557 (math-div (math-ln-raw (math-float x
))
1558 (math-log-base-raw b
))))
1559 ((and (eq (car-safe x
) 'sdev
)
1561 (math-make-sdev (calcFunc-log (nth 1 x
) b
)
1564 (math-log-base-raw b
)))))
1565 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1566 (not (math-intv-constp x
)))
1568 (math-make-intv (nth 1 x
)
1569 (calcFunc-log (nth 2 x
) b
)
1570 (calcFunc-log (nth 3 x
) b
)))
1571 ((or (eq (car-safe x
) 'intv
) (eq (car-safe b
) 'intv
))
1572 (math-div (calcFunc-ln x
) (calcFunc-ln b
)))
1573 ((or (math-infinitep x
)
1575 (math-div (calcFunc-ln x
) (calcFunc-ln b
)))
1576 (t (if (Math-numberp b
)
1577 (calc-record-why 'numberp x
)
1578 (calc-record-why 'numberp b
))
1579 (list 'calcFunc-log x b
))))
1581 (defun calcFunc-alog (x &optional b
)
1582 (cond ((or (null b
) (equal b
'(var e var-e
)))
1583 (math-normalize (list 'calcFunc-exp x
)))
1584 (t (math-pow b x
))))
1586 (defun calcFunc-ilog (x b
)
1587 (if (and (math-natnump x
) (not (eq x
0))
1588 (math-natnump b
) (not (eq b
0)))
1590 (math-reject-arg x
"*Logarithm base one")
1591 (if (Math-natnum-lessp x b
)
1593 (cdr (math-integer-log x b
))))
1594 (math-floor (calcFunc-log x b
))))
1596 (defun math-integer-log (x b
)
1597 (let ((pows (list b
))
1601 (while (not (Math-lessp x pow
))
1602 (setq pows
(cons pow pows
)
1603 pow
(math-sqr pow
)))
1604 (setq n
(lsh 1 (1- (length pows
)))
1607 (while (and (setq pows
(cdr pows
))
1610 next
(math-mul pow
(car pows
)))
1611 (or (Math-lessp x next
)
1614 (cons (equal pow x
) sum
)))
1617 (defvar math-log-base-cache nil
)
1618 (defun math-log-base-raw (b) ; [N N]
1619 (if (not (and (equal (car math-log-base-cache
) b
)
1620 (eq (nth 1 math-log-base-cache
) calc-internal-prec
)))
1621 (setq math-log-base-cache
(list b calc-internal-prec
1622 (math-ln-raw (math-float b
)))))
1623 (nth 2 math-log-base-cache
))
1625 (defun calcFunc-lnp1 (x) ; [N N] [Public]
1626 (cond ((Math-equal-int x -
1)
1627 (if calc-infinite-mode
1628 '(neg (var inf var-inf
))
1629 (math-reject-arg x
"*Logarithm of zero")))
1631 ((math-zerop x
) '(float 0 0))
1632 (calc-symbolic-mode (signal 'inexact-result nil
))
1634 (math-with-extra-prec 2
1635 (let ((x (math-float x
)))
1636 (if (and (eq (car x
) 'float
)
1637 (math-lessp-float x
'(float 5 -
1))
1638 (math-lessp-float '(float -
5 -
1) x
))
1639 (math-ln-plus-1-raw x
)
1640 (math-ln-raw (math-add-float x
'(float 1 0)))))))
1641 ((eq (car-safe x
) 'sdev
)
1642 (math-make-sdev (calcFunc-lnp1 (nth 1 x
))
1643 (math-div (nth 2 x
) (math-add (nth 1 x
) 1))))
1644 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1645 (not (math-intv-constp x
))))
1646 (math-make-intv (nth 1 x
)
1647 (calcFunc-lnp1 (nth 2 x
))
1648 (calcFunc-lnp1 (nth 3 x
))))
1650 (if (equal x
'(var nan var-nan
))
1652 '(var inf var-inf
)))
1653 (t (calc-record-why 'numberp x
)
1654 (list 'calcFunc-lnp1 x
))))
1656 (defun math-ln-raw (x) ; [N N] --- must be float format!
1657 (cond ((eq (car-safe x
) 'cplx
)
1659 (math-mul-float (math-ln-raw
1660 (math-add-float (math-sqr-float (nth 1 x
))
1661 (math-sqr-float (nth 2 x
))))
1663 (math-arctan2-raw (nth 2 x
) (nth 1 x
))))
1664 ((eq (car x
) 'polar
)
1665 (math-polar (list 'cplx
1666 (math-ln-raw (nth 1 x
))
1667 (math-to-radians (nth 2 x
)))))
1668 ((Math-equal-int x
1)
1670 (calc-symbolic-mode (signal 'inexact-result nil
))
1671 ((math-posp (nth 1 x
)) ; positive and real
1672 (let ((xdigs (1- (math-numdigs (nth 1 x
)))))
1673 (math-add-float (math-ln-raw-2 (list 'float
(nth 1 x
) (- xdigs
)))
1674 (math-mul-float (math-float (+ (nth 2 x
) xdigs
))
1677 (math-reject-arg x
"*Logarithm of zero"))
1678 ((eq calc-complex-mode
'polar
) ; negative and real
1680 (list 'cplx
; negative and real
1681 (math-ln-raw (math-neg-float x
))
1683 (t (list 'cplx
; negative and real
1684 (math-ln-raw (math-neg-float x
))
1687 (defun math-ln-raw-2 (x) ; [F F]
1688 (cond ((math-lessp-float '(float 14 -
1) x
)
1689 (math-add-float (math-ln-raw-2 (math-mul-float x
'(float 5 -
1)))
1691 (t ; now .7 < x <= 1.4
1692 (math-ln-raw-3 (math-div-float (math-sub-float x
'(float 1 0))
1693 (math-add-float x
'(float 1 0)))))))
1695 (defun math-ln-raw-3 (x) ; [F F]
1696 (math-mul-float (math-ln-raw-series x
3 x
(math-sqr-float x
))
1699 ;;; Compute ln((1+x)/(1-x))
1700 (defun math-ln-raw-series (sum n x xsqr
)
1701 (math-working "log" sum
)
1702 (let* ((nextx (math-mul-float x xsqr
))
1703 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1704 (if (math-nearly-equal-float sum nextsum
)
1706 (math-ln-raw-series nextsum
(+ n
2) nextx xsqr
))))
1708 (defun math-ln-plus-1-raw (x)
1709 (math-lnp1-series x
2 x
(math-neg x
)))
1711 (defun math-lnp1-series (sum n xpow x
)
1712 (math-working "lnp1" sum
)
1713 (let* ((nextx (math-mul-float xpow x
))
1714 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1715 (if (math-nearly-equal-float sum nextsum
)
1717 (math-lnp1-series nextsum
(1+ n
) nextx x
))))
1719 (math-defcache math-ln-10
(float (bigpos 018 684 045 994 092 585 302 2) -
21)
1720 (math-ln-raw-2 '(float 1 1)))
1722 (math-defcache math-ln-2
(float (bigpos 417 309 945 559 180 147 693) -
21)
1723 (math-ln-raw-3 (math-float '(frac 1 3))))
1727 ;;; Hyperbolic functions.
1729 (defun calcFunc-sinh (x) ; [N N] [Public]
1731 (math-expand-formulas
1733 (list '/ (list '-
(list 'calcFunc-exp x
)
1734 (list 'calcFunc-exp
(list 'neg x
))) 2)))
1736 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1737 (math-with-extra-prec 2
1738 (let ((expx (math-exp-raw (math-float x
))))
1739 (math-mul (math-add expx
(math-div -
1 expx
)) '(float 5 -
1)))))
1740 ((eq (car-safe x
) 'sdev
)
1741 (math-make-sdev (calcFunc-sinh (nth 1 x
))
1742 (math-mul (nth 2 x
) (calcFunc-cosh (nth 1 x
)))))
1744 (math-sort-intv (nth 1 x
)
1745 (calcFunc-sinh (nth 2 x
))
1746 (calcFunc-sinh (nth 3 x
))))
1747 ((or (equal x
'(var inf var-inf
))
1748 (equal x
'(neg (var inf var-inf
)))
1749 (equal x
'(var nan var-nan
)))
1751 (t (calc-record-why 'numberp x
)
1752 (list 'calcFunc-sinh x
))))
1753 (put 'calcFunc-sinh
'math-expandable t
)
1755 (defun calcFunc-cosh (x) ; [N N] [Public]
1757 (math-expand-formulas
1759 (list '/ (list '+ (list 'calcFunc-exp x
)
1760 (list 'calcFunc-exp
(list 'neg x
))) 2)))
1762 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1763 (math-with-extra-prec 2
1764 (let ((expx (math-exp-raw (math-float x
))))
1765 (math-mul (math-add expx
(math-div 1 expx
)) '(float 5 -
1)))))
1766 ((eq (car-safe x
) 'sdev
)
1767 (math-make-sdev (calcFunc-cosh (nth 1 x
))
1769 (calcFunc-sinh (nth 1 x
)))))
1770 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1771 (setq x
(math-abs x
))
1772 (math-sort-intv (nth 1 x
)
1773 (calcFunc-cosh (nth 2 x
))
1774 (calcFunc-cosh (nth 3 x
))))
1775 ((or (equal x
'(var inf var-inf
))
1776 (equal x
'(neg (var inf var-inf
)))
1777 (equal x
'(var nan var-nan
)))
1779 (t (calc-record-why 'numberp x
)
1780 (list 'calcFunc-cosh x
))))
1781 (put 'calcFunc-cosh
'math-expandable t
)
1783 (defun calcFunc-tanh (x) ; [N N] [Public]
1785 (math-expand-formulas
1787 (let ((expx (list 'calcFunc-exp x
))
1788 (expmx (list 'calcFunc-exp
(list 'neg x
))))
1790 (list '/ (list '- expx expmx
) (list '+ expx expmx
))))))
1792 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1793 (math-with-extra-prec 2
1794 (let* ((expx (calcFunc-exp (math-float x
)))
1795 (expmx (math-div 1 expx
)))
1796 (math-div (math-sub expx expmx
)
1797 (math-add expx expmx
)))))
1798 ((eq (car-safe x
) 'sdev
)
1799 (math-make-sdev (calcFunc-tanh (nth 1 x
))
1801 (math-sqr (calcFunc-cosh (nth 1 x
))))))
1803 (math-sort-intv (nth 1 x
)
1804 (calcFunc-tanh (nth 2 x
))
1805 (calcFunc-tanh (nth 3 x
))))
1806 ((equal x
'(var inf var-inf
))
1808 ((equal x
'(neg (var inf var-inf
)))
1810 ((equal x
'(var nan var-nan
))
1812 (t (calc-record-why 'numberp x
)
1813 (list 'calcFunc-tanh x
))))
1814 (put 'calcFunc-tanh
'math-expandable t
)
1816 (defun calcFunc-sech (x) ; [N N] [Public]
1818 (math-expand-formulas
1820 (list '/ 2 (list '+ (list 'calcFunc-exp x
)
1821 (list 'calcFunc-exp
(list 'neg x
))))))
1823 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1824 (math-with-extra-prec 2
1825 (let ((expx (math-exp-raw (math-float x
))))
1826 (math-div '(float 2 0) (math-add expx
(math-div 1 expx
))))))
1827 ((eq (car-safe x
) 'sdev
)
1828 (math-make-sdev (calcFunc-sech (nth 1 x
))
1830 (math-mul (calcFunc-sech (nth 1 x
))
1831 (calcFunc-tanh (nth 1 x
))))))
1832 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1833 (setq x
(math-abs x
))
1834 (math-sort-intv (nth 1 x
)
1835 (calcFunc-sech (nth 2 x
))
1836 (calcFunc-sech (nth 3 x
))))
1837 ((or (equal x
'(var inf var-inf
))
1838 (equal x
'(neg (var inf var-inf
))))
1840 ((equal x
'(var nan var-nan
))
1842 (t (calc-record-why 'numberp x
)
1843 (list 'calcFunc-sech x
))))
1844 (put 'calcFunc-sech
'math-expandable t
)
1846 (defun calcFunc-csch (x) ; [N N] [Public]
1847 (cond ((eq x
0) (math-div 1 0))
1848 (math-expand-formulas
1850 (list '/ 2 (list '-
(list 'calcFunc-exp x
)
1851 (list 'calcFunc-exp
(list 'neg x
))))))
1853 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1854 (math-with-extra-prec 2
1855 (let ((expx (math-exp-raw (math-float x
))))
1856 (math-div '(float 2 0) (math-add expx
(math-div -
1 expx
))))))
1857 ((eq (car-safe x
) 'sdev
)
1858 (math-make-sdev (calcFunc-csch (nth 1 x
))
1860 (math-mul (calcFunc-csch (nth 1 x
))
1861 (calcFunc-coth (nth 1 x
))))))
1863 (if (and (Math-negp (nth 2 x
))
1864 (Math-posp (nth 3 x
)))
1865 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
1866 (math-sort-intv (nth 1 x
)
1867 (calcFunc-csch (nth 2 x
))
1868 (calcFunc-csch (nth 3 x
)))))
1869 ((or (equal x
'(var inf var-inf
))
1870 (equal x
'(neg (var inf var-inf
))))
1872 ((equal x
'(var nan var-nan
))
1874 (t (calc-record-why 'numberp x
)
1875 (list 'calcFunc-csch x
))))
1876 (put 'calcFunc-csch
'math-expandable t
)
1878 (defun calcFunc-coth (x) ; [N N] [Public]
1879 (cond ((eq x
0) (math-div 1 0))
1880 (math-expand-formulas
1882 (let ((expx (list 'calcFunc-exp x
))
1883 (expmx (list 'calcFunc-exp
(list 'neg x
))))
1885 (list '/ (list '+ expx expmx
) (list '- expx expmx
))))))
1887 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1888 (math-with-extra-prec 2
1889 (let* ((expx (calcFunc-exp (math-float x
)))
1890 (expmx (math-div 1 expx
)))
1891 (math-div (math-add expx expmx
)
1892 (math-sub expx expmx
)))))
1893 ((eq (car-safe x
) 'sdev
)
1894 (math-make-sdev (calcFunc-coth (nth 1 x
))
1896 (math-sqr (calcFunc-sinh (nth 1 x
))))))
1898 (if (and (Math-negp (nth 2 x
))
1899 (Math-posp (nth 3 x
)))
1900 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
1901 (math-sort-intv (nth 1 x
)
1902 (calcFunc-coth (nth 2 x
))
1903 (calcFunc-coth (nth 3 x
)))))
1904 ((equal x
'(var inf var-inf
))
1906 ((equal x
'(neg (var inf var-inf
)))
1908 ((equal x
'(var nan var-nan
))
1910 (t (calc-record-why 'numberp x
)
1911 (list 'calcFunc-coth x
))))
1912 (put 'calcFunc-coth
'math-expandable t
)
1914 (defun calcFunc-arcsinh (x) ; [N N] [Public]
1916 (math-expand-formulas
1918 (list 'calcFunc-ln
(list '+ x
(list 'calcFunc-sqrt
1919 (list '+ (list '^ x
2) 1))))))
1921 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1922 (math-with-extra-prec 2
1923 (math-ln-raw (math-add x
(math-sqrt-raw (math-add (math-sqr x
)
1925 ((eq (car-safe x
) 'sdev
)
1926 (math-make-sdev (calcFunc-arcsinh (nth 1 x
))
1929 (math-add (math-sqr (nth 1 x
)) 1)))))
1931 (math-sort-intv (nth 1 x
)
1932 (calcFunc-arcsinh (nth 2 x
))
1933 (calcFunc-arcsinh (nth 3 x
))))
1934 ((or (equal x
'(var inf var-inf
))
1935 (equal x
'(neg (var inf var-inf
)))
1936 (equal x
'(var nan var-nan
)))
1938 (t (calc-record-why 'numberp x
)
1939 (list 'calcFunc-arcsinh x
))))
1940 (put 'calcFunc-arcsinh
'math-expandable t
)
1942 (defun calcFunc-arccosh (x) ; [N N] [Public]
1944 ((and (eq x -
1) calc-symbolic-mode
)
1946 ((and (eq x
0) calc-symbolic-mode
)
1947 (math-div (math-mul '(var pi var-pi
) '(var i var-i
)) 2))
1948 (math-expand-formulas
1950 (list 'calcFunc-ln
(list '+ x
(list 'calcFunc-sqrt
1951 (list '-
(list '^ x
2) 1))))))
1953 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1954 (if (Math-equal-int x -
1)
1955 (math-imaginary (math-pi))
1956 (math-with-extra-prec 2
1957 (if (or t
; need to do this even in the real case!
1958 (memq (car-safe x
) '(cplx polar
)))
1959 (let ((xp1 (math-add 1 x
))) ; this gets the branch cuts right
1961 (math-add x
(math-mul xp1
1968 (math-add x
(math-sqrt-raw (math-add (math-sqr x
)
1969 '(float -
1 0)))))))))
1970 ((eq (car-safe x
) 'sdev
)
1971 (math-make-sdev (calcFunc-arccosh (nth 1 x
))
1974 (math-add (math-sqr (nth 1 x
)) -
1)))))
1976 (math-sort-intv (nth 1 x
)
1977 (calcFunc-arccosh (nth 2 x
))
1978 (calcFunc-arccosh (nth 3 x
))))
1979 ((or (equal x
'(var inf var-inf
))
1980 (equal x
'(neg (var inf var-inf
)))
1981 (equal x
'(var nan var-nan
)))
1983 (t (calc-record-why 'numberp x
)
1984 (list 'calcFunc-arccosh x
))))
1985 (put 'calcFunc-arccosh
'math-expandable t
)
1987 (defun calcFunc-arctanh (x) ; [N N] [Public]
1989 ((and (Math-equal-int x
1) calc-infinite-mode
)
1991 ((and (Math-equal-int x -
1) calc-infinite-mode
)
1992 '(neg (var inf var-inf
)))
1993 (math-expand-formulas
1995 (list 'calcFunc-ln
(list '+ 1 x
))
1996 (list 'calcFunc-ln
(list '-
1 x
))) 2))
1998 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1999 (math-with-extra-prec 2
2000 (if (or (memq (car-safe x
) '(cplx polar
))
2002 (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x
))
2003 (math-ln-raw (math-sub '(float 1 0) x
)))
2005 (if (and (math-equal-int x
1) calc-infinite-mode
)
2007 (if (and (math-equal-int x -
1) calc-infinite-mode
)
2008 '(neg (var inf var-inf
))
2009 (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x
)
2012 ((eq (car-safe x
) 'sdev
)
2013 (math-make-sdev (calcFunc-arctanh (nth 1 x
))
2015 (math-sub 1 (math-sqr (nth 1 x
))))))
2017 (math-sort-intv (nth 1 x
)
2018 (calcFunc-arctanh (nth 2 x
))
2019 (calcFunc-arctanh (nth 3 x
))))
2020 ((equal x
'(var nan var-nan
))
2022 (t (calc-record-why 'numberp x
)
2023 (list 'calcFunc-arctanh x
))))
2024 (put 'calcFunc-arctanh
'math-expandable t
)
2027 ;;; Convert A from HMS or degrees to radians.
2028 (defun calcFunc-rad (a) ; [R R] [Public]
2029 (cond ((or (Math-numberp a
)
2031 (math-with-extra-prec 2
2032 (math-mul a
(math-pi-over-180))))
2034 (math-from-hms a
'rad
))
2036 (math-make-sdev (calcFunc-rad (nth 1 a
))
2037 (calcFunc-rad (nth 2 a
))))
2038 (math-expand-formulas
2039 (math-div (math-mul a
'(var pi var-pi
)) 180))
2040 ((math-infinitep a
) a
)
2041 (t (list 'calcFunc-rad a
))))
2042 (put 'calcFunc-rad
'math-expandable t
)
2044 ;;; Convert A from HMS or radians to degrees.
2045 (defun calcFunc-deg (a) ; [R R] [Public]
2046 (cond ((or (Math-numberp a
)
2048 (math-with-extra-prec 2
2049 (math-div a
(math-pi-over-180))))
2051 (math-from-hms a
'deg
))
2053 (math-make-sdev (calcFunc-deg (nth 1 a
))
2054 (calcFunc-deg (nth 2 a
))))
2055 (math-expand-formulas
2056 (math-div (math-mul 180 a
) '(var pi var-pi
)))
2057 ((math-infinitep a
) a
)
2058 (t (list 'calcFunc-deg a
))))
2059 (put 'calcFunc-deg
'math-expandable t
)
2061 (provide 'calc-math
)
2063 ;;; arch-tag: c7367e8e-d0b8-4f70-8577-2fb3f31dbb4c
2064 ;;; calc-math.el ends here