1 ;;; calc-math.el --- mathematical functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
4 ;; 2005, 2006, 2007, 2008, 2009 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 of the License, or
14 ;; (at your option) any later version.
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. If not, see <http://www.gnu.org/licenses/>.
28 ;; This file is autoloaded from calc-ext.el.
34 ;;; Find out how many 9s in 9.9999... will give distinct Emacs floats,
35 ;;; then back off by one.
37 (defvar math-emacs-precision
40 (xx (+ x
(* 9 (expt 10 (- n
))))))
45 (setq xx
(+ x
(* 9 (expt 10 (- n
)))))))
47 "The number of digits in an Emacs float.")
49 ;;; Find the largest power of 10 which is an Emacs float,
50 ;;; then back off by one so that any float d.dddd...eN
51 ;;; is an Emacs float, for acceptable d.dddd....
53 (defvar math-largest-emacs-expt
56 ;; The following loop is for efficiency; it should stop when
57 ;; 10^(2x) is too large. This could be indicated by a range
58 ;; error when computing 10^(2x) or an infinite value for 10^(2x).
63 (setq pow
(condition-case nil
66 ;; The following loop should stop when 10^(x+1) is too large.
67 (setq pow
(condition-case nil
74 (setq pow
(condition-case nil
78 "The largest exponent which Calc will convert to an Emacs float.")
80 (defvar math-smallest-emacs-expt
82 (while (condition-case nil
87 (while (condition-case nil
92 "The smallest exponent which Calc will convert to an Emacs float.")
94 (defun math-use-emacs-fn (fn x
)
95 "Use the native Emacs function FN to evaluate the Calc number X.
96 If this can't be done, return NIL."
98 (<= calc-internal-prec math-emacs-precision
)
100 (let* ((fx (math-float x
))
101 (xpon (+ (nth 2 x
) (1- (math-numdigs (nth 1 x
))))))
102 (and (<= math-smallest-emacs-expt xpon
)
103 (<= xpon math-largest-emacs-expt
)
109 (let ((calc-number-radix 10))
110 (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 &optional orgx
) ; [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
) (if orgx orgx 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)) (if orgx orgx x
)))
1095 (t (math-sin-raw-2 x
(if orgx orgx 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
) 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-with-extra-prec -
1 (math-nearly-zerop-float x orgx
))
1221 ((math-use-emacs-fn 'sin x
))
1222 (calc-symbolic-mode (signal 'inexact-result nil
))
1223 (t (math-sin-series x
6 4 x
(math-neg-float (math-sqr-float x
)))))))
1225 (defun math-cos-raw-2 (x orgx
) ; [F F]
1226 (cond ((math-with-extra-prec -
1 (math-nearly-zerop-float x orgx
))
1228 ((math-use-emacs-fn 'cos x
))
1229 (calc-symbolic-mode (signal 'inexact-result nil
))
1230 (t (let ((xnegsqr (math-neg-float (math-sqr-float x
))))
1232 (math-add-float '(float 1 0)
1233 (math-mul-float xnegsqr
'(float 5 -
1)))
1234 24 5 xnegsqr xnegsqr
)))))
1236 (defun math-sin-series (sum nfac n x xnegsqr
)
1237 (math-working "sin" sum
)
1238 (let* ((nextx (math-mul-float x xnegsqr
))
1239 (nextsum (math-add-float sum
(math-div-float nextx
1240 (math-float nfac
)))))
1241 (if (math-nearly-equal-float sum nextsum
)
1243 (math-sin-series nextsum
(math-mul nfac
(* n
(1+ n
)))
1244 (+ n
2) nextx xnegsqr
))))
1247 ;;; Inverse sine, cosine, tangent.
1249 (defun calcFunc-arcsin (x) ; [N N] [Public]
1251 ((and (eq x
1) (eq calc-angle-mode
'deg
)) 90)
1252 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) -
90)
1253 (calc-symbolic-mode (signal 'inexact-result nil
))
1255 (math-with-extra-prec 2
1256 (math-from-radians (math-arcsin-raw (math-float x
)))))
1258 (math-make-sdev (calcFunc-arcsin (nth 1 x
))
1262 (math-sub 1 (math-sqr (nth 1 x
))))))))
1264 (math-sort-intv (nth 1 x
)
1265 (calcFunc-arcsin (nth 2 x
))
1266 (calcFunc-arcsin (nth 3 x
))))
1267 ((equal x
'(var nan var-nan
))
1269 (t (calc-record-why 'numberp x
)
1270 (list 'calcFunc-arcsin x
))))
1272 (defun calcFunc-arccos (x) ; [N N] [Public]
1274 ((and (eq x
0) (eq calc-angle-mode
'deg
)) 90)
1275 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) 180)
1276 (calc-symbolic-mode (signal 'inexact-result nil
))
1278 (math-with-extra-prec 2
1279 (math-from-radians (math-arccos-raw (math-float x
)))))
1281 (math-make-sdev (calcFunc-arccos (nth 1 x
))
1285 (math-sub 1 (math-sqr (nth 1 x
))))))))
1287 (math-sort-intv (nth 1 x
)
1288 (calcFunc-arccos (nth 2 x
))
1289 (calcFunc-arccos (nth 3 x
))))
1290 ((equal x
'(var nan var-nan
))
1292 (t (calc-record-why 'numberp x
)
1293 (list 'calcFunc-arccos x
))))
1295 (defun calcFunc-arctan (x) ; [N N] [Public]
1297 ((and (eq x
1) (eq calc-angle-mode
'deg
)) 45)
1298 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) -
45)
1300 (math-with-extra-prec 2
1301 (math-from-radians (math-arctan-raw (math-float x
)))))
1303 (math-make-sdev (calcFunc-arctan (nth 1 x
))
1306 (math-add 1 (math-sqr (nth 1 x
)))))))
1308 (math-sort-intv (nth 1 x
)
1309 (calcFunc-arctan (nth 2 x
))
1310 (calcFunc-arctan (nth 3 x
))))
1311 ((equal x
'(var inf var-inf
))
1312 (math-quarter-circle t
))
1313 ((equal x
'(neg (var inf var-inf
)))
1314 (math-neg (math-quarter-circle t
)))
1315 ((equal x
'(var nan var-nan
))
1317 (t (calc-record-why 'numberp x
)
1318 (list 'calcFunc-arctan x
))))
1320 (defun math-arcsin-raw (x) ; [N N]
1321 (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x
)))))
1322 (if (or (memq (car x
) '(cplx polar
))
1323 (memq (car a
) '(cplx polar
)))
1324 (math-with-extra-prec 2 ; use extra precision for difficult case
1325 (math-mul '(cplx 0 -
1)
1326 (math-ln-raw (math-add (math-mul '(cplx 0 1) x
) a
))))
1327 (math-arctan2-raw x a
))))
1329 (defun math-arccos-raw (x) ; [N N]
1330 (math-sub (math-pi-over-2) (math-arcsin-raw x
)))
1332 (defun math-arctan-raw (x) ; [N N]
1333 (cond ((memq (car x
) '(cplx polar
))
1334 (math-with-extra-prec 2 ; extra-extra
1336 (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x
)))
1337 (math-ln-raw (math-add 1 (math-mul '(cplx 0 -
1) x
))))
1339 ((Math-integer-negp (nth 1 x
))
1340 (math-neg-float (math-arctan-raw (math-neg-float x
))))
1342 ((math-use-emacs-fn 'atan x
))
1343 (calc-symbolic-mode (signal 'inexact-result nil
))
1344 ((math-equal-int x
1) (math-pi-over-4))
1345 ((math-equal-int x -
1) (math-neg (math-pi-over-4)))
1346 ((math-lessp-float '(float 414214 -
6) x
) ; if x > sqrt(2) - 1, reduce
1347 (if (math-lessp-float '(float 1 0) x
)
1348 (math-sub-float (math-mul-float (math-pi) '(float 5 -
1))
1349 (math-arctan-raw (math-div-float '(float 1 0) x
)))
1350 (math-sub-float (math-mul-float (math-pi) '(float 25 -
2))
1351 (math-arctan-raw (math-div-float
1352 (math-sub-float '(float 1 0) x
)
1353 (math-add-float '(float 1 0)
1355 (t (math-arctan-series x
3 x
(math-neg-float (math-sqr-float x
))))))
1357 (defun math-arctan-series (sum n x xnegsqr
)
1358 (math-working "arctan" sum
)
1359 (let* ((nextx (math-mul-float x xnegsqr
))
1360 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1361 (if (math-nearly-equal-float sum nextsum
)
1363 (math-arctan-series nextsum
(+ n
2) nextx xnegsqr
))))
1365 (defun calcFunc-arctan2 (y x
) ; [F R R] [Public]
1368 (math-with-extra-prec 2
1369 (math-from-radians (math-arctan2-raw (math-float y
)
1371 (calc-record-why 'anglep x
)
1372 (list 'calcFunc-arctan2 y x
))
1373 (if (and (or (math-infinitep x
) (math-anglep x
))
1374 (or (math-infinitep y
) (math-anglep y
)))
1389 (calcFunc-arctan2 y x
)
1390 '(var nan var-nan
)))
1391 (calc-record-why 'anglep y
)
1392 (list 'calcFunc-arctan2 y x
))))
1394 (defun math-arctan2-raw (y x
) ; [F R R]
1395 (cond ((math-zerop y
)
1396 (if (math-negp x
) (math-pi)
1397 (if (or (math-floatp x
) (math-floatp y
)) '(float 0 0) 0)))
1401 (math-neg (math-pi-over-2))))
1403 (math-arctan-raw (math-div-float y x
)))
1405 (math-add-float (math-arctan-raw (math-div-float y x
))
1408 (math-sub-float (math-arctan-raw (math-div-float y x
))
1411 (defun calcFunc-arcsincos (x) ; [V N] [Public]
1412 (if (and (Math-vectorp x
)
1414 (calcFunc-arctan2 (nth 2 x
) (nth 1 x
))
1415 (math-reject-arg x
"*Two-element vector expected")))
1419 ;;; Exponential function.
1421 (defun calcFunc-exp (x) ; [N N] [Public]
1423 ((and (memq x
'(1 -
1)) calc-symbolic-mode
)
1424 (if (eq x
1) '(var e var-e
) (math-div 1 '(var e var-e
))))
1426 (math-with-extra-prec 2 (math-exp-raw (math-float x
))))
1427 ((eq (car-safe x
) 'sdev
)
1428 (let ((ex (calcFunc-exp (nth 1 x
))))
1429 (math-make-sdev ex
(math-mul (nth 2 x
) ex
))))
1430 ((eq (car-safe x
) 'intv
)
1431 (math-make-intv (nth 1 x
) (calcFunc-exp (nth 2 x
))
1432 (calcFunc-exp (nth 3 x
))))
1433 ((equal x
'(var inf var-inf
))
1435 ((equal x
'(neg (var inf var-inf
)))
1437 ((equal x
'(var nan var-nan
))
1439 (t (calc-record-why 'numberp x
)
1440 (list 'calcFunc-exp x
))))
1442 (defun calcFunc-expm1 (x) ; [N N] [Public]
1444 ((math-zerop x
) '(float 0 0))
1445 (calc-symbolic-mode (signal 'inexact-result nil
))
1447 (math-with-extra-prec 2
1448 (let ((x (math-float x
)))
1449 (if (and (eq (car x
) 'float
)
1450 (math-lessp-float x
'(float 1 0))
1451 (math-lessp-float '(float -
1 0) x
))
1452 (math-exp-minus-1-raw x
)
1453 (math-add (math-exp-raw x
) -
1)))))
1454 ((eq (car-safe x
) 'sdev
)
1456 (let ((ex (calcFunc-expm1 (nth 1 x
))))
1457 (math-make-sdev ex
(math-mul (nth 2 x
) (math-add ex
1))))
1458 (math-make-sdev (calcFunc-expm1 (nth 1 x
))
1459 (math-mul (nth 2 x
) (calcFunc-exp (nth 1 x
))))))
1460 ((eq (car-safe x
) 'intv
)
1461 (math-make-intv (nth 1 x
)
1462 (calcFunc-expm1 (nth 2 x
))
1463 (calcFunc-expm1 (nth 3 x
))))
1464 ((equal x
'(var inf var-inf
))
1466 ((equal x
'(neg (var inf var-inf
)))
1468 ((equal x
'(var nan var-nan
))
1470 (t (calc-record-why 'numberp x
)
1471 (list 'calcFunc-expm1 x
))))
1473 (defun calcFunc-exp10 (x) ; [N N] [Public]
1476 (math-pow '(float 1 1) x
)))
1478 (defun math-exp-raw (x) ; [N N]
1479 (cond ((math-zerop x
) '(float 1 0))
1480 (calc-symbolic-mode (signal 'inexact-result nil
))
1482 (let ((expx (math-exp-raw (nth 1 x
)))
1483 (sc (math-sin-cos-raw (nth 2 x
))))
1485 (math-mul-float expx
(cdr sc
))
1486 (math-mul-float expx
(car sc
)))))
1487 ((eq (car x
) 'polar
)
1488 (let ((xc (math-complex x
)))
1490 (math-exp-raw (nth 1 xc
))
1491 (math-from-radians (nth 2 xc
)))))
1492 ((math-use-emacs-fn 'exp x
))
1493 ((or (math-lessp-float '(float 5 -
1) x
)
1494 (math-lessp-float x
'(float -
5 -
1)))
1495 (if (math-lessp-float '(float 921035 1) x
)
1497 (if (math-lessp-float x
'(float -
921035 1))
1499 (let* ((two-x (math-mul-float x
'(float 2 0)))
1500 (hint (math-scale-int (nth 1 two-x
) (nth 2 two-x
)))
1501 (hfrac (math-sub-float x
(math-mul-float (math-float hint
)
1503 (math-mul-float (math-ipow (math-sqrt-e) hint
)
1504 (math-add-float '(float 1 0)
1505 (math-exp-minus-1-raw hfrac
)))))
1506 (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x
)))))
1508 (defun math-exp-minus-1-raw (x) ; [F F]
1509 (math-exp-series x
2 3 x x
))
1511 (defun math-exp-series (sum nfac n xpow x
)
1512 (math-working "exp" sum
)
1513 (let* ((nextx (math-mul-float xpow x
))
1514 (nextsum (math-add-float sum
(math-div-float nextx
1515 (math-float nfac
)))))
1516 (if (math-nearly-equal-float sum nextsum
)
1518 (math-exp-series nextsum
(math-mul nfac n
) (1+ n
) nextx x
))))
1524 (defun calcFunc-ln (x) ; [N N] [Public]
1525 (cond ((math-zerop x
)
1526 (if calc-infinite-mode
1527 '(neg (var inf var-inf
))
1528 (math-reject-arg x
"*Logarithm of zero")))
1531 (math-with-extra-prec 2 (math-ln-raw (math-float x
))))
1532 ((eq (car-safe x
) 'sdev
)
1533 (math-make-sdev (calcFunc-ln (nth 1 x
))
1534 (math-div (nth 2 x
) (nth 1 x
))))
1535 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1536 (Math-zerop (nth 2 x
))
1537 (not (math-intv-constp x
))))
1538 (let ((calc-infinite-mode t
))
1539 (math-make-intv (nth 1 x
) (calcFunc-ln (nth 2 x
))
1540 (calcFunc-ln (nth 3 x
)))))
1541 ((equal x
'(var e var-e
))
1543 ((and (eq (car-safe x
) '^
)
1544 (equal (nth 1 x
) '(var e var-e
))
1545 (math-known-realp (nth 2 x
)))
1548 (if (equal x
'(var nan var-nan
))
1550 '(var inf var-inf
)))
1551 (t (calc-record-why 'numberp x
)
1552 (list 'calcFunc-ln x
))))
1554 (defun calcFunc-log10 (x) ; [N N] [Public]
1555 (cond ((math-equal-int x
1)
1556 (if (math-floatp x
) '(float 0 0) 0))
1557 ((and (Math-integerp x
)
1559 (let ((res (math-integer-log x
10)))
1561 (setq x
(cdr res
)))))
1563 ((and (eq (car-safe x
) 'frac
)
1565 (let ((res (math-integer-log (nth 2 x
) 10)))
1567 (setq x
(- (cdr res
))))))
1570 (if calc-infinite-mode
1571 '(neg (var inf var-inf
))
1572 (math-reject-arg x
"*Logarithm of zero")))
1573 (calc-symbolic-mode (signal 'inexact-result nil
))
1575 (math-with-extra-prec 2
1576 (let ((xf (math-float x
)))
1577 (if (eq (nth 1 xf
) 0)
1578 (math-reject-arg x
"*Logarithm of zero"))
1579 (if (Math-integer-posp (nth 1 xf
))
1580 (if (eq (nth 1 xf
) 1) ; log10(1*10^n) = n
1581 (math-float (nth 2 xf
))
1582 (let ((xdigs (1- (math-numdigs (nth 1 xf
)))))
1584 (math-div-float (math-ln-raw-2
1585 (list 'float
(nth 1 xf
) (- xdigs
)))
1587 (math-float (+ (nth 2 xf
) xdigs
)))))
1588 (math-div (calcFunc-ln xf
) (math-ln-10))))))
1589 ((eq (car-safe x
) 'sdev
)
1590 (math-make-sdev (calcFunc-log10 (nth 1 x
))
1592 (math-mul (nth 1 x
) (math-ln-10)))))
1593 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1594 (not (math-intv-constp x
))))
1595 (math-make-intv (nth 1 x
)
1596 (calcFunc-log10 (nth 2 x
))
1597 (calcFunc-log10 (nth 3 x
))))
1599 (if (equal x
'(var nan var-nan
))
1601 '(var inf var-inf
)))
1602 (t (calc-record-why 'numberp x
)
1603 (list 'calcFunc-log10 x
))))
1605 (defun calcFunc-log (x &optional b
) ; [N N N] [Public]
1606 (cond ((or (null b
) (equal b
'(var e var-e
)))
1607 (math-normalize (list 'calcFunc-ln x
)))
1608 ((or (eq b
10) (equal b
'(float 1 1)))
1609 (math-normalize (list 'calcFunc-log10 x
)))
1611 (if calc-infinite-mode
1612 (math-div (calcFunc-ln x
) (calcFunc-ln b
))
1613 (math-reject-arg x
"*Logarithm of zero")))
1615 (if calc-infinite-mode
1616 (math-div (calcFunc-ln x
) (calcFunc-ln b
))
1617 (math-reject-arg b
"*Logarithm of zero")))
1618 ((math-equal-int b
1)
1619 (if calc-infinite-mode
1620 (math-div (calcFunc-ln x
) 0)
1621 (math-reject-arg b
"*Logarithm base one")))
1622 ((math-equal-int x
1)
1623 (if (math-floatp b
) '(float 0 0) 0))
1624 ((and (Math-ratp x
) (Math-ratp b
)
1625 (math-posp x
) (math-posp b
)
1626 (let* ((sign 1) (inv nil
)
1627 (xx (if (Math-lessp 1 x
)
1631 (bb (if (Math-lessp 1 b
)
1633 (setq sign
(- sign
))
1635 (res (if (Math-lessp xx bb
)
1636 (setq inv
(math-integer-log bb xx
))
1637 (math-integer-log xx bb
))))
1640 (math-div 1 (* sign
(cdr res
)))
1641 (* sign
(cdr res
)))))))
1643 (calc-symbolic-mode (signal 'inexact-result nil
))
1644 ((and (Math-numberp x
) (Math-numberp b
))
1645 (math-with-extra-prec 2
1646 (math-div (math-ln-raw (math-float x
))
1647 (math-log-base-raw b
))))
1648 ((and (eq (car-safe x
) 'sdev
)
1650 (math-make-sdev (calcFunc-log (nth 1 x
) b
)
1653 (math-log-base-raw b
)))))
1654 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1655 (not (math-intv-constp x
)))
1657 (math-make-intv (nth 1 x
)
1658 (calcFunc-log (nth 2 x
) b
)
1659 (calcFunc-log (nth 3 x
) b
)))
1660 ((or (eq (car-safe x
) 'intv
) (eq (car-safe b
) 'intv
))
1661 (math-div (calcFunc-ln x
) (calcFunc-ln b
)))
1662 ((or (math-infinitep x
)
1664 (math-div (calcFunc-ln x
) (calcFunc-ln b
)))
1665 (t (if (Math-numberp b
)
1666 (calc-record-why 'numberp x
)
1667 (calc-record-why 'numberp b
))
1668 (list 'calcFunc-log x b
))))
1670 (defun calcFunc-alog (x &optional b
)
1671 (cond ((or (null b
) (equal b
'(var e var-e
)))
1672 (math-normalize (list 'calcFunc-exp x
)))
1673 (t (math-pow b x
))))
1675 (defun calcFunc-ilog (x b
)
1676 (if (and (math-natnump x
) (not (eq x
0))
1677 (math-natnump b
) (not (eq b
0)))
1679 (math-reject-arg x
"*Logarithm base one")
1680 (if (Math-natnum-lessp x b
)
1682 (cdr (math-integer-log x b
))))
1683 (math-floor (calcFunc-log x b
))))
1685 (defun math-integer-log (x b
)
1686 (let ((pows (list b
))
1690 (while (not (Math-lessp x pow
))
1691 (setq pows
(cons pow pows
)
1692 pow
(math-sqr pow
)))
1693 (setq n
(lsh 1 (1- (length pows
)))
1696 (while (and (setq pows
(cdr pows
))
1699 next
(math-mul pow
(car pows
)))
1700 (or (Math-lessp x next
)
1703 (cons (equal pow x
) sum
)))
1706 (defvar math-log-base-cache nil
)
1707 (defun math-log-base-raw (b) ; [N N]
1708 (if (not (and (equal (car math-log-base-cache
) b
)
1709 (eq (nth 1 math-log-base-cache
) calc-internal-prec
)))
1710 (setq math-log-base-cache
(list b calc-internal-prec
1711 (math-ln-raw (math-float b
)))))
1712 (nth 2 math-log-base-cache
))
1714 (defun calcFunc-lnp1 (x) ; [N N] [Public]
1715 (cond ((Math-equal-int x -
1)
1716 (if calc-infinite-mode
1717 '(neg (var inf var-inf
))
1718 (math-reject-arg x
"*Logarithm of zero")))
1720 ((math-zerop x
) '(float 0 0))
1721 (calc-symbolic-mode (signal 'inexact-result nil
))
1723 (math-with-extra-prec 2
1724 (let ((x (math-float x
)))
1725 (if (and (eq (car x
) 'float
)
1726 (math-lessp-float x
'(float 5 -
1))
1727 (math-lessp-float '(float -
5 -
1) x
))
1728 (math-ln-plus-1-raw x
)
1729 (math-ln-raw (math-add-float x
'(float 1 0)))))))
1730 ((eq (car-safe x
) 'sdev
)
1731 (math-make-sdev (calcFunc-lnp1 (nth 1 x
))
1732 (math-div (nth 2 x
) (math-add (nth 1 x
) 1))))
1733 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1734 (not (math-intv-constp x
))))
1735 (math-make-intv (nth 1 x
)
1736 (calcFunc-lnp1 (nth 2 x
))
1737 (calcFunc-lnp1 (nth 3 x
))))
1739 (if (equal x
'(var nan var-nan
))
1741 '(var inf var-inf
)))
1742 (t (calc-record-why 'numberp x
)
1743 (list 'calcFunc-lnp1 x
))))
1745 (defun math-ln-raw (x) ; [N N] --- must be float format!
1746 (cond ((eq (car-safe x
) 'cplx
)
1748 (math-mul-float (math-ln-raw
1749 (math-add-float (math-sqr-float (nth 1 x
))
1750 (math-sqr-float (nth 2 x
))))
1752 (math-arctan2-raw (nth 2 x
) (nth 1 x
))))
1753 ((eq (car x
) 'polar
)
1754 (math-polar (list 'cplx
1755 (math-ln-raw (nth 1 x
))
1756 (math-to-radians (nth 2 x
)))))
1757 ((Math-equal-int x
1)
1759 (calc-symbolic-mode (signal 'inexact-result nil
))
1760 ((math-posp (nth 1 x
)) ; positive and real
1762 ((math-use-emacs-fn 'log x
))
1764 (let ((xdigs (1- (math-numdigs (nth 1 x
)))))
1765 (math-add-float (math-ln-raw-2 (list 'float
(nth 1 x
) (- xdigs
)))
1766 (math-mul-float (math-float (+ (nth 2 x
) xdigs
))
1769 (math-reject-arg x
"*Logarithm of zero"))
1770 ((eq calc-complex-mode
'polar
) ; negative and real
1772 (list 'cplx
; negative and real
1773 (math-ln-raw (math-neg-float x
))
1775 (t (list 'cplx
; negative and real
1776 (math-ln-raw (math-neg-float x
))
1779 (defun math-ln-raw-2 (x) ; [F F]
1780 (cond ((math-lessp-float '(float 14 -
1) x
)
1781 (math-add-float (math-ln-raw-2 (math-mul-float x
'(float 5 -
1)))
1783 (t ; now .7 < x <= 1.4
1784 (math-ln-raw-3 (math-div-float (math-sub-float x
'(float 1 0))
1785 (math-add-float x
'(float 1 0)))))))
1787 (defun math-ln-raw-3 (x) ; [F F]
1788 (math-mul-float (math-ln-raw-series x
3 x
(math-sqr-float x
))
1791 ;;; Compute ln((1+x)/(1-x))
1792 (defun math-ln-raw-series (sum n x xsqr
)
1793 (math-working "log" sum
)
1794 (let* ((nextx (math-mul-float x xsqr
))
1795 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1796 (if (math-nearly-equal-float sum nextsum
)
1798 (math-ln-raw-series nextsum
(+ n
2) nextx xsqr
))))
1800 (defun math-ln-plus-1-raw (x)
1801 (math-lnp1-series x
2 x
(math-neg x
)))
1803 (defun math-lnp1-series (sum n xpow x
)
1804 (math-working "lnp1" sum
)
1805 (let* ((nextx (math-mul-float xpow x
))
1806 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1807 (if (math-nearly-equal-float sum nextsum
)
1809 (math-lnp1-series nextsum
(1+ n
) nextx x
))))
1811 (defconst math-approx-ln-10
1812 (math-read-number-simple "2.302585092994045684018")
1813 "An approximation for ln(10).")
1815 (math-defcache math-ln-10 math-approx-ln-10
1816 (math-ln-raw-2 '(float 1 1)))
1818 (defconst math-approx-ln-2
1819 (math-read-number-simple "0.693147180559945309417")
1820 "An approximation for ln(2).")
1822 (math-defcache math-ln-2 math-approx-ln-2
1823 (math-ln-raw-3 (math-float '(frac 1 3))))
1827 ;;; Hyperbolic functions.
1829 (defun calcFunc-sinh (x) ; [N N] [Public]
1831 (math-expand-formulas
1833 (list '/ (list '-
(list 'calcFunc-exp x
)
1834 (list 'calcFunc-exp
(list 'neg x
))) 2)))
1836 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1837 (math-with-extra-prec 2
1838 (let ((expx (math-exp-raw (math-float x
))))
1839 (math-mul (math-add expx
(math-div -
1 expx
)) '(float 5 -
1)))))
1840 ((eq (car-safe x
) 'sdev
)
1841 (math-make-sdev (calcFunc-sinh (nth 1 x
))
1842 (math-mul (nth 2 x
) (calcFunc-cosh (nth 1 x
)))))
1844 (math-sort-intv (nth 1 x
)
1845 (calcFunc-sinh (nth 2 x
))
1846 (calcFunc-sinh (nth 3 x
))))
1847 ((or (equal x
'(var inf var-inf
))
1848 (equal x
'(neg (var inf var-inf
)))
1849 (equal x
'(var nan var-nan
)))
1851 (t (calc-record-why 'numberp x
)
1852 (list 'calcFunc-sinh x
))))
1853 (put 'calcFunc-sinh
'math-expandable t
)
1855 (defun calcFunc-cosh (x) ; [N N] [Public]
1857 (math-expand-formulas
1859 (list '/ (list '+ (list 'calcFunc-exp x
)
1860 (list 'calcFunc-exp
(list 'neg x
))) 2)))
1862 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1863 (math-with-extra-prec 2
1864 (let ((expx (math-exp-raw (math-float x
))))
1865 (math-mul (math-add expx
(math-div 1 expx
)) '(float 5 -
1)))))
1866 ((eq (car-safe x
) 'sdev
)
1867 (math-make-sdev (calcFunc-cosh (nth 1 x
))
1869 (calcFunc-sinh (nth 1 x
)))))
1870 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1871 (setq x
(math-abs x
))
1872 (math-sort-intv (nth 1 x
)
1873 (calcFunc-cosh (nth 2 x
))
1874 (calcFunc-cosh (nth 3 x
))))
1875 ((or (equal x
'(var inf var-inf
))
1876 (equal x
'(neg (var inf var-inf
)))
1877 (equal x
'(var nan var-nan
)))
1879 (t (calc-record-why 'numberp x
)
1880 (list 'calcFunc-cosh x
))))
1881 (put 'calcFunc-cosh
'math-expandable t
)
1883 (defun calcFunc-tanh (x) ; [N N] [Public]
1885 (math-expand-formulas
1887 (let ((expx (list 'calcFunc-exp x
))
1888 (expmx (list 'calcFunc-exp
(list 'neg x
))))
1890 (list '/ (list '- expx expmx
) (list '+ expx expmx
))))))
1892 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1893 (math-with-extra-prec 2
1894 (let* ((expx (calcFunc-exp (math-float x
)))
1895 (expmx (math-div 1 expx
)))
1896 (math-div (math-sub expx expmx
)
1897 (math-add expx expmx
)))))
1898 ((eq (car-safe x
) 'sdev
)
1899 (math-make-sdev (calcFunc-tanh (nth 1 x
))
1901 (math-sqr (calcFunc-cosh (nth 1 x
))))))
1903 (math-sort-intv (nth 1 x
)
1904 (calcFunc-tanh (nth 2 x
))
1905 (calcFunc-tanh (nth 3 x
))))
1906 ((equal x
'(var inf var-inf
))
1908 ((equal x
'(neg (var inf var-inf
)))
1910 ((equal x
'(var nan var-nan
))
1912 (t (calc-record-why 'numberp x
)
1913 (list 'calcFunc-tanh x
))))
1914 (put 'calcFunc-tanh
'math-expandable t
)
1916 (defun calcFunc-sech (x) ; [N N] [Public]
1918 (math-expand-formulas
1920 (list '/ 2 (list '+ (list 'calcFunc-exp x
)
1921 (list 'calcFunc-exp
(list 'neg x
))))))
1923 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1924 (math-with-extra-prec 2
1925 (let ((expx (math-exp-raw (math-float x
))))
1926 (math-div '(float 2 0) (math-add expx
(math-div 1 expx
))))))
1927 ((eq (car-safe x
) 'sdev
)
1928 (math-make-sdev (calcFunc-sech (nth 1 x
))
1930 (math-mul (calcFunc-sech (nth 1 x
))
1931 (calcFunc-tanh (nth 1 x
))))))
1932 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1933 (setq x
(math-abs x
))
1934 (math-sort-intv (nth 1 x
)
1935 (calcFunc-sech (nth 2 x
))
1936 (calcFunc-sech (nth 3 x
))))
1937 ((or (equal x
'(var inf var-inf
))
1938 (equal x
'(neg (var inf var-inf
))))
1940 ((equal x
'(var nan var-nan
))
1942 (t (calc-record-why 'numberp x
)
1943 (list 'calcFunc-sech x
))))
1944 (put 'calcFunc-sech
'math-expandable t
)
1946 (defun calcFunc-csch (x) ; [N N] [Public]
1947 (cond ((eq x
0) (math-div 1 0))
1948 (math-expand-formulas
1950 (list '/ 2 (list '-
(list 'calcFunc-exp x
)
1951 (list 'calcFunc-exp
(list 'neg x
))))))
1953 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1954 (math-with-extra-prec 2
1955 (let ((expx (math-exp-raw (math-float x
))))
1956 (math-div '(float 2 0) (math-add expx
(math-div -
1 expx
))))))
1957 ((eq (car-safe x
) 'sdev
)
1958 (math-make-sdev (calcFunc-csch (nth 1 x
))
1960 (math-mul (calcFunc-csch (nth 1 x
))
1961 (calcFunc-coth (nth 1 x
))))))
1963 (if (and (Math-negp (nth 2 x
))
1964 (Math-posp (nth 3 x
)))
1965 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
1966 (math-sort-intv (nth 1 x
)
1967 (calcFunc-csch (nth 2 x
))
1968 (calcFunc-csch (nth 3 x
)))))
1969 ((or (equal x
'(var inf var-inf
))
1970 (equal x
'(neg (var inf var-inf
))))
1972 ((equal x
'(var nan var-nan
))
1974 (t (calc-record-why 'numberp x
)
1975 (list 'calcFunc-csch x
))))
1976 (put 'calcFunc-csch
'math-expandable t
)
1978 (defun calcFunc-coth (x) ; [N N] [Public]
1979 (cond ((eq x
0) (math-div 1 0))
1980 (math-expand-formulas
1982 (let ((expx (list 'calcFunc-exp x
))
1983 (expmx (list 'calcFunc-exp
(list 'neg x
))))
1985 (list '/ (list '+ expx expmx
) (list '- expx expmx
))))))
1987 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1988 (math-with-extra-prec 2
1989 (let* ((expx (calcFunc-exp (math-float x
)))
1990 (expmx (math-div 1 expx
)))
1991 (math-div (math-add expx expmx
)
1992 (math-sub expx expmx
)))))
1993 ((eq (car-safe x
) 'sdev
)
1994 (math-make-sdev (calcFunc-coth (nth 1 x
))
1996 (math-sqr (calcFunc-sinh (nth 1 x
))))))
1998 (if (and (Math-negp (nth 2 x
))
1999 (Math-posp (nth 3 x
)))
2000 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
2001 (math-sort-intv (nth 1 x
)
2002 (calcFunc-coth (nth 2 x
))
2003 (calcFunc-coth (nth 3 x
)))))
2004 ((equal x
'(var inf var-inf
))
2006 ((equal x
'(neg (var inf var-inf
)))
2008 ((equal x
'(var nan var-nan
))
2010 (t (calc-record-why 'numberp x
)
2011 (list 'calcFunc-coth x
))))
2012 (put 'calcFunc-coth
'math-expandable t
)
2014 (defun calcFunc-arcsinh (x) ; [N N] [Public]
2016 (math-expand-formulas
2018 (list 'calcFunc-ln
(list '+ x
(list 'calcFunc-sqrt
2019 (list '+ (list '^ x
2) 1))))))
2021 (if calc-symbolic-mode
(signal 'inexact-result nil
))
2022 (math-with-extra-prec 2
2023 (math-ln-raw (math-add x
(math-sqrt-raw (math-add (math-sqr x
)
2025 ((eq (car-safe x
) 'sdev
)
2026 (math-make-sdev (calcFunc-arcsinh (nth 1 x
))
2029 (math-add (math-sqr (nth 1 x
)) 1)))))
2031 (math-sort-intv (nth 1 x
)
2032 (calcFunc-arcsinh (nth 2 x
))
2033 (calcFunc-arcsinh (nth 3 x
))))
2034 ((or (equal x
'(var inf var-inf
))
2035 (equal x
'(neg (var inf var-inf
)))
2036 (equal x
'(var nan var-nan
)))
2038 (t (calc-record-why 'numberp x
)
2039 (list 'calcFunc-arcsinh x
))))
2040 (put 'calcFunc-arcsinh
'math-expandable t
)
2042 (defun calcFunc-arccosh (x) ; [N N] [Public]
2044 ((and (eq x -
1) calc-symbolic-mode
)
2046 ((and (eq x
0) calc-symbolic-mode
)
2047 (math-div (math-mul '(var pi var-pi
) '(var i var-i
)) 2))
2048 (math-expand-formulas
2050 (list 'calcFunc-ln
(list '+ x
(list 'calcFunc-sqrt
2051 (list '-
(list '^ x
2) 1))))))
2053 (if calc-symbolic-mode
(signal 'inexact-result nil
))
2054 (if (Math-equal-int x -
1)
2055 (math-imaginary (math-pi))
2056 (math-with-extra-prec 2
2057 (if (or t
; need to do this even in the real case!
2058 (memq (car-safe x
) '(cplx polar
)))
2059 (let ((xp1 (math-add 1 x
))) ; this gets the branch cuts right
2061 (math-add x
(math-mul xp1
2068 (math-add x
(math-sqrt-raw (math-add (math-sqr x
)
2069 '(float -
1 0)))))))))
2070 ((eq (car-safe x
) 'sdev
)
2071 (math-make-sdev (calcFunc-arccosh (nth 1 x
))
2074 (math-add (math-sqr (nth 1 x
)) -
1)))))
2076 (math-sort-intv (nth 1 x
)
2077 (calcFunc-arccosh (nth 2 x
))
2078 (calcFunc-arccosh (nth 3 x
))))
2079 ((or (equal x
'(var inf var-inf
))
2080 (equal x
'(neg (var inf var-inf
)))
2081 (equal x
'(var nan var-nan
)))
2083 (t (calc-record-why 'numberp x
)
2084 (list 'calcFunc-arccosh x
))))
2085 (put 'calcFunc-arccosh
'math-expandable t
)
2087 (defun calcFunc-arctanh (x) ; [N N] [Public]
2089 ((and (Math-equal-int x
1) calc-infinite-mode
)
2091 ((and (Math-equal-int x -
1) calc-infinite-mode
)
2092 '(neg (var inf var-inf
)))
2093 (math-expand-formulas
2095 (list 'calcFunc-ln
(list '+ 1 x
))
2096 (list 'calcFunc-ln
(list '-
1 x
))) 2))
2098 (if calc-symbolic-mode
(signal 'inexact-result nil
))
2099 (math-with-extra-prec 2
2100 (if (or (memq (car-safe x
) '(cplx polar
))
2102 (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x
))
2103 (math-ln-raw (math-sub '(float 1 0) x
)))
2105 (if (and (math-equal-int x
1) calc-infinite-mode
)
2107 (if (and (math-equal-int x -
1) calc-infinite-mode
)
2108 '(neg (var inf var-inf
))
2109 (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x
)
2112 ((eq (car-safe x
) 'sdev
)
2113 (math-make-sdev (calcFunc-arctanh (nth 1 x
))
2115 (math-sub 1 (math-sqr (nth 1 x
))))))
2117 (math-sort-intv (nth 1 x
)
2118 (calcFunc-arctanh (nth 2 x
))
2119 (calcFunc-arctanh (nth 3 x
))))
2120 ((equal x
'(var nan var-nan
))
2122 (t (calc-record-why 'numberp x
)
2123 (list 'calcFunc-arctanh x
))))
2124 (put 'calcFunc-arctanh
'math-expandable t
)
2127 ;;; Convert A from HMS or degrees to radians.
2128 (defun calcFunc-rad (a) ; [R R] [Public]
2129 (cond ((or (Math-numberp a
)
2131 (math-with-extra-prec 2
2132 (math-mul a
(math-pi-over-180))))
2134 (math-from-hms a
'rad
))
2136 (math-make-sdev (calcFunc-rad (nth 1 a
))
2137 (calcFunc-rad (nth 2 a
))))
2138 (math-expand-formulas
2139 (math-div (math-mul a
'(var pi var-pi
)) 180))
2140 ((math-infinitep a
) a
)
2141 (t (list 'calcFunc-rad a
))))
2142 (put 'calcFunc-rad
'math-expandable t
)
2144 ;;; Convert A from HMS or radians to degrees.
2145 (defun calcFunc-deg (a) ; [R R] [Public]
2146 (cond ((or (Math-numberp a
)
2148 (math-with-extra-prec 2
2149 (math-div a
(math-pi-over-180))))
2151 (math-from-hms a
'deg
))
2153 (math-make-sdev (calcFunc-deg (nth 1 a
))
2154 (calcFunc-deg (nth 2 a
))))
2155 (math-expand-formulas
2156 (math-div (math-mul 180 a
) '(var pi var-pi
)))
2157 ((math-infinitep a
) a
)
2158 (t (list 'calcFunc-deg a
))))
2159 (put 'calcFunc-deg
'math-expandable t
)
2161 (provide 'calc-math
)
2163 ;; arch-tag: c7367e8e-d0b8-4f70-8577-2fb3f31dbb4c
2164 ;;; calc-math.el ends here