1 ;;; calc-math.el --- mathematical functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
4 ;; 2005, 2006, 2007 Free Software Foundation, Inc.
6 ;; Author: David Gillespie <daveg@synaptics.com>
7 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
9 ;; This file is part of GNU Emacs.
11 ;; GNU Emacs is free software; you can redistribute it and/or modify
12 ;; it under the terms of the GNU General Public License as published by
13 ;; the Free Software Foundation; either version 3, or (at your option)
16 ;; GNU Emacs is distributed in the hope that it will be useful,
17 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
18 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 ;; GNU General Public License for more details.
21 ;; You should have received a copy of the GNU General Public License
22 ;; along with GNU Emacs; see the file COPYING. If not, write to the
23 ;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
24 ;; Boston, MA 02110-1301, USA.
30 ;; This file is autoloaded from calc-ext.el.
36 ;;; Find out how many 9s in 9.9999... will give distinct Emacs floats,
37 ;;; then back off by one.
39 (defvar math-emacs-precision
42 (xx (+ x (* 9 (expt 10 (- n))))))
47 (setq xx (+ x (* 9 (expt 10 (- n)))))))
49 "The number of digits in an Emacs float.")
51 ;;; Find the largest power of 10 which is an Emacs float,
52 ;;; then back off by one so that any float d.dddd...eN
53 ;;; is an Emacs float, for acceptable d.dddd....
55 (defvar math-largest-emacs-expt
58 ;; The following loop is for efficiency; it should stop when
59 ;; 10^(2x) is too large. This could be indicated by a range
60 ;; error when computing 10^(2x) or an infinite value for 10^(2x).
65 (setq pow (condition-case nil
68 ;; The following loop should stop when 10^(x+1) is too large.
69 (setq pow (condition-case nil
76 (setq pow (condition-case nil
80 "The largest exponent which Calc will convert to an Emacs float.")
82 (defvar math-smallest-emacs-expt
84 (while (condition-case nil
89 (while (condition-case nil
94 "The smallest exponent which Calc will convert to an Emacs float.")
96 (defun math-use-emacs-fn (fn x)
97 "Use the native Emacs function FN to evaluate the Calc number X.
98 If this can't be done, return NIL."
100 (<= calc-internal-prec math-emacs-precision)
102 (let* ((fx (math-float x))
103 (xpon (+ (nth 2 x) (1- (math-numdigs (nth 1 x))))))
104 (and (<= math-smallest-emacs-expt xpon)
105 (<= xpon math-largest-emacs-expt)
110 (string-to-number (math-format-number (math-float x))))))
113 (defun calc-sqrt (arg)
116 (if (calc-is-inverse)
117 (calc-unary-op "^2" 'calcFunc-sqr arg)
118 (calc-unary-op "sqrt" 'calcFunc-sqrt arg))))
120 (defun calc-isqrt (arg)
123 (if (calc-is-inverse)
124 (calc-unary-op "^2" 'calcFunc-sqr arg)
125 (calc-unary-op "isqt" 'calcFunc-isqrt arg))))
128 (defun calc-hypot (arg)
131 (calc-binary-op "hypt" 'calcFunc-hypot arg)))
138 (defun calc-log10 (arg)
140 (calc-hyperbolic-func)
143 (defun calc-log (arg)
146 (if (calc-is-inverse)
147 (calc-binary-op "alog" 'calcFunc-alog arg)
148 (calc-binary-op "log" 'calcFunc-log arg))))
150 (defun calc-ilog (arg)
153 (if (calc-is-inverse)
154 (calc-binary-op "alog" 'calcFunc-alog arg)
155 (calc-binary-op "ilog" 'calcFunc-ilog arg))))
157 (defun calc-lnp1 (arg)
162 (defun calc-exp (arg)
165 (if (calc-is-hyperbolic)
166 (if (calc-is-inverse)
167 (calc-unary-op "lg10" 'calcFunc-log10 arg)
168 (calc-unary-op "10^" 'calcFunc-exp10 arg))
169 (if (calc-is-inverse)
170 (calc-unary-op "ln" 'calcFunc-ln arg)
171 (calc-unary-op "exp" 'calcFunc-exp arg)))))
173 (defun calc-expm1 (arg)
176 (if (calc-is-inverse)
177 (calc-unary-op "ln+1" 'calcFunc-lnp1 arg)
178 (calc-unary-op "ex-1" 'calcFunc-expm1 arg))))
183 (if (calc-is-inverse)
184 (if (calc-is-hyperbolic)
185 (if calc-symbolic-mode
186 (calc-pop-push-record 0 "phi" '(var phi var-phi))
187 (calc-pop-push-record 0 "phi" (math-phi)))
188 (if calc-symbolic-mode
189 (calc-pop-push-record 0 "gmma" '(var gamma var-gamma))
190 (calc-pop-push-record 0 "gmma" (math-gamma-const))))
191 (if (calc-is-hyperbolic)
192 (if calc-symbolic-mode
193 (calc-pop-push-record 0 "e" '(var e var-e))
194 (calc-pop-push-record 0 "e" (math-e)))
195 (if calc-symbolic-mode
196 (calc-pop-push-record 0 "pi" '(var pi var-pi))
197 (calc-pop-push-record 0 "pi" (math-pi)))))))
199 (defun calc-sin (arg)
202 (if (calc-is-hyperbolic)
203 (if (calc-is-inverse)
204 (calc-unary-op "asnh" 'calcFunc-arcsinh arg)
205 (calc-unary-op "sinh" 'calcFunc-sinh arg))
206 (if (calc-is-inverse)
207 (calc-unary-op "asin" 'calcFunc-arcsin arg)
208 (calc-unary-op "sin" 'calcFunc-sin arg)))))
210 (defun calc-arcsin (arg)
215 (defun calc-sinh (arg)
217 (calc-hyperbolic-func)
220 (defun calc-arcsinh (arg)
223 (calc-hyperbolic-func)
226 (defun calc-sec (arg)
229 (if (calc-is-hyperbolic)
230 (calc-unary-op "sech" 'calcFunc-sech arg)
231 (calc-unary-op "sec" 'calcFunc-sec arg))))
233 (defun calc-sech (arg)
235 (calc-hyperbolic-func)
238 (defun calc-cos (arg)
241 (if (calc-is-hyperbolic)
242 (if (calc-is-inverse)
243 (calc-unary-op "acsh" 'calcFunc-arccosh arg)
244 (calc-unary-op "cosh" 'calcFunc-cosh arg))
245 (if (calc-is-inverse)
246 (calc-unary-op "acos" 'calcFunc-arccos arg)
247 (calc-unary-op "cos" 'calcFunc-cos arg)))))
249 (defun calc-arccos (arg)
254 (defun calc-cosh (arg)
256 (calc-hyperbolic-func)
259 (defun calc-arccosh (arg)
262 (calc-hyperbolic-func)
265 (defun calc-csc (arg)
268 (if (calc-is-hyperbolic)
269 (calc-unary-op "csch" 'calcFunc-csch arg)
270 (calc-unary-op "csc" 'calcFunc-csc arg))))
272 (defun calc-csch (arg)
274 (calc-hyperbolic-func)
277 (defun calc-sincos ()
280 (if (calc-is-inverse)
281 (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1)))
282 (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1))))))
284 (defun calc-tan (arg)
287 (if (calc-is-hyperbolic)
288 (if (calc-is-inverse)
289 (calc-unary-op "atnh" 'calcFunc-arctanh arg)
290 (calc-unary-op "tanh" 'calcFunc-tanh arg))
291 (if (calc-is-inverse)
292 (calc-unary-op "atan" 'calcFunc-arctan arg)
293 (calc-unary-op "tan" 'calcFunc-tan arg)))))
295 (defun calc-arctan (arg)
300 (defun calc-tanh (arg)
302 (calc-hyperbolic-func)
305 (defun calc-arctanh (arg)
308 (calc-hyperbolic-func)
311 (defun calc-cot (arg)
314 (if (calc-is-hyperbolic)
315 (calc-unary-op "coth" 'calcFunc-coth arg)
316 (calc-unary-op "cot" 'calcFunc-cot arg))))
318 (defun calc-coth (arg)
320 (calc-hyperbolic-func)
323 (defun calc-arctan2 ()
326 (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2)))))
328 (defun calc-conj (arg)
331 (calc-unary-op "conj" 'calcFunc-conj arg)))
333 (defun calc-imaginary ()
336 (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))))
338 (defun calc-to-degrees (arg)
341 (calc-unary-op ">deg" 'calcFunc-deg arg)))
343 (defun calc-to-radians (arg)
346 (calc-unary-op ">rad" 'calcFunc-rad arg)))
349 (defun calc-degrees-mode (arg)
353 (calc-change-mode 'calc-angle-mode 'deg)
354 (message "Angles measured in degrees")))
355 ((= arg 2) (calc-radians-mode))
356 ((= arg 3) (calc-hms-mode))
357 (t (error "Prefix argument out of range"))))
359 (defun calc-radians-mode ()
362 (calc-change-mode 'calc-angle-mode 'rad)
363 (message "Angles measured in radians")))
366 ;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public]
367 ;;; This method takes advantage of the fact that Newton's method starting
368 ;;; with an overestimate always works, even using truncating integer division!
369 (defun math-isqrt (a)
370 (cond ((Math-zerop a) a)
371 ((not (math-natnump a))
372 (math-reject-arg a 'natnump))
374 (math-isqrt-small a))
376 (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a))))))))
378 (defun calcFunc-isqrt (a)
380 (math-isqrt (math-floor a))
381 (math-floor (math-sqrt a))))
384 ;;; This returns (flag . result) where the flag is t if A is a perfect square.
385 (defun math-isqrt-bignum (a) ; [P.l L]
386 (let ((len (length a)))
388 (let* ((top (nthcdr (- len 2) a)))
389 (math-isqrt-bignum-iter
391 (math-scale-bignum-digit-size
393 (1+ (math-isqrt-small
394 (+ (* (nth 1 top) math-bignum-digit-size) (car top)))))
396 (let* ((top (nth (1- len) a)))
397 (math-isqrt-bignum-iter
399 (math-scale-bignum-digit-size
400 (list (1+ (math-isqrt-small top)))
403 (defun math-isqrt-bignum-iter (a guess) ; [l L l]
404 (math-working "isqrt" (cons 'bigpos guess))
405 (let* ((q (math-div-bignum a guess))
406 (s (math-add-bignum (car q) guess))
407 (g2 (math-div2-bignum s))
408 (comp (math-compare-bignum g2 guess)))
410 (math-isqrt-bignum-iter a g2)
411 (cons (and (= comp 0)
412 (math-zerop-bignum (cdr q))
416 (defun math-zerop-bignum (a)
419 (while (eq (car (setq a (cdr a))) 0))
422 (defun math-scale-bignum-digit-size (a n) ; [L L S]
428 (defun math-isqrt-small (a) ; A > 0. [S S]
429 (let ((g (cond ((>= a 1000000) 10000)
434 (while (< (setq g2 (/ (+ g (/ a g)) 2)) g)
441 ;;; Compute the square root of a number.
442 ;;; [T N] if possible, else [F N] if possible, else [C N]. [Public]
445 (and (Math-zerop a) a)
446 (and (math-known-nonposp a)
447 (math-imaginary (math-sqrt (math-neg a))))
449 (let ((sqrt (math-isqrt-small a)))
450 (if (= (* sqrt sqrt) a)
452 (if calc-symbolic-mode
453 (list 'calcFunc-sqrt a)
454 (math-sqrt-float (math-float a) (math-float sqrt))))))
455 (and (eq (car-safe a) 'bigpos)
456 (let* ((res (math-isqrt-bignum (cdr a)))
457 (sqrt (math-normalize (cons 'bigpos (cdr res)))))
460 (if calc-symbolic-mode
461 (list 'calcFunc-sqrt a)
462 (math-sqrt-float (math-float a) (math-float sqrt))))))
463 (and (eq (car-safe a) 'frac)
464 (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a)))))
465 (num-sqrt (math-normalize (cons 'bigpos (cdr num-res))))
466 (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a)))))
467 (den-sqrt (math-normalize (cons 'bigpos (cdr den-res)))))
468 (if (and (car num-res) (car den-res))
469 (list 'frac num-sqrt den-sqrt)
470 (if calc-symbolic-mode
471 (if (or (car num-res) (car den-res))
472 (math-div (if (car num-res)
473 num-sqrt (list 'calcFunc-sqrt (nth 1 a)))
475 den-sqrt (list 'calcFunc-sqrt (nth 2 a))))
476 (list 'calcFunc-sqrt a))
477 (math-sqrt-float (math-float a)
478 (math-div (math-float num-sqrt) den-sqrt))))))
479 (and (eq (car-safe a) 'float)
480 (if calc-symbolic-mode
481 (if (= (% (nth 2 a) 2) 0)
482 (let ((res (math-isqrt-bignum
483 (cdr (Math-bignum-test (nth 1 a))))))
485 (math-make-float (math-normalize
486 (cons 'bigpos (cdr res)))
488 (signal 'inexact-result nil)))
489 (signal 'inexact-result nil))
490 (math-sqrt-float a)))
491 (and (eq (car-safe a) 'cplx)
492 (math-with-extra-prec 2
493 (let* ((d (math-abs a))
494 (imag (math-sqrt (math-mul (math-sub d (nth 1 a))
497 (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1)))
498 (if (math-negp (nth 2 a)) (math-neg imag) imag)))))
499 (and (eq (car-safe a) 'polar)
501 (math-sqrt (nth 1 a))
502 (math-mul (nth 2 a) '(float 5 -1))))
503 (and (eq (car-safe a) 'sdev)
504 (let ((sqrt (math-sqrt (nth 1 a))))
506 (math-div (nth 2 a) (math-mul sqrt 2)))))
507 (and (eq (car-safe a) 'intv)
508 (not (math-negp (nth 2 a)))
509 (math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a))))
510 (and (eq (car-safe a) '*)
511 (or (math-known-nonnegp (nth 1 a))
512 (math-known-nonnegp (nth 2 a)))
513 (math-mul (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
514 (and (eq (car-safe a) '/)
515 (or (and (math-known-nonnegp (nth 2 a))
516 (math-div (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
517 (and (math-known-nonnegp (nth 1 a))
518 (not (math-equal-int (nth 1 a) 1))
519 (math-mul (math-sqrt (nth 1 a))
520 (math-sqrt (math-div 1 (nth 2 a)))))))
521 (and (eq (car-safe a) '^)
522 (math-known-evenp (nth 2 a))
523 (math-known-realp (nth 1 a))
524 (math-abs (math-pow (nth 1 a) (math-div (nth 2 a) 2))))
525 (let ((inf (math-infinitep a)))
527 (math-mul (math-sqrt (math-infinite-dir a inf)) inf)))
529 (calc-record-why 'numberp a)
530 (list 'calcFunc-sqrt a))))
531 (defalias 'calcFunc-sqrt 'math-sqrt)
533 (defun math-infinite-dir (a &optional inf)
534 (or inf (setq inf (math-infinitep a)))
535 (math-normalize (math-expr-subst a inf 1)))
537 (defun math-sqrt-float (a &optional guess) ; [F F F]
538 (if calc-symbolic-mode
539 (signal 'inexact-result nil)
540 (math-with-extra-prec 1 (math-sqrt-raw a guess))))
542 (defun math-sqrt-raw (a &optional guess) ; [F F F]
543 (if (not (Math-posp a))
546 ((math-use-emacs-fn 'sqrt a))
549 (let ((ldiff (- (math-numdigs (nth 1 a)) 6)))
550 (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff)))
551 (setq guess (math-make-float (math-isqrt-small
552 (math-scale-int (nth 1 a) (- ldiff)))
553 (/ (+ (nth 2 a) ldiff) 2)))))
554 (math-sqrt-float-iter a guess)))))
556 (defun math-sqrt-float-iter (a guess) ; [F F F]
557 (math-working "sqrt" guess)
558 (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess))
560 (if (math-nearly-equal-float g2 guess)
562 (math-sqrt-float-iter a g2))))
564 ;;; True if A and B differ only in the last digit of precision. [P F F]
565 (defun math-nearly-equal-float (a b)
566 (let ((ediff (- (nth 2 a) (nth 2 b))))
567 (cond ((= ediff 0) ;; Expanded out for speed
568 (setq ediff (math-add (Math-integer-neg (nth 1 a)) (nth 1 b)))
570 (and (not (consp ediff))
573 (= (math-numdigs (nth 1 a)) calc-internal-prec))))
575 (setq ediff (math-add (Math-integer-neg (nth 1 b))
576 (math-scale-int (nth 1 a) 1)))
577 (and (not (consp ediff))
580 (= (math-numdigs (nth 1 b)) calc-internal-prec)))
582 (setq ediff (math-add (Math-integer-neg (nth 1 a))
583 (math-scale-int (nth 1 b) 1)))
584 (and (not (consp ediff))
587 (= (math-numdigs (nth 1 a)) calc-internal-prec))))))
589 (defun math-nearly-equal (a b) ; [P N N] [Public]
590 (setq a (math-float a))
591 (setq b (math-float b))
592 (if (eq (car a) 'polar) (setq a (math-complex a)))
593 (if (eq (car b) 'polar) (setq b (math-complex b)))
594 (if (eq (car a) 'cplx)
595 (if (eq (car b) 'cplx)
596 (and (or (math-nearly-equal-float (nth 1 a) (nth 1 b))
597 (and (math-nearly-zerop-float (nth 1 a) (nth 2 a))
598 (math-nearly-zerop-float (nth 1 b) (nth 2 b))))
599 (or (math-nearly-equal-float (nth 2 a) (nth 2 b))
600 (and (math-nearly-zerop-float (nth 2 a) (nth 1 a))
601 (math-nearly-zerop-float (nth 2 b) (nth 1 b)))))
602 (and (math-nearly-equal-float (nth 1 a) b)
603 (math-nearly-zerop-float (nth 2 a) b)))
604 (if (eq (car b) 'cplx)
605 (and (math-nearly-equal-float a (nth 1 b))
606 (math-nearly-zerop-float a (nth 2 b)))
607 (math-nearly-equal-float a b))))
609 ;;; True if A is nearly zero compared to B. [P F F]
610 (defun math-nearly-zerop-float (a b)
612 (<= (+ (math-numdigs (nth 1 a)) (nth 2 a))
613 (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec)))))
615 (defun math-nearly-zerop (a b) ; [P N R] [Public]
616 (setq a (math-float a))
617 (setq b (math-float b))
618 (if (eq (car a) 'cplx)
619 (and (math-nearly-zerop-float (nth 1 a) b)
620 (math-nearly-zerop-float (nth 2 a) b))
621 (if (eq (car a) 'polar)
622 (math-nearly-zerop-float (nth 1 a) b)
623 (math-nearly-zerop-float a b))))
625 ;;; This implementation could be improved, accuracy-wise.
626 (defun math-hypot (a b)
627 (cond ((Math-zerop a) (math-abs b))
628 ((Math-zerop b) (math-abs a))
629 ((not (Math-scalarp a))
630 (if (math-infinitep a)
631 (if (math-infinitep b)
636 (calc-record-why 'scalarp a)
637 (list 'calcFunc-hypot a b)))
638 ((not (Math-scalarp b))
639 (if (math-infinitep b)
641 (calc-record-why 'scalarp b)
642 (list 'calcFunc-hypot a b)))
643 ((and (Math-numberp a) (Math-numberp b))
644 (math-with-extra-prec 1
645 (math-sqrt (math-add (calcFunc-abssqr a) (calcFunc-abssqr b)))))
646 ((eq (car-safe a) 'hms)
647 (if (eq (car-safe b) 'hms) ; this helps sdev's of hms forms
648 (math-to-hms (math-hypot (math-from-hms a 'deg)
649 (math-from-hms b 'deg)))
650 (math-to-hms (math-hypot (math-from-hms a 'deg) b))))
651 ((eq (car-safe b) 'hms)
652 (math-to-hms (math-hypot a (math-from-hms b 'deg))))
654 (defalias 'calcFunc-hypot 'math-hypot)
656 (defun calcFunc-sqr (x)
661 (defun math-nth-root (a n)
662 (cond ((= n 2) (math-sqrt a))
666 (let ((root (math-nth-root-integer a n)))
669 (and (not calc-symbolic-mode)
670 (math-nth-root-float (math-float a) n
671 (math-float (cdr root)))))))
672 ((eq (car-safe a) 'frac)
673 (let* ((num-root (math-nth-root-integer (nth 1 a) n))
674 (den-root (math-nth-root-integer (nth 2 a) n)))
675 (if (and (car num-root) (car den-root))
676 (list 'frac (cdr num-root) (cdr den-root))
677 (and (not calc-symbolic-mode)
680 (math-div-float (math-float (cdr num-root))
681 (math-float (cdr den-root))))))))
682 ((eq (car-safe a) 'float)
683 (and (not calc-symbolic-mode)
684 (math-nth-root-float a n)))
685 ((eq (car-safe a) 'polar)
686 (let ((root (math-nth-root (nth 1 a) n)))
687 (and root (list 'polar root (math-div (nth 2 a) n)))))
690 ;; The variables math-nrf-n, math-nrf-nf and math-nrf-nfm1 are local
691 ;; to math-nth-root-float, but are used by math-nth-root-float-iter,
692 ;; which is called by math-nth-root-float.
695 (defvar math-nrf-nfm1)
697 (defun math-nth-root-float (a math-nrf-n &optional guess)
698 (math-inexact-result)
699 (math-with-extra-prec 1
700 (let ((math-nrf-nf (math-float math-nrf-n))
701 (math-nrf-nfm1 (math-float (1- math-nrf-n))))
702 (math-nth-root-float-iter a (or guess
704 1 (/ (+ (math-numdigs (nth 1 a))
709 (defun math-nth-root-float-iter (a guess)
710 (math-working "root" guess)
711 (let ((g2 (math-div-float (math-add-float (math-mul math-nrf-nfm1 guess)
713 a (math-ipow guess (1- math-nrf-n))))
715 (if (math-nearly-equal-float g2 guess)
717 (math-nth-root-float-iter a g2))))
719 ;; The variable math-nri-n is local to math-nth-root-integer, but
720 ;; is used by math-nth-root-int-iter, which is called by
721 ;; math-nth-root-int.
724 (defun math-nth-root-integer (a math-nri-n &optional guess) ; [I I S]
725 (math-nth-root-int-iter a (or guess
726 (math-scale-int 1 (/ (+ (math-numdigs a)
730 (defun math-nth-root-int-iter (a guess)
731 (math-working "root" guess)
732 (let* ((q (math-idivmod a (math-ipow guess (1- math-nri-n))))
733 (s (math-add (car q) (math-mul (1- math-nri-n) guess)))
734 (g2 (math-idivmod s math-nri-n)))
735 (if (Math-natnum-lessp (car g2) guess)
736 (math-nth-root-int-iter a (car g2))
737 (cons (and (equal (car g2) guess)
742 (defun calcFunc-nroot (x n)
743 (calcFunc-pow x (if (integerp n)
750 ;;;; Transcendental functions.
752 ;;; All of these functions are defined on the complex plane.
753 ;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
755 ;;; Most functions increase calc-internal-prec by 2 digits, then round
756 ;;; down afterward. "-raw" functions use the current precision, require
757 ;;; their arguments to be in float (or complex float) format, and always
758 ;;; work in radians (where applicable).
760 (defun math-to-radians (a) ; [N N]
761 (cond ((eq (car-safe a) 'hms)
762 (math-from-hms a 'rad))
763 ((memq calc-angle-mode '(deg hms))
764 (math-mul a (math-pi-over-180)))
767 (defun math-from-radians (a) ; [N N]
768 (cond ((eq calc-angle-mode 'deg)
770 (math-div a (math-pi-over-180))
771 (list 'calcFunc-deg a)))
772 ((eq calc-angle-mode 'hms)
773 (math-to-hms a 'rad))
776 (defun math-to-radians-2 (a) ; [N N]
777 (cond ((eq (car-safe a) 'hms)
778 (math-from-hms a 'rad))
779 ((memq calc-angle-mode '(deg hms))
780 (if calc-symbolic-mode
781 (math-div (math-mul a '(var pi var-pi)) 180)
782 (math-mul a (math-pi-over-180))))
785 (defun math-from-radians-2 (a) ; [N N]
786 (cond ((memq calc-angle-mode '(deg hms))
787 (if calc-symbolic-mode
788 (math-div (math-mul 180 a) '(var pi var-pi))
789 (math-div a (math-pi-over-180))))
794 ;;; Sine, cosine, and tangent.
796 (defun calcFunc-sin (x) ; [N N] [Public]
797 (cond ((and (integerp x)
798 (if (eq calc-angle-mode 'deg)
801 (aref [0 1 0 -1] (math-mod (/ x 90) 4)))
803 (math-with-extra-prec 2
804 (math-sin-raw (math-to-radians (math-float x)))))
807 (math-with-extra-prec 2
808 (let* ((xx (math-to-radians (math-float (nth 1 x))))
809 (xs (math-to-radians (math-float (nth 2 x))))
810 (sc (math-sin-cos-raw xx)))
811 (math-make-sdev (car sc) (math-mul xs (cdr sc)))))
812 (math-make-sdev (calcFunc-sin (nth 1 x))
813 (math-mul (nth 2 x) (calcFunc-cos (nth 1 x))))))
814 ((and (eq (car x) 'intv) (math-intv-constp x))
815 (calcFunc-cos (math-sub x (math-quarter-circle nil))))
816 ((equal x '(var nan var-nan))
818 (t (calc-record-why 'scalarp x)
819 (list 'calcFunc-sin x))))
821 (defun calcFunc-cos (x) ; [N N] [Public]
822 (cond ((and (integerp x)
823 (if (eq calc-angle-mode 'deg)
826 (aref [1 0 -1 0] (math-mod (/ x 90) 4)))
828 (math-with-extra-prec 2
829 (math-cos-raw (math-to-radians (math-float x)))))
832 (math-with-extra-prec 2
833 (let* ((xx (math-to-radians (math-float (nth 1 x))))
834 (xs (math-to-radians (math-float (nth 2 x))))
835 (sc (math-sin-cos-raw xx)))
836 (math-make-sdev (cdr sc) (math-mul xs (car sc)))))
837 (math-make-sdev (calcFunc-cos (nth 1 x))
838 (math-mul (nth 2 x) (calcFunc-sin (nth 1 x))))))
839 ((and (eq (car x) 'intv) (math-intv-constp x))
840 (math-with-extra-prec 2
841 (let* ((xx (math-to-radians (math-float x)))
842 (na (math-floor (math-div (nth 2 xx) (math-pi))))
843 (nb (math-floor (math-div (nth 3 xx) (math-pi))))
844 (span (math-sub nb na)))
845 (if (memq span '(0 1))
846 (let ((int (math-sort-intv (nth 1 x)
847 (math-cos-raw (nth 2 xx))
848 (math-cos-raw (nth 3 xx)))))
851 (math-make-intv (logior (nth 1 x) 2)
854 (math-make-intv (logior (nth 1 x) 1)
858 (list 'intv 3 -1 1)))))
859 ((equal x '(var nan var-nan))
861 (t (calc-record-why 'scalarp x)
862 (list 'calcFunc-cos x))))
864 (defun calcFunc-sincos (x) ; [V N] [Public]
866 (math-with-extra-prec 2
867 (let ((sc (math-sin-cos-raw (math-to-radians (math-float x)))))
868 (list 'vec (cdr sc) (car sc)))) ; the vector [cos, sin]
869 (list 'vec (calcFunc-sin x) (calcFunc-cos x))))
871 (defun calcFunc-tan (x) ; [N N] [Public]
872 (cond ((and (integerp x)
873 (if (eq calc-angle-mode 'deg)
878 (math-with-extra-prec 2
879 (math-tan-raw (math-to-radians (math-float x)))))
882 (math-with-extra-prec 2
883 (let* ((xx (math-to-radians (math-float (nth 1 x))))
884 (xs (math-to-radians (math-float (nth 2 x))))
885 (sc (math-sin-cos-raw xx)))
886 (if (and (math-zerop (cdr sc)) (not calc-infinite-mode))
888 (calc-record-why "*Division by zero")
889 (list 'calcFunc-tan x))
890 (math-make-sdev (math-div-float (car sc) (cdr sc))
891 (math-div-float xs (math-sqr (cdr sc)))))))
892 (math-make-sdev (calcFunc-tan (nth 1 x))
894 (math-sqr (calcFunc-cos (nth 1 x)))))))
895 ((and (eq (car x) 'intv) (math-intv-constp x))
896 (or (math-with-extra-prec 2
897 (let* ((xx (math-to-radians (math-float x)))
898 (na (math-floor (math-div (math-sub (nth 2 xx)
901 (nb (math-floor (math-div (math-sub (nth 3 xx)
905 (math-sort-intv (nth 1 x)
906 (math-tan-raw (nth 2 xx))
907 (math-tan-raw (nth 3 xx))))))
908 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
909 ((equal x '(var nan var-nan))
911 (t (calc-record-why 'scalarp x)
912 (list 'calcFunc-tan x))))
914 (defun calcFunc-sec (x)
915 (cond ((and (integerp x)
916 (eq calc-angle-mode 'deg)
922 (eq calc-angle-mode 'rad)
926 (math-with-extra-prec 2
927 (math-sec-raw (math-to-radians (math-float x)))))
930 (math-with-extra-prec 2
931 (let* ((xx (math-to-radians (math-float (nth 1 x))))
932 (xs (math-to-radians (math-float (nth 2 x))))
933 (sc (math-sin-cos-raw xx)))
934 (if (and (math-zerop (cdr sc))
935 (not calc-infinite-mode))
937 (calc-record-why "*Division by zero")
938 (list 'calcFunc-sec x))
939 (math-make-sdev (math-div-float '(float 1 0) (cdr sc))
941 (math-mul xs (car sc))
942 (math-sqr (cdr sc)))))))
943 (math-make-sdev (calcFunc-sec (nth 1 x))
946 (calcFunc-sin (nth 1 x)))
947 (math-sqr (calcFunc-cos (nth 1 x)))))))
948 ((and (eq (car x) 'intv)
949 (math-intv-constp x))
950 (math-with-extra-prec 2
951 (let* ((xx (math-to-radians (math-float x)))
952 (na (math-floor (math-div (math-sub (nth 2 xx)
955 (nb (math-floor (math-div (math-sub (nth 3 xx)
958 (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2))))
959 (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2))))
960 (span (math-sub nbb naa)))
961 (if (not (equal na nb))
962 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
963 (let ((int (math-sort-intv (nth 1 x)
964 (math-sec-raw (nth 2 xx))
965 (math-sec-raw (nth 3 xx)))))
967 (if (math-evenp (math-div (math-add naa 1) 2))
968 (math-make-intv (logior (nth 1 int) 2)
971 (math-make-intv (logior (nth 1 int) 1)
975 ((equal x '(var nan var-nan))
977 (t (calc-record-why 'scalarp x)
978 (list 'calcFunc-sec x))))
980 (defun calcFunc-csc (x)
981 (cond ((and (integerp x)
982 (eq calc-angle-mode 'deg)
983 (= (% (- x 90) 180) 0))
984 (if (= (% (- x 90) 360) 0)
988 (math-with-extra-prec 2
989 (math-csc-raw (math-to-radians (math-float x)))))
992 (math-with-extra-prec 2
993 (let* ((xx (math-to-radians (math-float (nth 1 x))))
994 (xs (math-to-radians (math-float (nth 2 x))))
995 (sc (math-sin-cos-raw xx)))
996 (if (and (math-zerop (car sc))
997 (not calc-infinite-mode))
999 (calc-record-why "*Division by zero")
1000 (list 'calcFunc-csc x))
1001 (math-make-sdev (math-div-float '(float 1 0) (car sc))
1003 (math-mul xs (cdr sc))
1004 (math-sqr (car sc)))))))
1005 (math-make-sdev (calcFunc-csc (nth 1 x))
1008 (calcFunc-cos (nth 1 x)))
1009 (math-sqr (calcFunc-sin (nth 1 x)))))))
1010 ((and (eq (car x) 'intv)
1011 (math-intv-constp x))
1012 (math-with-extra-prec 2
1013 (let* ((xx (math-to-radians (math-float x)))
1014 (na (math-floor (math-div (nth 2 xx) (math-pi))))
1015 (nb (math-floor (math-div (nth 3 xx) (math-pi))))
1016 (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2))))
1017 (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2))))
1018 (span (math-sub nbb naa)))
1019 (if (not (equal na nb))
1020 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1021 (let ((int (math-sort-intv (nth 1 x)
1022 (math-csc-raw (nth 2 xx))
1023 (math-csc-raw (nth 3 xx)))))
1025 (if (math-evenp (math-div naa 2))
1026 (math-make-intv (logior (nth 1 int) 2)
1029 (math-make-intv (logior (nth 1 int) 1)
1033 ((equal x '(var nan var-nan))
1035 (t (calc-record-why 'scalarp x)
1036 (list 'calcFunc-csc x))))
1038 (defun calcFunc-cot (x) ; [N N] [Public]
1039 (cond ((and (integerp x)
1040 (if (eq calc-angle-mode 'deg)
1041 (= (% (- x 90) 180) 0)
1045 (math-with-extra-prec 2
1046 (math-cot-raw (math-to-radians (math-float x)))))
1049 (math-with-extra-prec 2
1050 (let* ((xx (math-to-radians (math-float (nth 1 x))))
1051 (xs (math-to-radians (math-float (nth 2 x))))
1052 (sc (math-sin-cos-raw xx)))
1053 (if (and (math-zerop (car sc)) (not calc-infinite-mode))
1055 (calc-record-why "*Division by zero")
1056 (list 'calcFunc-cot x))
1057 (math-make-sdev (math-div-float (cdr sc) (car sc))
1058 (math-div-float xs (math-sqr (car sc)))))))
1059 (math-make-sdev (calcFunc-cot (nth 1 x))
1061 (math-sqr (calcFunc-sin (nth 1 x)))))))
1062 ((and (eq (car x) 'intv) (math-intv-constp x))
1063 (or (math-with-extra-prec 2
1064 (let* ((xx (math-to-radians (math-float x)))
1065 (na (math-floor (math-div (nth 2 xx) (math-pi))))
1066 (nb (math-floor (math-div (nth 3 xx) (math-pi)))))
1068 (math-sort-intv (nth 1 x)
1069 (math-cot-raw (nth 2 xx))
1070 (math-cot-raw (nth 3 xx))))))
1071 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
1072 ((equal x '(var nan var-nan))
1074 (t (calc-record-why 'scalarp x)
1075 (list 'calcFunc-cot x))))
1077 (defun math-sin-raw (x) ; [N N]
1078 (cond ((eq (car x) 'cplx)
1079 (let* ((expx (math-exp-raw (nth 2 x)))
1080 (expmx (math-div-float '(float 1 0) expx))
1081 (sc (math-sin-cos-raw (nth 1 x))))
1083 (math-mul-float (car sc)
1084 (math-mul-float (math-add-float expx expmx)
1086 (math-mul-float (cdr sc)
1087 (math-mul-float (math-sub-float expx expmx)
1089 ((eq (car x) 'polar)
1090 (math-polar (math-sin-raw (math-complex x))))
1091 ((Math-integer-negp (nth 1 x))
1092 (math-neg-float (math-sin-raw (math-neg-float x))))
1093 ((math-lessp-float '(float 7 0) x) ; avoid inf loops due to roundoff
1094 (math-sin-raw (math-mod x (math-two-pi))))
1095 (t (math-sin-raw-2 x x))))
1097 (defun math-cos-raw (x) ; [N N]
1098 (if (eq (car-safe x) 'polar)
1099 (math-polar (math-cos-raw (math-complex x)))
1100 (math-sin-raw (math-sub (math-pi-over-2) x))))
1102 (defun math-sec-raw (x) ; [N N]
1103 (cond ((eq (car x) 'cplx)
1104 (let* ((x (math-mul x '(float 1 0)))
1105 (expx (math-exp-raw (nth 2 x)))
1106 (expmx (math-div-float '(float 1 0) expx))
1107 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
1108 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
1109 (sc (math-sin-cos-raw (nth 1 x)))
1111 (math-mul-float (math-sqr (car sc))
1113 (math-mul-float (math-sqr (cdr sc))
1115 (and (not (eq (nth 1 d) 0))
1117 (math-div-float (math-mul-float (cdr sc) ch) d)
1118 (math-div-float (math-mul-float (car sc) sh) d)))))
1119 ((eq (car x) 'polar)
1120 (math-polar (math-sec-raw (math-complex x))))
1122 (let ((cs (math-cos-raw x)))
1125 (math-div-float '(float 1 0) cs))))))
1127 (defun math-csc-raw (x) ; [N N]
1128 (cond ((eq (car x) 'cplx)
1129 (let* ((x (math-mul x '(float 1 0)))
1130 (expx (math-exp-raw (nth 2 x)))
1131 (expmx (math-div-float '(float 1 0) expx))
1132 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
1133 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
1134 (sc (math-sin-cos-raw (nth 1 x)))
1136 (math-mul-float (math-sqr (car sc))
1138 (math-mul-float (math-sqr (cdr sc))
1140 (and (not (eq (nth 1 d) 0))
1142 (math-div-float (math-mul-float (car sc) ch) d)
1143 (math-div-float (math-mul-float (cdr sc) sh) d)))))
1144 ((eq (car x) 'polar)
1145 (math-polar (math-csc-raw (math-complex x))))
1147 (let ((sn (math-sin-raw x)))
1150 (math-div-float '(float 1 0) sn))))))
1152 (defun math-cot-raw (x) ; [N N]
1153 (cond ((eq (car x) 'cplx)
1154 (let* ((x (math-mul x '(float 1 0)))
1155 (expx (math-exp-raw (nth 2 x)))
1156 (expmx (math-div-float '(float 1 0) expx))
1157 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
1158 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
1159 (sc (math-sin-cos-raw (nth 1 x)))
1163 (and (not (eq (nth 1 d) 0))
1166 (math-mul-float (car sc) (cdr sc))
1170 (math-mul-float sh ch)
1172 ((eq (car x) 'polar)
1173 (math-polar (math-cot-raw (math-complex x))))
1175 (let ((sc (math-sin-cos-raw x)))
1176 (if (eq (nth 1 (car sc)) 0)
1177 (math-div (cdr sc) 0)
1178 (math-div-float (cdr sc) (car sc)))))))
1181 ;;; This could use a smarter method: Reduce x as in math-sin-raw, then
1182 ;;; compute either sin(x) or cos(x), whichever is smaller, and compute
1183 ;;; the other using the identity sin(x)^2 + cos(x)^2 = 1.
1184 (defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x))
1185 (cons (math-sin-raw x) (math-cos-raw x)))
1187 (defun math-tan-raw (x) ; [N N]
1188 (cond ((eq (car x) 'cplx)
1189 (let* ((x (math-mul x '(float 2 0)))
1190 (expx (math-exp-raw (nth 2 x)))
1191 (expmx (math-div-float '(float 1 0) expx))
1192 (sc (math-sin-cos-raw (nth 1 x)))
1193 (d (math-add-float (cdr sc)
1194 (math-mul-float (math-add-float expx expmx)
1196 (and (not (eq (nth 1 d) 0))
1198 (math-div-float (car sc) d)
1199 (math-div-float (math-mul-float (math-sub-float expx
1201 '(float 5 -1)) d)))))
1202 ((eq (car x) 'polar)
1203 (math-polar (math-tan-raw (math-complex x))))
1205 (let ((sc (math-sin-cos-raw x)))
1206 (if (eq (nth 1 (cdr sc)) 0)
1207 (math-div (car sc) 0)
1208 (math-div-float (car sc) (cdr sc)))))))
1210 (defun math-sin-raw-2 (x orgx) ; This avoids poss of inf recursion. [F F]
1211 (let ((xmpo2 (math-sub-float (math-pi-over-2) x)))
1212 (cond ((Math-integer-negp (nth 1 xmpo2))
1213 (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi))
1215 ((math-lessp-float (math-pi-over-4) x)
1216 (math-cos-raw-2 xmpo2 orgx))
1217 ((math-lessp-float x (math-neg (math-pi-over-4)))
1218 (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx)))
1219 ((math-nearly-zerop-float x orgx) '(float 0 0))
1220 ((math-use-emacs-fn 'sin x))
1221 (calc-symbolic-mode (signal 'inexact-result nil))
1222 (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x)))))))
1224 (defun math-cos-raw-2 (x orgx) ; [F F]
1225 (cond ((math-nearly-zerop-float x orgx) '(float 1 0))
1226 ((math-use-emacs-fn 'cos x))
1227 (calc-symbolic-mode (signal 'inexact-result nil))
1228 (t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
1230 (math-add-float '(float 1 0)
1231 (math-mul-float xnegsqr '(float 5 -1)))
1232 24 5 xnegsqr xnegsqr)))))
1234 (defun math-sin-series (sum nfac n x xnegsqr)
1235 (math-working "sin" sum)
1236 (let* ((nextx (math-mul-float x xnegsqr))
1237 (nextsum (math-add-float sum (math-div-float nextx
1238 (math-float nfac)))))
1239 (if (math-nearly-equal-float sum nextsum)
1241 (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
1242 (+ n 2) nextx xnegsqr))))
1245 ;;; Inverse sine, cosine, tangent.
1247 (defun calcFunc-arcsin (x) ; [N N] [Public]
1249 ((and (eq x 1) (eq calc-angle-mode 'deg)) 90)
1250 ((and (eq x -1) (eq calc-angle-mode 'deg)) -90)
1251 (calc-symbolic-mode (signal 'inexact-result nil))
1253 (math-with-extra-prec 2
1254 (math-from-radians (math-arcsin-raw (math-float x)))))
1256 (math-make-sdev (calcFunc-arcsin (nth 1 x))
1260 (math-sub 1 (math-sqr (nth 1 x))))))))
1262 (math-sort-intv (nth 1 x)
1263 (calcFunc-arcsin (nth 2 x))
1264 (calcFunc-arcsin (nth 3 x))))
1265 ((equal x '(var nan var-nan))
1267 (t (calc-record-why 'numberp x)
1268 (list 'calcFunc-arcsin x))))
1270 (defun calcFunc-arccos (x) ; [N N] [Public]
1272 ((and (eq x 0) (eq calc-angle-mode 'deg)) 90)
1273 ((and (eq x -1) (eq calc-angle-mode 'deg)) 180)
1274 (calc-symbolic-mode (signal 'inexact-result nil))
1276 (math-with-extra-prec 2
1277 (math-from-radians (math-arccos-raw (math-float x)))))
1279 (math-make-sdev (calcFunc-arccos (nth 1 x))
1283 (math-sub 1 (math-sqr (nth 1 x))))))))
1285 (math-sort-intv (nth 1 x)
1286 (calcFunc-arccos (nth 2 x))
1287 (calcFunc-arccos (nth 3 x))))
1288 ((equal x '(var nan var-nan))
1290 (t (calc-record-why 'numberp x)
1291 (list 'calcFunc-arccos x))))
1293 (defun calcFunc-arctan (x) ; [N N] [Public]
1295 ((and (eq x 1) (eq calc-angle-mode 'deg)) 45)
1296 ((and (eq x -1) (eq calc-angle-mode 'deg)) -45)
1298 (math-with-extra-prec 2
1299 (math-from-radians (math-arctan-raw (math-float x)))))
1301 (math-make-sdev (calcFunc-arctan (nth 1 x))
1304 (math-add 1 (math-sqr (nth 1 x)))))))
1306 (math-sort-intv (nth 1 x)
1307 (calcFunc-arctan (nth 2 x))
1308 (calcFunc-arctan (nth 3 x))))
1309 ((equal x '(var inf var-inf))
1310 (math-quarter-circle t))
1311 ((equal x '(neg (var inf var-inf)))
1312 (math-neg (math-quarter-circle t)))
1313 ((equal x '(var nan var-nan))
1315 (t (calc-record-why 'numberp x)
1316 (list 'calcFunc-arctan x))))
1318 (defun math-arcsin-raw (x) ; [N N]
1319 (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x)))))
1320 (if (or (memq (car x) '(cplx polar))
1321 (memq (car a) '(cplx polar)))
1322 (math-with-extra-prec 2 ; use extra precision for difficult case
1323 (math-mul '(cplx 0 -1)
1324 (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a))))
1325 (math-arctan2-raw x a))))
1327 (defun math-arccos-raw (x) ; [N N]
1328 (math-sub (math-pi-over-2) (math-arcsin-raw x)))
1330 (defun math-arctan-raw (x) ; [N N]
1331 (cond ((memq (car x) '(cplx polar))
1332 (math-with-extra-prec 2 ; extra-extra
1334 (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x)))
1335 (math-ln-raw (math-add 1 (math-mul '(cplx 0 -1) x))))
1337 ((Math-integer-negp (nth 1 x))
1338 (math-neg-float (math-arctan-raw (math-neg-float x))))
1340 ((math-use-emacs-fn 'atan x))
1341 (calc-symbolic-mode (signal 'inexact-result nil))
1342 ((math-equal-int x 1) (math-pi-over-4))
1343 ((math-equal-int x -1) (math-neg (math-pi-over-4)))
1344 ((math-lessp-float '(float 414214 -6) x) ; if x > sqrt(2) - 1, reduce
1345 (if (math-lessp-float '(float 1 0) x)
1346 (math-sub-float (math-mul-float (math-pi) '(float 5 -1))
1347 (math-arctan-raw (math-div-float '(float 1 0) x)))
1348 (math-sub-float (math-mul-float (math-pi) '(float 25 -2))
1349 (math-arctan-raw (math-div-float
1350 (math-sub-float '(float 1 0) x)
1351 (math-add-float '(float 1 0)
1353 (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x))))))
1355 (defun math-arctan-series (sum n x xnegsqr)
1356 (math-working "arctan" sum)
1357 (let* ((nextx (math-mul-float x xnegsqr))
1358 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1359 (if (math-nearly-equal-float sum nextsum)
1361 (math-arctan-series nextsum (+ n 2) nextx xnegsqr))))
1363 (defun calcFunc-arctan2 (y x) ; [F R R] [Public]
1366 (math-with-extra-prec 2
1367 (math-from-radians (math-arctan2-raw (math-float y)
1369 (calc-record-why 'anglep x)
1370 (list 'calcFunc-arctan2 y x))
1371 (if (and (or (math-infinitep x) (math-anglep x))
1372 (or (math-infinitep y) (math-anglep y)))
1387 (calcFunc-arctan2 y x)
1388 '(var nan var-nan)))
1389 (calc-record-why 'anglep y)
1390 (list 'calcFunc-arctan2 y x))))
1392 (defun math-arctan2-raw (y x) ; [F R R]
1393 (cond ((math-zerop y)
1394 (if (math-negp x) (math-pi)
1395 (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0)))
1399 (math-neg (math-pi-over-2))))
1401 (math-arctan-raw (math-div-float y x)))
1403 (math-add-float (math-arctan-raw (math-div-float y x))
1406 (math-sub-float (math-arctan-raw (math-div-float y x))
1409 (defun calcFunc-arcsincos (x) ; [V N] [Public]
1410 (if (and (Math-vectorp x)
1412 (calcFunc-arctan2 (nth 2 x) (nth 1 x))
1413 (math-reject-arg x "*Two-element vector expected")))
1417 ;;; Exponential function.
1419 (defun calcFunc-exp (x) ; [N N] [Public]
1421 ((and (memq x '(1 -1)) calc-symbolic-mode)
1422 (if (eq x 1) '(var e var-e) (math-div 1 '(var e var-e))))
1424 (math-with-extra-prec 2 (math-exp-raw (math-float x))))
1425 ((eq (car-safe x) 'sdev)
1426 (let ((ex (calcFunc-exp (nth 1 x))))
1427 (math-make-sdev ex (math-mul (nth 2 x) ex))))
1428 ((eq (car-safe x) 'intv)
1429 (math-make-intv (nth 1 x) (calcFunc-exp (nth 2 x))
1430 (calcFunc-exp (nth 3 x))))
1431 ((equal x '(var inf var-inf))
1433 ((equal x '(neg (var inf var-inf)))
1435 ((equal x '(var nan var-nan))
1437 (t (calc-record-why 'numberp x)
1438 (list 'calcFunc-exp x))))
1440 (defun calcFunc-expm1 (x) ; [N N] [Public]
1442 ((math-zerop x) '(float 0 0))
1443 (calc-symbolic-mode (signal 'inexact-result nil))
1445 (math-with-extra-prec 2
1446 (let ((x (math-float x)))
1447 (if (and (eq (car x) 'float)
1448 (math-lessp-float x '(float 1 0))
1449 (math-lessp-float '(float -1 0) x))
1450 (math-exp-minus-1-raw x)
1451 (math-add (math-exp-raw x) -1)))))
1452 ((eq (car-safe x) 'sdev)
1454 (let ((ex (calcFunc-expm1 (nth 1 x))))
1455 (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1))))
1456 (math-make-sdev (calcFunc-expm1 (nth 1 x))
1457 (math-mul (nth 2 x) (calcFunc-exp (nth 1 x))))))
1458 ((eq (car-safe x) 'intv)
1459 (math-make-intv (nth 1 x)
1460 (calcFunc-expm1 (nth 2 x))
1461 (calcFunc-expm1 (nth 3 x))))
1462 ((equal x '(var inf var-inf))
1464 ((equal x '(neg (var inf var-inf)))
1466 ((equal x '(var nan var-nan))
1468 (t (calc-record-why 'numberp x)
1469 (list 'calcFunc-expm1 x))))
1471 (defun calcFunc-exp10 (x) ; [N N] [Public]
1474 (math-pow '(float 1 1) x)))
1476 (defun math-exp-raw (x) ; [N N]
1477 (cond ((math-zerop x) '(float 1 0))
1478 (calc-symbolic-mode (signal 'inexact-result nil))
1480 (let ((expx (math-exp-raw (nth 1 x)))
1481 (sc (math-sin-cos-raw (nth 2 x))))
1483 (math-mul-float expx (cdr sc))
1484 (math-mul-float expx (car sc)))))
1485 ((eq (car x) 'polar)
1486 (let ((xc (math-complex x)))
1488 (math-exp-raw (nth 1 xc))
1489 (math-from-radians (nth 2 xc)))))
1490 ((math-use-emacs-fn 'exp x))
1491 ((or (math-lessp-float '(float 5 -1) x)
1492 (math-lessp-float x '(float -5 -1)))
1493 (if (math-lessp-float '(float 921035 1) x)
1495 (if (math-lessp-float x '(float -921035 1))
1497 (let* ((two-x (math-mul-float x '(float 2 0)))
1498 (hint (math-scale-int (nth 1 two-x) (nth 2 two-x)))
1499 (hfrac (math-sub-float x (math-mul-float (math-float hint)
1501 (math-mul-float (math-ipow (math-sqrt-e) hint)
1502 (math-add-float '(float 1 0)
1503 (math-exp-minus-1-raw hfrac)))))
1504 (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x)))))
1506 (defun math-exp-minus-1-raw (x) ; [F F]
1507 (math-exp-series x 2 3 x x))
1509 (defun math-exp-series (sum nfac n xpow x)
1510 (math-working "exp" sum)
1511 (let* ((nextx (math-mul-float xpow x))
1512 (nextsum (math-add-float sum (math-div-float nextx
1513 (math-float nfac)))))
1514 (if (math-nearly-equal-float sum nextsum)
1516 (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x))))
1522 (defun calcFunc-ln (x) ; [N N] [Public]
1523 (cond ((math-zerop x)
1524 (if calc-infinite-mode
1525 '(neg (var inf var-inf))
1526 (math-reject-arg x "*Logarithm of zero")))
1529 (math-with-extra-prec 2 (math-ln-raw (math-float x))))
1530 ((eq (car-safe x) 'sdev)
1531 (math-make-sdev (calcFunc-ln (nth 1 x))
1532 (math-div (nth 2 x) (nth 1 x))))
1533 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1534 (Math-zerop (nth 2 x))
1535 (not (math-intv-constp x))))
1536 (let ((calc-infinite-mode t))
1537 (math-make-intv (nth 1 x) (calcFunc-ln (nth 2 x))
1538 (calcFunc-ln (nth 3 x)))))
1539 ((equal x '(var e var-e))
1541 ((and (eq (car-safe x) '^)
1542 (equal (nth 1 x) '(var e var-e))
1543 (math-known-realp (nth 2 x)))
1546 (if (equal x '(var nan var-nan))
1548 '(var inf var-inf)))
1549 (t (calc-record-why 'numberp x)
1550 (list 'calcFunc-ln x))))
1552 (defun calcFunc-log10 (x) ; [N N] [Public]
1553 (cond ((math-equal-int x 1)
1554 (if (math-floatp x) '(float 0 0) 0))
1555 ((and (Math-integerp x)
1557 (let ((res (math-integer-log x 10)))
1559 (setq x (cdr res)))))
1561 ((and (eq (car-safe x) 'frac)
1563 (let ((res (math-integer-log (nth 2 x) 10)))
1565 (setq x (- (cdr res))))))
1568 (if calc-infinite-mode
1569 '(neg (var inf var-inf))
1570 (math-reject-arg x "*Logarithm of zero")))
1571 (calc-symbolic-mode (signal 'inexact-result nil))
1573 (math-with-extra-prec 2
1574 (let ((xf (math-float x)))
1575 (if (eq (nth 1 xf) 0)
1576 (math-reject-arg x "*Logarithm of zero"))
1577 (if (Math-integer-posp (nth 1 xf))
1578 (if (eq (nth 1 xf) 1) ; log10(1*10^n) = n
1579 (math-float (nth 2 xf))
1580 (let ((xdigs (1- (math-numdigs (nth 1 xf)))))
1582 (math-div-float (math-ln-raw-2
1583 (list 'float (nth 1 xf) (- xdigs)))
1585 (math-float (+ (nth 2 xf) xdigs)))))
1586 (math-div (calcFunc-ln xf) (math-ln-10))))))
1587 ((eq (car-safe x) 'sdev)
1588 (math-make-sdev (calcFunc-log10 (nth 1 x))
1590 (math-mul (nth 1 x) (math-ln-10)))))
1591 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1592 (not (math-intv-constp x))))
1593 (math-make-intv (nth 1 x)
1594 (calcFunc-log10 (nth 2 x))
1595 (calcFunc-log10 (nth 3 x))))
1597 (if (equal x '(var nan var-nan))
1599 '(var inf var-inf)))
1600 (t (calc-record-why 'numberp x)
1601 (list 'calcFunc-log10 x))))
1603 (defun calcFunc-log (x &optional b) ; [N N N] [Public]
1604 (cond ((or (null b) (equal b '(var e var-e)))
1605 (math-normalize (list 'calcFunc-ln x)))
1606 ((or (eq b 10) (equal b '(float 1 1)))
1607 (math-normalize (list 'calcFunc-log10 x)))
1609 (if calc-infinite-mode
1610 (math-div (calcFunc-ln x) (calcFunc-ln b))
1611 (math-reject-arg x "*Logarithm of zero")))
1613 (if calc-infinite-mode
1614 (math-div (calcFunc-ln x) (calcFunc-ln b))
1615 (math-reject-arg b "*Logarithm of zero")))
1616 ((math-equal-int b 1)
1617 (if calc-infinite-mode
1618 (math-div (calcFunc-ln x) 0)
1619 (math-reject-arg b "*Logarithm base one")))
1620 ((math-equal-int x 1)
1621 (if (math-floatp b) '(float 0 0) 0))
1622 ((and (Math-ratp x) (Math-ratp b)
1623 (math-posp x) (math-posp b)
1624 (let* ((sign 1) (inv nil)
1625 (xx (if (Math-lessp 1 x)
1629 (bb (if (Math-lessp 1 b)
1631 (setq sign (- sign))
1633 (res (if (Math-lessp xx bb)
1634 (setq inv (math-integer-log bb xx))
1635 (math-integer-log xx bb))))
1638 (math-div 1 (* sign (cdr res)))
1639 (* sign (cdr res)))))))
1641 (calc-symbolic-mode (signal 'inexact-result nil))
1642 ((and (Math-numberp x) (Math-numberp b))
1643 (math-with-extra-prec 2
1644 (math-div (math-ln-raw (math-float x))
1645 (math-log-base-raw b))))
1646 ((and (eq (car-safe x) 'sdev)
1648 (math-make-sdev (calcFunc-log (nth 1 x) b)
1651 (math-log-base-raw b)))))
1652 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1653 (not (math-intv-constp x)))
1655 (math-make-intv (nth 1 x)
1656 (calcFunc-log (nth 2 x) b)
1657 (calcFunc-log (nth 3 x) b)))
1658 ((or (eq (car-safe x) 'intv) (eq (car-safe b) 'intv))
1659 (math-div (calcFunc-ln x) (calcFunc-ln b)))
1660 ((or (math-infinitep x)
1662 (math-div (calcFunc-ln x) (calcFunc-ln b)))
1663 (t (if (Math-numberp b)
1664 (calc-record-why 'numberp x)
1665 (calc-record-why 'numberp b))
1666 (list 'calcFunc-log x b))))
1668 (defun calcFunc-alog (x &optional b)
1669 (cond ((or (null b) (equal b '(var e var-e)))
1670 (math-normalize (list 'calcFunc-exp x)))
1671 (t (math-pow b x))))
1673 (defun calcFunc-ilog (x b)
1674 (if (and (math-natnump x) (not (eq x 0))
1675 (math-natnump b) (not (eq b 0)))
1677 (math-reject-arg x "*Logarithm base one")
1678 (if (Math-natnum-lessp x b)
1680 (cdr (math-integer-log x b))))
1681 (math-floor (calcFunc-log x b))))
1683 (defun math-integer-log (x b)
1684 (let ((pows (list b))
1688 (while (not (Math-lessp x pow))
1689 (setq pows (cons pow pows)
1690 pow (math-sqr pow)))
1691 (setq n (lsh 1 (1- (length pows)))
1694 (while (and (setq pows (cdr pows))
1697 next (math-mul pow (car pows)))
1698 (or (Math-lessp x next)
1701 (cons (equal pow x) sum)))
1704 (defvar math-log-base-cache nil)
1705 (defun math-log-base-raw (b) ; [N N]
1706 (if (not (and (equal (car math-log-base-cache) b)
1707 (eq (nth 1 math-log-base-cache) calc-internal-prec)))
1708 (setq math-log-base-cache (list b calc-internal-prec
1709 (math-ln-raw (math-float b)))))
1710 (nth 2 math-log-base-cache))
1712 (defun calcFunc-lnp1 (x) ; [N N] [Public]
1713 (cond ((Math-equal-int x -1)
1714 (if calc-infinite-mode
1715 '(neg (var inf var-inf))
1716 (math-reject-arg x "*Logarithm of zero")))
1718 ((math-zerop x) '(float 0 0))
1719 (calc-symbolic-mode (signal 'inexact-result nil))
1721 (math-with-extra-prec 2
1722 (let ((x (math-float x)))
1723 (if (and (eq (car x) 'float)
1724 (math-lessp-float x '(float 5 -1))
1725 (math-lessp-float '(float -5 -1) x))
1726 (math-ln-plus-1-raw x)
1727 (math-ln-raw (math-add-float x '(float 1 0)))))))
1728 ((eq (car-safe x) 'sdev)
1729 (math-make-sdev (calcFunc-lnp1 (nth 1 x))
1730 (math-div (nth 2 x) (math-add (nth 1 x) 1))))
1731 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1732 (not (math-intv-constp x))))
1733 (math-make-intv (nth 1 x)
1734 (calcFunc-lnp1 (nth 2 x))
1735 (calcFunc-lnp1 (nth 3 x))))
1737 (if (equal x '(var nan var-nan))
1739 '(var inf var-inf)))
1740 (t (calc-record-why 'numberp x)
1741 (list 'calcFunc-lnp1 x))))
1743 (defun math-ln-raw (x) ; [N N] --- must be float format!
1744 (cond ((eq (car-safe x) 'cplx)
1746 (math-mul-float (math-ln-raw
1747 (math-add-float (math-sqr-float (nth 1 x))
1748 (math-sqr-float (nth 2 x))))
1750 (math-arctan2-raw (nth 2 x) (nth 1 x))))
1751 ((eq (car x) 'polar)
1752 (math-polar (list 'cplx
1753 (math-ln-raw (nth 1 x))
1754 (math-to-radians (nth 2 x)))))
1755 ((Math-equal-int x 1)
1757 (calc-symbolic-mode (signal 'inexact-result nil))
1758 ((math-posp (nth 1 x)) ; positive and real
1760 ((math-use-emacs-fn 'log x))
1762 (let ((xdigs (1- (math-numdigs (nth 1 x)))))
1763 (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs)))
1764 (math-mul-float (math-float (+ (nth 2 x) xdigs))
1767 (math-reject-arg x "*Logarithm of zero"))
1768 ((eq calc-complex-mode 'polar) ; negative and real
1770 (list 'cplx ; negative and real
1771 (math-ln-raw (math-neg-float x))
1773 (t (list 'cplx ; negative and real
1774 (math-ln-raw (math-neg-float x))
1777 (defun math-ln-raw-2 (x) ; [F F]
1778 (cond ((math-lessp-float '(float 14 -1) x)
1779 (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1)))
1781 (t ; now .7 < x <= 1.4
1782 (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
1783 (math-add-float x '(float 1 0)))))))
1785 (defun math-ln-raw-3 (x) ; [F F]
1786 (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
1789 ;;; Compute ln((1+x)/(1-x))
1790 (defun math-ln-raw-series (sum n x xsqr)
1791 (math-working "log" sum)
1792 (let* ((nextx (math-mul-float x xsqr))
1793 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1794 (if (math-nearly-equal-float sum nextsum)
1796 (math-ln-raw-series nextsum (+ n 2) nextx xsqr))))
1798 (defun math-ln-plus-1-raw (x)
1799 (math-lnp1-series x 2 x (math-neg x)))
1801 (defun math-lnp1-series (sum n xpow x)
1802 (math-working "lnp1" sum)
1803 (let* ((nextx (math-mul-float xpow x))
1804 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1805 (if (math-nearly-equal-float sum nextsum)
1807 (math-lnp1-series nextsum (1+ n) nextx x))))
1809 (defconst math-approx-ln-10
1810 (math-read-number-simple "2.302585092994045684018")
1811 "An approximation for ln(10).")
1813 (math-defcache math-ln-10 math-approx-ln-10
1814 (math-ln-raw-2 '(float 1 1)))
1816 (defconst math-approx-ln-2
1817 (math-read-number-simple "0.693147180559945309417")
1818 "An approximation for ln(2).")
1820 (math-defcache math-ln-2 math-approx-ln-2
1821 (math-ln-raw-3 (math-float '(frac 1 3))))
1825 ;;; Hyperbolic functions.
1827 (defun calcFunc-sinh (x) ; [N N] [Public]
1829 (math-expand-formulas
1831 (list '/ (list '- (list 'calcFunc-exp x)
1832 (list 'calcFunc-exp (list 'neg x))) 2)))
1834 (if calc-symbolic-mode (signal 'inexact-result nil))
1835 (math-with-extra-prec 2
1836 (let ((expx (math-exp-raw (math-float x))))
1837 (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1)))))
1838 ((eq (car-safe x) 'sdev)
1839 (math-make-sdev (calcFunc-sinh (nth 1 x))
1840 (math-mul (nth 2 x) (calcFunc-cosh (nth 1 x)))))
1842 (math-sort-intv (nth 1 x)
1843 (calcFunc-sinh (nth 2 x))
1844 (calcFunc-sinh (nth 3 x))))
1845 ((or (equal x '(var inf var-inf))
1846 (equal x '(neg (var inf var-inf)))
1847 (equal x '(var nan var-nan)))
1849 (t (calc-record-why 'numberp x)
1850 (list 'calcFunc-sinh x))))
1851 (put 'calcFunc-sinh 'math-expandable t)
1853 (defun calcFunc-cosh (x) ; [N N] [Public]
1855 (math-expand-formulas
1857 (list '/ (list '+ (list 'calcFunc-exp x)
1858 (list 'calcFunc-exp (list 'neg x))) 2)))
1860 (if calc-symbolic-mode (signal 'inexact-result nil))
1861 (math-with-extra-prec 2
1862 (let ((expx (math-exp-raw (math-float x))))
1863 (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1)))))
1864 ((eq (car-safe x) 'sdev)
1865 (math-make-sdev (calcFunc-cosh (nth 1 x))
1867 (calcFunc-sinh (nth 1 x)))))
1868 ((and (eq (car x) 'intv) (math-intv-constp x))
1869 (setq x (math-abs x))
1870 (math-sort-intv (nth 1 x)
1871 (calcFunc-cosh (nth 2 x))
1872 (calcFunc-cosh (nth 3 x))))
1873 ((or (equal x '(var inf var-inf))
1874 (equal x '(neg (var inf var-inf)))
1875 (equal x '(var nan var-nan)))
1877 (t (calc-record-why 'numberp x)
1878 (list 'calcFunc-cosh x))))
1879 (put 'calcFunc-cosh 'math-expandable t)
1881 (defun calcFunc-tanh (x) ; [N N] [Public]
1883 (math-expand-formulas
1885 (let ((expx (list 'calcFunc-exp x))
1886 (expmx (list 'calcFunc-exp (list 'neg x))))
1888 (list '/ (list '- expx expmx) (list '+ expx expmx))))))
1890 (if calc-symbolic-mode (signal 'inexact-result nil))
1891 (math-with-extra-prec 2
1892 (let* ((expx (calcFunc-exp (math-float x)))
1893 (expmx (math-div 1 expx)))
1894 (math-div (math-sub expx expmx)
1895 (math-add expx expmx)))))
1896 ((eq (car-safe x) 'sdev)
1897 (math-make-sdev (calcFunc-tanh (nth 1 x))
1899 (math-sqr (calcFunc-cosh (nth 1 x))))))
1901 (math-sort-intv (nth 1 x)
1902 (calcFunc-tanh (nth 2 x))
1903 (calcFunc-tanh (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-tanh x))))
1912 (put 'calcFunc-tanh 'math-expandable t)
1914 (defun calcFunc-sech (x) ; [N N] [Public]
1916 (math-expand-formulas
1918 (list '/ 2 (list '+ (list 'calcFunc-exp x)
1919 (list 'calcFunc-exp (list 'neg x))))))
1921 (if calc-symbolic-mode (signal 'inexact-result nil))
1922 (math-with-extra-prec 2
1923 (let ((expx (math-exp-raw (math-float x))))
1924 (math-div '(float 2 0) (math-add expx (math-div 1 expx))))))
1925 ((eq (car-safe x) 'sdev)
1926 (math-make-sdev (calcFunc-sech (nth 1 x))
1928 (math-mul (calcFunc-sech (nth 1 x))
1929 (calcFunc-tanh (nth 1 x))))))
1930 ((and (eq (car x) 'intv) (math-intv-constp x))
1931 (setq x (math-abs x))
1932 (math-sort-intv (nth 1 x)
1933 (calcFunc-sech (nth 2 x))
1934 (calcFunc-sech (nth 3 x))))
1935 ((or (equal x '(var inf var-inf))
1936 (equal x '(neg (var inf var-inf))))
1938 ((equal x '(var nan var-nan))
1940 (t (calc-record-why 'numberp x)
1941 (list 'calcFunc-sech x))))
1942 (put 'calcFunc-sech 'math-expandable t)
1944 (defun calcFunc-csch (x) ; [N N] [Public]
1945 (cond ((eq x 0) (math-div 1 0))
1946 (math-expand-formulas
1948 (list '/ 2 (list '- (list 'calcFunc-exp x)
1949 (list 'calcFunc-exp (list 'neg x))))))
1951 (if calc-symbolic-mode (signal 'inexact-result nil))
1952 (math-with-extra-prec 2
1953 (let ((expx (math-exp-raw (math-float x))))
1954 (math-div '(float 2 0) (math-add expx (math-div -1 expx))))))
1955 ((eq (car-safe x) 'sdev)
1956 (math-make-sdev (calcFunc-csch (nth 1 x))
1958 (math-mul (calcFunc-csch (nth 1 x))
1959 (calcFunc-coth (nth 1 x))))))
1961 (if (and (Math-negp (nth 2 x))
1962 (Math-posp (nth 3 x)))
1963 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1964 (math-sort-intv (nth 1 x)
1965 (calcFunc-csch (nth 2 x))
1966 (calcFunc-csch (nth 3 x)))))
1967 ((or (equal x '(var inf var-inf))
1968 (equal x '(neg (var inf var-inf))))
1970 ((equal x '(var nan var-nan))
1972 (t (calc-record-why 'numberp x)
1973 (list 'calcFunc-csch x))))
1974 (put 'calcFunc-csch 'math-expandable t)
1976 (defun calcFunc-coth (x) ; [N N] [Public]
1977 (cond ((eq x 0) (math-div 1 0))
1978 (math-expand-formulas
1980 (let ((expx (list 'calcFunc-exp x))
1981 (expmx (list 'calcFunc-exp (list 'neg x))))
1983 (list '/ (list '+ expx expmx) (list '- expx expmx))))))
1985 (if calc-symbolic-mode (signal 'inexact-result nil))
1986 (math-with-extra-prec 2
1987 (let* ((expx (calcFunc-exp (math-float x)))
1988 (expmx (math-div 1 expx)))
1989 (math-div (math-add expx expmx)
1990 (math-sub expx expmx)))))
1991 ((eq (car-safe x) 'sdev)
1992 (math-make-sdev (calcFunc-coth (nth 1 x))
1994 (math-sqr (calcFunc-sinh (nth 1 x))))))
1996 (if (and (Math-negp (nth 2 x))
1997 (Math-posp (nth 3 x)))
1998 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1999 (math-sort-intv (nth 1 x)
2000 (calcFunc-coth (nth 2 x))
2001 (calcFunc-coth (nth 3 x)))))
2002 ((equal x '(var inf var-inf))
2004 ((equal x '(neg (var inf var-inf)))
2006 ((equal x '(var nan var-nan))
2008 (t (calc-record-why 'numberp x)
2009 (list 'calcFunc-coth x))))
2010 (put 'calcFunc-coth 'math-expandable t)
2012 (defun calcFunc-arcsinh (x) ; [N N] [Public]
2014 (math-expand-formulas
2016 (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
2017 (list '+ (list '^ x 2) 1))))))
2019 (if calc-symbolic-mode (signal 'inexact-result nil))
2020 (math-with-extra-prec 2
2021 (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x)
2023 ((eq (car-safe x) 'sdev)
2024 (math-make-sdev (calcFunc-arcsinh (nth 1 x))
2027 (math-add (math-sqr (nth 1 x)) 1)))))
2029 (math-sort-intv (nth 1 x)
2030 (calcFunc-arcsinh (nth 2 x))
2031 (calcFunc-arcsinh (nth 3 x))))
2032 ((or (equal x '(var inf var-inf))
2033 (equal x '(neg (var inf var-inf)))
2034 (equal x '(var nan var-nan)))
2036 (t (calc-record-why 'numberp x)
2037 (list 'calcFunc-arcsinh x))))
2038 (put 'calcFunc-arcsinh 'math-expandable t)
2040 (defun calcFunc-arccosh (x) ; [N N] [Public]
2042 ((and (eq x -1) calc-symbolic-mode)
2044 ((and (eq x 0) calc-symbolic-mode)
2045 (math-div (math-mul '(var pi var-pi) '(var i var-i)) 2))
2046 (math-expand-formulas
2048 (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
2049 (list '- (list '^ x 2) 1))))))
2051 (if calc-symbolic-mode (signal 'inexact-result nil))
2052 (if (Math-equal-int x -1)
2053 (math-imaginary (math-pi))
2054 (math-with-extra-prec 2
2055 (if (or t ; need to do this even in the real case!
2056 (memq (car-safe x) '(cplx polar)))
2057 (let ((xp1 (math-add 1 x))) ; this gets the branch cuts right
2059 (math-add x (math-mul xp1
2066 (math-add x (math-sqrt-raw (math-add (math-sqr x)
2067 '(float -1 0)))))))))
2068 ((eq (car-safe x) 'sdev)
2069 (math-make-sdev (calcFunc-arccosh (nth 1 x))
2072 (math-add (math-sqr (nth 1 x)) -1)))))
2074 (math-sort-intv (nth 1 x)
2075 (calcFunc-arccosh (nth 2 x))
2076 (calcFunc-arccosh (nth 3 x))))
2077 ((or (equal x '(var inf var-inf))
2078 (equal x '(neg (var inf var-inf)))
2079 (equal x '(var nan var-nan)))
2081 (t (calc-record-why 'numberp x)
2082 (list 'calcFunc-arccosh x))))
2083 (put 'calcFunc-arccosh 'math-expandable t)
2085 (defun calcFunc-arctanh (x) ; [N N] [Public]
2087 ((and (Math-equal-int x 1) calc-infinite-mode)
2089 ((and (Math-equal-int x -1) calc-infinite-mode)
2090 '(neg (var inf var-inf)))
2091 (math-expand-formulas
2093 (list 'calcFunc-ln (list '+ 1 x))
2094 (list 'calcFunc-ln (list '- 1 x))) 2))
2096 (if calc-symbolic-mode (signal 'inexact-result nil))
2097 (math-with-extra-prec 2
2098 (if (or (memq (car-safe x) '(cplx polar))
2100 (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x))
2101 (math-ln-raw (math-sub '(float 1 0) x)))
2103 (if (and (math-equal-int x 1) calc-infinite-mode)
2105 (if (and (math-equal-int x -1) calc-infinite-mode)
2106 '(neg (var inf var-inf))
2107 (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x)
2110 ((eq (car-safe x) 'sdev)
2111 (math-make-sdev (calcFunc-arctanh (nth 1 x))
2113 (math-sub 1 (math-sqr (nth 1 x))))))
2115 (math-sort-intv (nth 1 x)
2116 (calcFunc-arctanh (nth 2 x))
2117 (calcFunc-arctanh (nth 3 x))))
2118 ((equal x '(var nan var-nan))
2120 (t (calc-record-why 'numberp x)
2121 (list 'calcFunc-arctanh x))))
2122 (put 'calcFunc-arctanh 'math-expandable t)
2125 ;;; Convert A from HMS or degrees to radians.
2126 (defun calcFunc-rad (a) ; [R R] [Public]
2127 (cond ((or (Math-numberp a)
2129 (math-with-extra-prec 2
2130 (math-mul a (math-pi-over-180))))
2132 (math-from-hms a 'rad))
2134 (math-make-sdev (calcFunc-rad (nth 1 a))
2135 (calcFunc-rad (nth 2 a))))
2136 (math-expand-formulas
2137 (math-div (math-mul a '(var pi var-pi)) 180))
2138 ((math-infinitep a) a)
2139 (t (list 'calcFunc-rad a))))
2140 (put 'calcFunc-rad 'math-expandable t)
2142 ;;; Convert A from HMS or radians to degrees.
2143 (defun calcFunc-deg (a) ; [R R] [Public]
2144 (cond ((or (Math-numberp a)
2146 (math-with-extra-prec 2
2147 (math-div a (math-pi-over-180))))
2149 (math-from-hms a 'deg))
2151 (math-make-sdev (calcFunc-deg (nth 1 a))
2152 (calcFunc-deg (nth 2 a))))
2153 (math-expand-formulas
2154 (math-div (math-mul 180 a) '(var pi var-pi)))
2155 ((math-infinitep a) a)
2156 (t (list 'calcFunc-deg a))))
2157 (put 'calcFunc-deg 'math-expandable t)
2159 (provide 'calc-math)
2161 ;;; arch-tag: c7367e8e-d0b8-4f70-8577-2fb3f31dbb4c
2162 ;;; calc-math.el ends here