1 ;;; calc-math.el --- mathematical functions for Calc
3 ;; Copyright (C) 1990-1993, 2001-2015 Free Software Foundation, Inc.
5 ;; Author: David Gillespie <daveg@synaptics.com>
6 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
8 ;; This file is part of GNU Emacs.
10 ;; GNU Emacs is free software: you can redistribute it and/or modify
11 ;; it under the terms of the GNU General Public License as published by
12 ;; the Free Software Foundation, either version 3 of the License, or
13 ;; (at your option) any later version.
15 ;; GNU Emacs is distributed in the hope that it will be useful,
16 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 ;; GNU General Public License for more details.
20 ;; You should have received a copy of the GNU General Public License
21 ;; along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
27 ;; This file is autoloaded from calc-ext.el.
33 ;;; Find out how many 9s in 9.9999... will give distinct Emacs floats,
34 ;;; then back off by one.
36 (defvar math-emacs-precision
39 (xx (+ x (* 9 (expt 10 (- n))))))
44 (setq xx (+ x (* 9 (expt 10 (- n)))))))
46 "The number of digits in an Emacs float.")
48 ;;; Find the largest power of 10 which is an Emacs float,
49 ;;; then back off by one so that any float d.dddd...eN
50 ;;; is an Emacs float, for acceptable d.dddd....
52 (defvar math-largest-emacs-expt
55 ;; The following loop is for efficiency; it should stop when
56 ;; 10^(2x) is too large. This could be indicated by a range
57 ;; error when computing 10^(2x) or an infinite value for 10^(2x).
62 (setq pow (condition-case nil
65 ;; The following loop should stop when 10^(x+1) is too large.
66 (setq pow (condition-case nil
73 (setq pow (condition-case nil
77 "The largest exponent which Calc will convert to an Emacs float.")
79 (defvar math-smallest-emacs-expt
81 (while (condition-case nil
86 (while (condition-case nil
91 "The smallest exponent which Calc will convert to an Emacs float.")
93 (defun math-use-emacs-fn (fn x)
94 "Use the native Emacs function FN to evaluate the Calc number X.
95 If this can't be done, return NIL."
97 (<= calc-internal-prec math-emacs-precision)
99 (let* ((fx (math-float x))
100 (xpon (+ (nth 2 x) (1- (math-numdigs (nth 1 x))))))
101 (and (<= math-smallest-emacs-expt xpon)
102 (<= xpon math-largest-emacs-expt)
109 ((calc-number-radix 10)
110 (calc-twos-complement-mode nil)
111 (calc-float-format (list 'float calc-internal-prec))
112 (calc-group-digits nil)
113 (calc-point-char "."))
114 (math-format-number (math-float x)))))))
117 (defun calc-sqrt (arg)
120 (if (calc-is-inverse)
121 (calc-unary-op "^2" 'calcFunc-sqr arg)
122 (calc-unary-op "sqrt" 'calcFunc-sqrt arg))))
124 (defun calc-isqrt (arg)
127 (if (calc-is-inverse)
128 (calc-unary-op "^2" 'calcFunc-sqr arg)
129 (calc-unary-op "isqt" 'calcFunc-isqrt arg))))
132 (defun calc-hypot (arg)
135 (calc-binary-op "hypt" 'calcFunc-hypot arg)))
142 (defun calc-log10 (arg)
144 (calc-hyperbolic-func)
147 (defun calc-log (arg)
150 (if (calc-is-inverse)
151 (calc-binary-op "alog" 'calcFunc-alog arg)
152 (calc-binary-op "log" 'calcFunc-log arg))))
154 (defun calc-ilog (arg)
157 (if (calc-is-inverse)
158 (calc-binary-op "alog" 'calcFunc-alog arg)
159 (calc-binary-op "ilog" 'calcFunc-ilog arg))))
161 (defun calc-lnp1 (arg)
166 (defun calc-exp (arg)
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)
175 (calc-unary-op "exp" 'calcFunc-exp arg)))))
177 (defun calc-expm1 (arg)
180 (if (calc-is-inverse)
181 (calc-unary-op "ln+1" 'calcFunc-lnp1 arg)
182 (calc-unary-op "ex-1" 'calcFunc-expm1 arg))))
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))
201 (calc-pop-push-record 0 "pi" (math-pi)))))))
203 (defun calc-sin (arg)
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)
212 (calc-unary-op "sin" 'calcFunc-sin arg)))))
214 (defun calc-arcsin (arg)
219 (defun calc-sinh (arg)
221 (calc-hyperbolic-func)
224 (defun calc-arcsinh (arg)
227 (calc-hyperbolic-func)
230 (defun calc-sec (arg)
233 (if (calc-is-hyperbolic)
234 (calc-unary-op "sech" 'calcFunc-sech arg)
235 (calc-unary-op "sec" 'calcFunc-sec arg))))
237 (defun calc-sech (arg)
239 (calc-hyperbolic-func)
242 (defun calc-cos (arg)
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)
251 (calc-unary-op "cos" 'calcFunc-cos arg)))))
253 (defun calc-arccos (arg)
258 (defun calc-cosh (arg)
260 (calc-hyperbolic-func)
263 (defun calc-arccosh (arg)
266 (calc-hyperbolic-func)
269 (defun calc-csc (arg)
272 (if (calc-is-hyperbolic)
273 (calc-unary-op "csch" 'calcFunc-csch arg)
274 (calc-unary-op "csc" 'calcFunc-csc arg))))
276 (defun calc-csch (arg)
278 (calc-hyperbolic-func)
281 (defun calc-sincos ()
284 (if (calc-is-inverse)
285 (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1)))
286 (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1))))))
288 (defun calc-tan (arg)
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)
297 (calc-unary-op "tan" 'calcFunc-tan arg)))))
299 (defun calc-arctan (arg)
304 (defun calc-tanh (arg)
306 (calc-hyperbolic-func)
309 (defun calc-arctanh (arg)
312 (calc-hyperbolic-func)
315 (defun calc-cot (arg)
318 (if (calc-is-hyperbolic)
319 (calc-unary-op "coth" 'calcFunc-coth arg)
320 (calc-unary-op "cot" 'calcFunc-cot arg))))
322 (defun calc-coth (arg)
324 (calc-hyperbolic-func)
327 (defun calc-arctan2 ()
330 (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2)))))
332 (defun calc-conj (arg)
335 (calc-unary-op "conj" 'calcFunc-conj arg)))
337 (defun calc-imaginary ()
340 (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))))
342 (defun calc-to-degrees (arg)
345 (calc-unary-op ">deg" 'calcFunc-deg arg)))
347 (defun calc-to-radians (arg)
350 (calc-unary-op ">rad" 'calcFunc-rad arg)))
353 (defun calc-degrees-mode (arg)
357 (calc-change-mode 'calc-angle-mode 'deg)
358 (message "Angles measured in degrees")))
359 ((= arg 2) (calc-radians-mode))
360 ((= arg 3) (calc-hms-mode))
361 (t (error "Prefix argument out of range"))))
363 (defun calc-radians-mode ()
366 (calc-change-mode 'calc-angle-mode 'rad)
367 (message "Angles measured in radians")))
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))
378 (math-isqrt-small a))
380 (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a))))))))
382 (defun calcFunc-isqrt (a)
384 (math-isqrt (math-floor a))
385 (math-floor (math-sqrt a))))
388 ;;; This returns (flag . result) where the flag is t if A is a perfect square.
389 (defun math-isqrt-bignum (a) ; [P.l L]
390 (let ((len (length a)))
392 (let* ((top (nthcdr (- len 2) a)))
393 (math-isqrt-bignum-iter
395 (math-scale-bignum-digit-size
397 (1+ (math-isqrt-small
398 (+ (* (nth 1 top) math-bignum-digit-size) (car top)))))
400 (let* ((top (nth (1- len) a)))
401 (math-isqrt-bignum-iter
403 (math-scale-bignum-digit-size
404 (list (1+ (math-isqrt-small top)))
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)))
414 (math-isqrt-bignum-iter a g2)
415 (cons (and (= comp 0)
416 (math-zerop-bignum (cdr q))
420 (defun math-zerop-bignum (a)
423 (while (eq (car (setq a (cdr a))) 0))
426 (defun math-scale-bignum-digit-size (a n) ; [L L S]
432 (defun math-isqrt-small (a) ; A > 0. [S S]
433 (let ((g (cond ((>= a 1000000) 10000)
438 (while (< (setq g2 (/ (+ g (/ a g)) 2)) g)
445 ;;; Compute the square root of a number.
446 ;;; [T N] if possible, else [F N] if possible, else [C N]. [Public]
449 (and (Math-zerop a) a)
450 (and (math-known-nonposp a)
451 (math-imaginary (math-sqrt (math-neg a))))
453 (let ((sqrt (math-isqrt-small a)))
454 (if (= (* sqrt sqrt) a)
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)))))
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)))
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))))))
489 (math-make-float (math-normalize
490 (cons 'bigpos (cdr res)))
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))
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)
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))))
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)))
531 (math-mul (math-sqrt (math-infinite-dir a inf)) inf)))
533 (calc-record-why 'numberp a)
534 (list 'calcFunc-sqrt a))))
535 (defalias 'calcFunc-sqrt 'math-sqrt)
537 (defun math-infinite-dir (a &optional inf)
538 (or inf (setq inf (math-infinitep a)))
539 (math-normalize (math-expr-subst a inf 1)))
541 (defun math-sqrt-float (a &optional guess) ; [F F F]
542 (if calc-symbolic-mode
543 (signal 'inexact-result nil)
544 (math-with-extra-prec 1 (math-sqrt-raw a guess))))
546 (defun math-sqrt-raw (a &optional guess) ; [F F F]
547 (if (not (Math-posp a))
550 ((math-use-emacs-fn 'sqrt a))
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)))))
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))
564 (if (math-nearly-equal-float g2 guess)
566 (math-sqrt-float-iter a g2))))
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)))
574 (and (not (consp ediff))
577 (= (math-numdigs (nth 1 a)) calc-internal-prec))))
579 (setq ediff (math-add (Math-integer-neg (nth 1 b))
580 (math-scale-int (nth 1 a) 1)))
581 (and (not (consp ediff))
584 (= (math-numdigs (nth 1 b)) calc-internal-prec)))
586 (setq ediff (math-add (Math-integer-neg (nth 1 a))
587 (math-scale-int (nth 1 b) 1)))
588 (and (not (consp ediff))
591 (= (math-numdigs (nth 1 a)) calc-internal-prec))))))
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)))
611 (math-nearly-equal-float a b))))
613 ;;; True if A is nearly zero compared to B. [P F F]
614 (defun math-nearly-zerop-float (a b)
616 (<= (+ (math-numdigs (nth 1 a)) (nth 2 a))
617 (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec)))))
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)
627 (math-nearly-zerop-float a b))))
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)
640 (calc-record-why 'scalarp a)
641 (list 'calcFunc-hypot a b)))
642 ((not (Math-scalarp b))
643 (if (math-infinitep 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))))
658 (defalias 'calcFunc-hypot 'math-hypot)
660 (defun calcFunc-sqr (x)
665 (defun math-nth-root (a n)
666 (cond ((= n 2) (math-sqrt a))
670 (let ((root (math-nth-root-integer a n)))
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)
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)))))
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.
699 (defvar math-nrf-nfm1)
701 (defun math-nth-root-float (a math-nrf-n &optional guess)
702 (math-inexact-result)
703 (math-with-extra-prec 1
704 (let ((math-nrf-nf (math-float math-nrf-n))
705 (math-nrf-nfm1 (math-float (1- math-nrf-n))))
706 (math-nth-root-float-iter a (or guess
708 1 (/ (+ (math-numdigs (nth 1 a))
713 (defun math-nth-root-float-iter (a guess)
714 (math-working "root" guess)
715 (let ((g2 (math-div-float (math-add-float (math-mul math-nrf-nfm1 guess)
717 a (math-ipow guess (1- math-nrf-n))))
719 (if (math-nearly-equal-float g2 guess)
721 (math-nth-root-float-iter a g2))))
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.
728 (defun math-nth-root-integer (a math-nri-n &optional guess) ; [I I S]
729 (math-nth-root-int-iter a (or guess
730 (math-scale-int 1 (/ (+ (math-numdigs a)
734 (defun math-nth-root-int-iter (a guess)
735 (math-working "root" guess)
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)))
739 (if (Math-natnum-lessp (car g2) guess)
740 (math-nth-root-int-iter a (car g2))
741 (cons (and (equal (car g2) guess)
746 (defun calcFunc-nroot (x n)
747 (calcFunc-pow x (if (integerp n)
754 ;;;; Transcendental functions.
756 ;;; All of these functions are defined on the complex plane.
757 ;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
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).
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)))
771 (defun math-from-radians (a) ; [N N]
772 (cond ((eq calc-angle-mode 'deg)
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))
780 (defun math-to-radians-2 (a &optional force-symbolic) ; [N N]
781 (cond ((eq (car-safe a) 'hms)
782 (math-from-hms a 'rad))
783 ((memq calc-angle-mode '(deg hms))
784 (if (or calc-symbolic-mode force-symbolic)
785 (math-div (math-mul a '(var pi var-pi)) 180)
786 (math-mul a (math-pi-over-180))))
789 (defun math-from-radians-2 (a &optional force-symbolic) ; [N N]
790 (cond ((memq calc-angle-mode '(deg hms))
791 (if (or calc-symbolic-mode force-symbolic)
792 (math-div (math-mul 180 a) '(var pi var-pi))
793 (math-div a (math-pi-over-180))))
798 ;;; Sine, cosine, and tangent.
800 (defun calcFunc-sin (x) ; [N N] [Public]
801 (cond ((and (integerp x)
802 (if (eq calc-angle-mode 'deg)
805 (aref [0 1 0 -1] (math-mod (/ x 90) 4)))
807 (math-with-extra-prec 2
808 (math-sin-raw (math-to-radians (math-float 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))
822 (t (calc-record-why 'scalarp x)
823 (list 'calcFunc-sin x))))
825 (defun calcFunc-cos (x) ; [N N] [Public]
826 (cond ((and (integerp x)
827 (if (eq calc-angle-mode 'deg)
830 (aref [1 0 -1 0] (math-mod (/ x 90) 4)))
832 (math-with-extra-prec 2
833 (math-cos-raw (math-to-radians (math-float 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)))))
855 (math-make-intv (logior (nth 1 x) 2)
858 (math-make-intv (logior (nth 1 x) 1)
862 (list 'intv 3 -1 1)))))
863 ((equal x '(var nan var-nan))
865 (t (calc-record-why 'scalarp x)
866 (list 'calcFunc-cos x))))
868 (defun calcFunc-sincos (x) ; [V N] [Public]
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]
873 (list 'vec (calcFunc-sin x) (calcFunc-cos x))))
875 (defun calcFunc-tan (x) ; [N N] [Public]
876 (cond ((and (integerp x)
877 (if (eq calc-angle-mode 'deg)
882 (math-with-extra-prec 2
883 (math-tan-raw (math-to-radians (math-float 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))
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))
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)
905 (nb (math-floor (math-div (math-sub (nth 3 xx)
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))
915 (t (calc-record-why 'scalarp x)
916 (list 'calcFunc-tan x))))
918 (defun calcFunc-sec (x)
919 (cond ((and (integerp x)
920 (eq calc-angle-mode 'deg)
926 (eq calc-angle-mode 'rad)
930 (math-with-extra-prec 2
931 (math-sec-raw (math-to-radians (math-float 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))
941 (calc-record-why "*Division by zero")
942 (list 'calcFunc-sec x))
943 (math-make-sdev (math-div-float '(float 1 0) (cdr sc))
945 (math-mul xs (car sc))
946 (math-sqr (cdr sc)))))))
947 (math-make-sdev (calcFunc-sec (nth 1 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)
959 (nb (math-floor (math-div (math-sub (nth 3 xx)
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)))))
971 (if (math-evenp (math-div (math-add naa 1) 2))
972 (math-make-intv (logior (nth 1 int) 2)
975 (math-make-intv (logior (nth 1 int) 1)
979 ((equal x '(var nan var-nan))
981 (t (calc-record-why 'scalarp x)
982 (list 'calcFunc-sec x))))
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)
992 (math-with-extra-prec 2
993 (math-csc-raw (math-to-radians (math-float 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))
1003 (calc-record-why "*Division by zero")
1004 (list 'calcFunc-csc x))
1005 (math-make-sdev (math-div-float '(float 1 0) (car sc))
1007 (math-mul xs (cdr sc))
1008 (math-sqr (car sc)))))))
1009 (math-make-sdev (calcFunc-csc (nth 1 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)))))
1029 (if (math-evenp (math-div naa 2))
1030 (math-make-intv (logior (nth 1 int) 2)
1033 (math-make-intv (logior (nth 1 int) 1)
1037 ((equal x '(var nan var-nan))
1039 (t (calc-record-why 'scalarp x)
1040 (list 'calcFunc-csc x))))
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)
1049 (math-with-extra-prec 2
1050 (math-cot-raw (math-to-radians (math-float 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))
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))
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))))
1070 (nb (math-floor (math-div (nth 3 xx) (math-pi)))))
1072 (math-sort-intv (nth 1 x)
1073 (math-cot-raw (nth 2 xx))
1074 (math-cot-raw (nth 3 xx))))))
1075 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
1076 ((equal x '(var nan var-nan))
1078 (t (calc-record-why 'scalarp x)
1079 (list 'calcFunc-cot x))))
1081 (defun math-sin-raw (x &optional orgx) ; [N N]
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))))
1087 (math-mul-float (car sc)
1088 (math-mul-float (math-add-float expx expmx)
1090 (math-mul-float (cdr sc)
1091 (math-mul-float (math-sub-float expx expmx)
1093 ((eq (car x) 'polar)
1094 (math-polar (math-sin-raw (math-complex x))))
1095 ((Math-integer-negp (nth 1 x))
1096 (math-neg-float (math-sin-raw (math-neg-float x) (if orgx orgx x))))
1097 ((math-lessp-float '(float 7 0) x) ; avoid inf loops due to roundoff
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)))))
1101 (defun math-cos-raw (x) ; [N N]
1102 (if (eq (car-safe x) 'polar)
1103 (math-polar (math-cos-raw (math-complex x)))
1104 (math-sin-raw (math-sub (math-pi-over-2) x) x)))
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)))
1115 (math-mul-float (math-sqr (car sc))
1117 (math-mul-float (math-sqr (cdr sc))
1119 (and (not (eq (nth 1 d) 0))
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))))
1126 (let ((cs (math-cos-raw x)))
1129 (math-div-float '(float 1 0) cs))))))
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)))
1140 (math-mul-float (math-sqr (car sc))
1142 (math-mul-float (math-sqr (cdr sc))
1144 (and (not (eq (nth 1 d) 0))
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)
1149 (math-polar (math-csc-raw (math-complex x))))
1151 (let ((sn (math-sin-raw x)))
1154 (math-div-float '(float 1 0) sn))))))
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)))
1167 (and (not (eq (nth 1 d) 0))
1170 (math-mul-float (car sc) (cdr sc))
1174 (math-mul-float sh ch)
1176 ((eq (car x) 'polar)
1177 (math-polar (math-cot-raw (math-complex x))))
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)))))))
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))
1189 (cons (math-sin-raw x) (math-cos-raw x)))
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)
1200 (and (not (eq (nth 1 d) 0))
1202 (math-div-float (car sc) d)
1203 (math-div-float (math-mul-float (math-sub-float expx
1205 '(float 5 -1)) d)))))
1206 ((eq (car x) 'polar)
1207 (math-polar (math-tan-raw (math-complex x))))
1209 (let ((sc (math-sin-cos-raw x)))
1210 (if (eq (nth 1 (cdr sc)) 0)
1211 (math-div (car sc) 0)
1212 (math-div-float (car sc) (cdr sc)))))))
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))
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)))
1223 ((math-with-extra-prec -1 (math-nearly-zerop-float x orgx))
1225 ((math-use-emacs-fn 'sin x))
1226 (calc-symbolic-mode (signal 'inexact-result nil))
1227 (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x)))))))
1229 (defun math-cos-raw-2 (x orgx) ; [F F]
1230 (cond ((math-with-extra-prec -1 (math-nearly-zerop-float x orgx))
1232 ((math-use-emacs-fn 'cos x))
1233 (calc-symbolic-mode (signal 'inexact-result nil))
1234 (t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
1236 (math-add-float '(float 1 0)
1237 (math-mul-float xnegsqr '(float 5 -1)))
1238 24 5 xnegsqr xnegsqr)))))
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)
1247 (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
1248 (+ n 2) nextx xnegsqr))))
1251 ;;; Inverse sine, cosine, tangent.
1253 (defun calcFunc-arcsin (x) ; [N N] [Public]
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))
1259 (math-with-extra-prec 2
1260 (math-from-radians (math-arcsin-raw (math-float x)))))
1262 (math-make-sdev (calcFunc-arcsin (nth 1 x))
1266 (math-sub 1 (math-sqr (nth 1 x))))))))
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))
1273 (t (calc-record-why 'numberp x)
1274 (list 'calcFunc-arcsin x))))
1276 (defun calcFunc-arccos (x) ; [N N] [Public]
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))
1282 (math-with-extra-prec 2
1283 (math-from-radians (math-arccos-raw (math-float x)))))
1285 (math-make-sdev (calcFunc-arccos (nth 1 x))
1289 (math-sub 1 (math-sqr (nth 1 x))))))))
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))
1296 (t (calc-record-why 'numberp x)
1297 (list 'calcFunc-arccos x))))
1299 (defun calcFunc-arctan (x) ; [N N] [Public]
1301 ((and (eq x 1) (eq calc-angle-mode 'deg)) 45)
1302 ((and (eq x -1) (eq calc-angle-mode 'deg)) -45)
1304 (math-with-extra-prec 2
1305 (math-from-radians (math-arctan-raw (math-float x)))))
1307 (math-make-sdev (calcFunc-arctan (nth 1 x))
1310 (math-add 1 (math-sqr (nth 1 x)))))))
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))
1321 (t (calc-record-why 'numberp x)
1322 (list 'calcFunc-arctan x))))
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))))
1331 (math-arctan2-raw x a))))
1333 (defun math-arccos-raw (x) ; [N N]
1334 (math-sub (math-pi-over-2) (math-arcsin-raw x)))
1336 (defun math-arctan-raw (x) ; [N N]
1337 (cond ((memq (car x) '(cplx polar))
1338 (math-with-extra-prec 2 ; extra-extra
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))))
1343 ((Math-integer-negp (nth 1 x))
1344 (math-neg-float (math-arctan-raw (math-neg-float x))))
1346 ((math-use-emacs-fn 'atan x))
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)
1359 (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x))))))
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)
1367 (math-arctan-series nextsum (+ n 2) nextx xnegsqr))))
1369 (defun calcFunc-arctan2 (y x) ; [F R R] [Public]
1372 (math-with-extra-prec 2
1373 (math-from-radians (math-arctan2-raw (math-float y)
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)))
1393 (calcFunc-arctan2 y x)
1394 '(var nan var-nan)))
1395 (calc-record-why 'anglep y)
1396 (list 'calcFunc-arctan2 y x))))
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)))
1405 (math-neg (math-pi-over-2))))
1407 (math-arctan-raw (math-div-float y x)))
1409 (math-add-float (math-arctan-raw (math-div-float y x))
1412 (math-sub-float (math-arctan-raw (math-div-float y x))
1415 (defun calcFunc-arcsincos (x) ; [V N] [Public]
1416 (if (and (Math-vectorp x)
1418 (calcFunc-arctan2 (nth 2 x) (nth 1 x))
1419 (math-reject-arg x "*Two-element vector expected")))
1423 ;;; Exponential function.
1425 (defun calcFunc-exp (x) ; [N N] [Public]
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))))
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))
1439 ((equal x '(neg (var inf var-inf)))
1441 ((equal x '(var nan var-nan))
1443 (t (calc-record-why 'numberp x)
1444 (list 'calcFunc-exp x))))
1446 (defun calcFunc-expm1 (x) ; [N N] [Public]
1448 ((math-zerop x) '(float 0 0))
1449 (calc-symbolic-mode (signal 'inexact-result nil))
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)
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))
1470 ((equal x '(neg (var inf var-inf)))
1472 ((equal x '(var nan var-nan))
1474 (t (calc-record-why 'numberp x)
1475 (list 'calcFunc-expm1 x))))
1477 (defun calcFunc-exp10 (x) ; [N N] [Public]
1480 (math-pow '(float 1 1) x)))
1482 (defun math-exp-raw (x) ; [N N]
1483 (cond ((math-zerop x) '(float 1 0))
1484 (calc-symbolic-mode (signal 'inexact-result nil))
1486 (let ((expx (math-exp-raw (nth 1 x)))
1487 (sc (math-sin-cos-raw (nth 2 x))))
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)))
1494 (math-exp-raw (nth 1 xc))
1495 (math-from-radians (nth 2 xc)))))
1496 ((math-use-emacs-fn 'exp x))
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)
1501 (if (math-lessp-float x '(float -921035 1))
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)
1507 (math-mul-float (math-ipow (math-sqrt-e) hint)
1508 (math-add-float '(float 1 0)
1509 (math-exp-minus-1-raw hfrac)))))
1510 (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x)))))
1512 (defun math-exp-minus-1-raw (x) ; [F F]
1513 (math-exp-series x 2 3 x x))
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)
1522 (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x))))
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")))
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))
1547 ((and (eq (car-safe x) '^)
1548 (equal (nth 1 x) '(var e var-e))
1549 (math-known-realp (nth 2 x)))
1552 (if (equal x '(var nan var-nan))
1554 '(var inf var-inf)))
1555 (t (calc-record-why 'numberp x)
1556 (list 'calcFunc-ln x))))
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)
1563 (let ((res (math-integer-log x 10)))
1565 (setq x (cdr res)))))
1567 ((and (eq (car-safe x) 'frac)
1569 (let ((res (math-integer-log (nth 2 x) 10)))
1571 (setq x (- (cdr res))))))
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))
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)))))
1588 (math-div-float (math-ln-raw-2
1589 (list 'float (nth 1 xf) (- xdigs)))
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))
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))))
1603 (if (equal x '(var nan var-nan))
1605 '(var inf var-inf)))
1606 (t (calc-record-why 'numberp x)
1607 (list 'calcFunc-log10 x))))
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)))
1615 (if calc-infinite-mode
1616 (math-div (calcFunc-ln x) (calcFunc-ln b))
1617 (math-reject-arg x "*Logarithm of zero")))
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)
1627 (if (math-floatp b) '(float 0 0) 0))
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)
1635 (bb (if (Math-lessp 1 b)
1637 (setq sign (- sign))
1639 (res (if (Math-lessp xx bb)
1640 (setq inv (math-integer-log bb xx))
1641 (math-integer-log xx bb))))
1644 (math-div 1 (* sign (cdr res)))
1645 (* sign (cdr res)))))))
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)
1654 (math-make-sdev (calcFunc-log (nth 1 x) b)
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)))
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)
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))
1672 (list 'calcFunc-log x b))))
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)))
1677 (t (math-pow b x))))
1679 (defun calcFunc-ilog (x b)
1680 (if (and (math-natnump x) (not (eq x 0))
1681 (math-natnump b) (not (eq b 0)))
1683 (math-reject-arg x "*Logarithm base one")
1684 (if (Math-natnum-lessp x b)
1686 (cdr (math-integer-log x b))))
1687 (math-floor (calcFunc-log x b))))
1689 (defun math-integer-log (x b)
1690 (let ((pows (list b))
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)))
1700 (while (and (setq pows (cdr pows))
1703 next (math-mul pow (car pows)))
1704 (or (Math-lessp x next)
1707 (cons (equal pow x) sum)))
1710 (defvar math-log-base-cache nil)
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)))))
1716 (nth 2 math-log-base-cache))
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")))
1724 ((math-zerop x) '(float 0 0))
1725 (calc-symbolic-mode (signal 'inexact-result nil))
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))))
1743 (if (equal x '(var nan var-nan))
1745 '(var inf var-inf)))
1746 (t (calc-record-why 'numberp x)
1747 (list 'calcFunc-lnp1 x))))
1749 (defun math-ln-raw (x) ; [N N] --- must be float format!
1750 (cond ((eq (car-safe x) '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))))
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)
1763 (calc-symbolic-mode (signal 'inexact-result nil))
1764 ((math-posp (nth 1 x)) ; positive and real
1766 ((math-use-emacs-fn 'log x))
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))
1773 (math-reject-arg x "*Logarithm of zero"))
1774 ((eq calc-complex-mode 'polar) ; negative and real
1776 (list 'cplx ; negative and real
1777 (math-ln-raw (math-neg-float x))
1779 (t (list 'cplx ; negative and real
1780 (math-ln-raw (math-neg-float x))
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)))
1787 (t ; now .7 < x <= 1.4
1788 (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
1789 (math-add-float x '(float 1 0)))))))
1791 (defun math-ln-raw-3 (x) ; [F F]
1792 (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
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)
1802 (math-ln-raw-series nextsum (+ n 2) nextx xsqr))))
1804 (defun math-ln-plus-1-raw (x)
1805 (math-lnp1-series x 2 x (math-neg x)))
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)
1813 (math-lnp1-series nextsum (1+ n) nextx x))))
1815 (defconst math-approx-ln-10
1816 (math-read-number-simple "2.302585092994045684018")
1817 "An approximation for ln(10).")
1819 (math-defcache math-ln-10 math-approx-ln-10
1820 (math-ln-raw-2 '(float 1 1)))
1822 (defconst math-approx-ln-2
1823 (math-read-number-simple "0.693147180559945309417")
1824 "An approximation for ln(2).")
1826 (math-defcache math-ln-2 math-approx-ln-2
1827 (math-ln-raw-3 (math-float '(frac 1 3))))
1831 ;;; Hyperbolic functions.
1833 (defun calcFunc-sinh (x) ; [N N] [Public]
1835 (math-expand-formulas
1837 (list '/ (list '- (list 'calcFunc-exp x)
1838 (list 'calcFunc-exp (list 'neg x))) 2)))
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)))))
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)))
1855 (t (calc-record-why 'numberp x)
1856 (list 'calcFunc-sinh x))))
1857 (put 'calcFunc-sinh 'math-expandable t)
1859 (defun calcFunc-cosh (x) ; [N N] [Public]
1861 (math-expand-formulas
1863 (list '/ (list '+ (list 'calcFunc-exp x)
1864 (list 'calcFunc-exp (list 'neg x))) 2)))
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))
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)))
1883 (t (calc-record-why 'numberp x)
1884 (list 'calcFunc-cosh x))))
1885 (put 'calcFunc-cosh 'math-expandable t)
1887 (defun calcFunc-tanh (x) ; [N N] [Public]
1889 (math-expand-formulas
1891 (let ((expx (list 'calcFunc-exp x))
1892 (expmx (list 'calcFunc-exp (list 'neg x))))
1894 (list '/ (list '- expx expmx) (list '+ expx expmx))))))
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))
1905 (math-sqr (calcFunc-cosh (nth 1 x))))))
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))
1912 ((equal x '(neg (var inf var-inf)))
1914 ((equal x '(var nan var-nan))
1916 (t (calc-record-why 'numberp x)
1917 (list 'calcFunc-tanh x))))
1918 (put 'calcFunc-tanh 'math-expandable t)
1920 (defun calcFunc-sech (x) ; [N N] [Public]
1922 (math-expand-formulas
1924 (list '/ 2 (list '+ (list 'calcFunc-exp x)
1925 (list 'calcFunc-exp (list 'neg 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))
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))))
1944 ((equal x '(var nan var-nan))
1946 (t (calc-record-why 'numberp x)
1947 (list 'calcFunc-sech x))))
1948 (put 'calcFunc-sech 'math-expandable t)
1950 (defun calcFunc-csch (x) ; [N N] [Public]
1951 (cond ((eq x 0) (math-div 1 0))
1952 (math-expand-formulas
1954 (list '/ 2 (list '- (list 'calcFunc-exp x)
1955 (list 'calcFunc-exp (list 'neg 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))
1964 (math-mul (calcFunc-csch (nth 1 x))
1965 (calcFunc-coth (nth 1 x))))))
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))))
1976 ((equal x '(var nan var-nan))
1978 (t (calc-record-why 'numberp x)
1979 (list 'calcFunc-csch x))))
1980 (put 'calcFunc-csch 'math-expandable t)
1982 (defun calcFunc-coth (x) ; [N N] [Public]
1983 (cond ((eq x 0) (math-div 1 0))
1984 (math-expand-formulas
1986 (let ((expx (list 'calcFunc-exp x))
1987 (expmx (list 'calcFunc-exp (list 'neg x))))
1989 (list '/ (list '+ expx expmx) (list '- expx expmx))))))
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))
2000 (math-sqr (calcFunc-sinh (nth 1 x))))))
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))
2010 ((equal x '(neg (var inf var-inf)))
2012 ((equal x '(var nan var-nan))
2014 (t (calc-record-why 'numberp x)
2015 (list 'calcFunc-coth x))))
2016 (put 'calcFunc-coth 'math-expandable t)
2018 (defun calcFunc-arcsinh (x) ; [N N] [Public]
2020 (math-expand-formulas
2022 (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
2023 (list '+ (list '^ x 2) 1))))))
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)
2029 ((eq (car-safe x) 'sdev)
2030 (math-make-sdev (calcFunc-arcsinh (nth 1 x))
2033 (math-add (math-sqr (nth 1 x)) 1)))))
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)))
2042 (t (calc-record-why 'numberp x)
2043 (list 'calcFunc-arcsinh x))))
2044 (put 'calcFunc-arcsinh 'math-expandable t)
2046 (defun calcFunc-arccosh (x) ; [N N] [Public]
2048 ((and (eq x -1) calc-symbolic-mode)
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
2054 (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
2055 (list '- (list '^ x 2) 1))))))
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
2065 (math-add x (math-mul xp1
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))
2078 (math-add (math-sqr (nth 1 x)) -1)))))
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)))
2087 (t (calc-record-why 'numberp x)
2088 (list 'calcFunc-arccosh x))))
2089 (put 'calcFunc-arccosh 'math-expandable t)
2091 (defun calcFunc-arctanh (x) ; [N N] [Public]
2093 ((and (Math-equal-int x 1) calc-infinite-mode)
2095 ((and (Math-equal-int x -1) calc-infinite-mode)
2096 '(neg (var inf var-inf)))
2097 (math-expand-formulas
2099 (list 'calcFunc-ln (list '+ 1 x))
2100 (list 'calcFunc-ln (list '- 1 x))) 2))
2102 (if calc-symbolic-mode (signal 'inexact-result nil))
2103 (math-with-extra-prec 2
2104 (if (or (memq (car-safe x) '(cplx polar))
2106 (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x))
2107 (math-ln-raw (math-sub '(float 1 0) x)))
2109 (if (and (math-equal-int x 1) calc-infinite-mode)
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)
2116 ((eq (car-safe x) 'sdev)
2117 (math-make-sdev (calcFunc-arctanh (nth 1 x))
2119 (math-sub 1 (math-sqr (nth 1 x))))))
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))
2126 (t (calc-record-why 'numberp x)
2127 (list 'calcFunc-arctanh x))))
2128 (put 'calcFunc-arctanh 'math-expandable t)
2131 ;;; Convert A from HMS or degrees to radians.
2132 (defun calcFunc-rad (a) ; [R R] [Public]
2133 (cond ((or (Math-numberp a)
2135 (math-with-extra-prec 2
2136 (math-mul a (math-pi-over-180))))
2138 (math-from-hms a 'rad))
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)
2145 (t (list 'calcFunc-rad a))))
2146 (put 'calcFunc-rad 'math-expandable t)
2148 ;;; Convert A from HMS or radians to degrees.
2149 (defun calcFunc-deg (a) ; [R R] [Public]
2150 (cond ((or (Math-numberp a)
2152 (math-with-extra-prec 2
2153 (math-div a (math-pi-over-180))))
2155 (math-from-hms a 'deg))
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)
2162 (t (list 'calcFunc-deg a))))
2163 (put 'calcFunc-deg 'math-expandable t)
2165 (provide 'calc-math)
2167 ;;; calc-math.el ends here