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, 2010 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
)
110 ((calc-number-radix 10)
111 (calc-twos-complement-mode nil
)
112 (calc-float-format (list 'float calc-internal-prec
))
113 (calc-group-digits nil
)
114 (calc-point-char "."))
115 (math-format-number (math-float x
)))))))
118 (defun calc-sqrt (arg)
121 (if (calc-is-inverse)
122 (calc-unary-op "^2" 'calcFunc-sqr arg
)
123 (calc-unary-op "sqrt" 'calcFunc-sqrt arg
))))
125 (defun calc-isqrt (arg)
128 (if (calc-is-inverse)
129 (calc-unary-op "^2" 'calcFunc-sqr arg
)
130 (calc-unary-op "isqt" 'calcFunc-isqrt arg
))))
133 (defun calc-hypot (arg)
136 (calc-binary-op "hypt" 'calcFunc-hypot arg
)))
143 (defun calc-log10 (arg)
145 (calc-hyperbolic-func)
148 (defun calc-log (arg)
151 (if (calc-is-inverse)
152 (calc-binary-op "alog" 'calcFunc-alog arg
)
153 (calc-binary-op "log" 'calcFunc-log arg
))))
155 (defun calc-ilog (arg)
158 (if (calc-is-inverse)
159 (calc-binary-op "alog" 'calcFunc-alog arg
)
160 (calc-binary-op "ilog" 'calcFunc-ilog arg
))))
162 (defun calc-lnp1 (arg)
167 (defun calc-exp (arg)
170 (if (calc-is-hyperbolic)
171 (if (calc-is-inverse)
172 (calc-unary-op "lg10" 'calcFunc-log10 arg
)
173 (calc-unary-op "10^" 'calcFunc-exp10 arg
))
174 (if (calc-is-inverse)
175 (calc-unary-op "ln" 'calcFunc-ln arg
)
176 (calc-unary-op "exp" 'calcFunc-exp arg
)))))
178 (defun calc-expm1 (arg)
181 (if (calc-is-inverse)
182 (calc-unary-op "ln+1" 'calcFunc-lnp1 arg
)
183 (calc-unary-op "ex-1" 'calcFunc-expm1 arg
))))
188 (if (calc-is-inverse)
189 (if (calc-is-hyperbolic)
190 (if calc-symbolic-mode
191 (calc-pop-push-record 0 "phi" '(var phi var-phi
))
192 (calc-pop-push-record 0 "phi" (math-phi)))
193 (if calc-symbolic-mode
194 (calc-pop-push-record 0 "gmma" '(var gamma var-gamma
))
195 (calc-pop-push-record 0 "gmma" (math-gamma-const))))
196 (if (calc-is-hyperbolic)
197 (if calc-symbolic-mode
198 (calc-pop-push-record 0 "e" '(var e var-e
))
199 (calc-pop-push-record 0 "e" (math-e)))
200 (if calc-symbolic-mode
201 (calc-pop-push-record 0 "pi" '(var pi var-pi
))
202 (calc-pop-push-record 0 "pi" (math-pi)))))))
204 (defun calc-sin (arg)
207 (if (calc-is-hyperbolic)
208 (if (calc-is-inverse)
209 (calc-unary-op "asnh" 'calcFunc-arcsinh arg
)
210 (calc-unary-op "sinh" 'calcFunc-sinh arg
))
211 (if (calc-is-inverse)
212 (calc-unary-op "asin" 'calcFunc-arcsin arg
)
213 (calc-unary-op "sin" 'calcFunc-sin arg
)))))
215 (defun calc-arcsin (arg)
220 (defun calc-sinh (arg)
222 (calc-hyperbolic-func)
225 (defun calc-arcsinh (arg)
228 (calc-hyperbolic-func)
231 (defun calc-sec (arg)
234 (if (calc-is-hyperbolic)
235 (calc-unary-op "sech" 'calcFunc-sech arg
)
236 (calc-unary-op "sec" 'calcFunc-sec arg
))))
238 (defun calc-sech (arg)
240 (calc-hyperbolic-func)
243 (defun calc-cos (arg)
246 (if (calc-is-hyperbolic)
247 (if (calc-is-inverse)
248 (calc-unary-op "acsh" 'calcFunc-arccosh arg
)
249 (calc-unary-op "cosh" 'calcFunc-cosh arg
))
250 (if (calc-is-inverse)
251 (calc-unary-op "acos" 'calcFunc-arccos arg
)
252 (calc-unary-op "cos" 'calcFunc-cos arg
)))))
254 (defun calc-arccos (arg)
259 (defun calc-cosh (arg)
261 (calc-hyperbolic-func)
264 (defun calc-arccosh (arg)
267 (calc-hyperbolic-func)
270 (defun calc-csc (arg)
273 (if (calc-is-hyperbolic)
274 (calc-unary-op "csch" 'calcFunc-csch arg
)
275 (calc-unary-op "csc" 'calcFunc-csc arg
))))
277 (defun calc-csch (arg)
279 (calc-hyperbolic-func)
282 (defun calc-sincos ()
285 (if (calc-is-inverse)
286 (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos
(calc-top-n 1)))
287 (calc-enter-result 1 "sncs" (list 'calcFunc-sincos
(calc-top-n 1))))))
289 (defun calc-tan (arg)
292 (if (calc-is-hyperbolic)
293 (if (calc-is-inverse)
294 (calc-unary-op "atnh" 'calcFunc-arctanh arg
)
295 (calc-unary-op "tanh" 'calcFunc-tanh arg
))
296 (if (calc-is-inverse)
297 (calc-unary-op "atan" 'calcFunc-arctan arg
)
298 (calc-unary-op "tan" 'calcFunc-tan arg
)))))
300 (defun calc-arctan (arg)
305 (defun calc-tanh (arg)
307 (calc-hyperbolic-func)
310 (defun calc-arctanh (arg)
313 (calc-hyperbolic-func)
316 (defun calc-cot (arg)
319 (if (calc-is-hyperbolic)
320 (calc-unary-op "coth" 'calcFunc-coth arg
)
321 (calc-unary-op "cot" 'calcFunc-cot arg
))))
323 (defun calc-coth (arg)
325 (calc-hyperbolic-func)
328 (defun calc-arctan2 ()
331 (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2
(calc-top-list-n 2)))))
333 (defun calc-conj (arg)
336 (calc-unary-op "conj" 'calcFunc-conj arg
)))
338 (defun calc-imaginary ()
341 (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))))
343 (defun calc-to-degrees (arg)
346 (calc-unary-op ">deg" 'calcFunc-deg arg
)))
348 (defun calc-to-radians (arg)
351 (calc-unary-op ">rad" 'calcFunc-rad arg
)))
354 (defun calc-degrees-mode (arg)
358 (calc-change-mode 'calc-angle-mode
'deg
)
359 (message "Angles measured in degrees")))
360 ((= arg
2) (calc-radians-mode))
361 ((= arg
3) (calc-hms-mode))
362 (t (error "Prefix argument out of range"))))
364 (defun calc-radians-mode ()
367 (calc-change-mode 'calc-angle-mode
'rad
)
368 (message "Angles measured in radians")))
371 ;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public]
372 ;;; This method takes advantage of the fact that Newton's method starting
373 ;;; with an overestimate always works, even using truncating integer division!
374 (defun math-isqrt (a)
375 (cond ((Math-zerop a
) a
)
376 ((not (math-natnump a
))
377 (math-reject-arg a
'natnump
))
379 (math-isqrt-small a
))
381 (math-normalize (cons 'bigpos
(cdr (math-isqrt-bignum (cdr a
))))))))
383 (defun calcFunc-isqrt (a)
385 (math-isqrt (math-floor a
))
386 (math-floor (math-sqrt a
))))
389 ;;; This returns (flag . result) where the flag is t if A is a perfect square.
390 (defun math-isqrt-bignum (a) ; [P.l L]
391 (let ((len (length a
)))
393 (let* ((top (nthcdr (- len
2) a
)))
394 (math-isqrt-bignum-iter
396 (math-scale-bignum-digit-size
398 (1+ (math-isqrt-small
399 (+ (* (nth 1 top
) math-bignum-digit-size
) (car top
)))))
401 (let* ((top (nth (1- len
) a
)))
402 (math-isqrt-bignum-iter
404 (math-scale-bignum-digit-size
405 (list (1+ (math-isqrt-small top
)))
408 (defun math-isqrt-bignum-iter (a guess
) ; [l L l]
409 (math-working "isqrt" (cons 'bigpos guess
))
410 (let* ((q (math-div-bignum a guess
))
411 (s (math-add-bignum (car q
) guess
))
412 (g2 (math-div2-bignum s
))
413 (comp (math-compare-bignum g2 guess
)))
415 (math-isqrt-bignum-iter a g2
)
416 (cons (and (= comp
0)
417 (math-zerop-bignum (cdr q
))
421 (defun math-zerop-bignum (a)
424 (while (eq (car (setq a
(cdr a
))) 0))
427 (defun math-scale-bignum-digit-size (a n
) ; [L L S]
433 (defun math-isqrt-small (a) ; A > 0. [S S]
434 (let ((g (cond ((>= a
1000000) 10000)
439 (while (< (setq g2
(/ (+ g
(/ a g
)) 2)) g
)
446 ;;; Compute the square root of a number.
447 ;;; [T N] if possible, else [F N] if possible, else [C N]. [Public]
450 (and (Math-zerop a
) a
)
451 (and (math-known-nonposp a
)
452 (math-imaginary (math-sqrt (math-neg a
))))
454 (let ((sqrt (math-isqrt-small a
)))
455 (if (= (* sqrt sqrt
) a
)
457 (if calc-symbolic-mode
458 (list 'calcFunc-sqrt a
)
459 (math-sqrt-float (math-float a
) (math-float sqrt
))))))
460 (and (eq (car-safe a
) 'bigpos
)
461 (let* ((res (math-isqrt-bignum (cdr a
)))
462 (sqrt (math-normalize (cons 'bigpos
(cdr res
)))))
465 (if calc-symbolic-mode
466 (list 'calcFunc-sqrt a
)
467 (math-sqrt-float (math-float a
) (math-float sqrt
))))))
468 (and (eq (car-safe a
) 'frac
)
469 (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a
)))))
470 (num-sqrt (math-normalize (cons 'bigpos
(cdr num-res
))))
471 (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a
)))))
472 (den-sqrt (math-normalize (cons 'bigpos
(cdr den-res
)))))
473 (if (and (car num-res
) (car den-res
))
474 (list 'frac num-sqrt den-sqrt
)
475 (if calc-symbolic-mode
476 (if (or (car num-res
) (car den-res
))
477 (math-div (if (car num-res
)
478 num-sqrt
(list 'calcFunc-sqrt
(nth 1 a
)))
480 den-sqrt
(list 'calcFunc-sqrt
(nth 2 a
))))
481 (list 'calcFunc-sqrt a
))
482 (math-sqrt-float (math-float a
)
483 (math-div (math-float num-sqrt
) den-sqrt
))))))
484 (and (eq (car-safe a
) 'float
)
485 (if calc-symbolic-mode
486 (if (= (%
(nth 2 a
) 2) 0)
487 (let ((res (math-isqrt-bignum
488 (cdr (Math-bignum-test (nth 1 a
))))))
490 (math-make-float (math-normalize
491 (cons 'bigpos
(cdr res
)))
493 (signal 'inexact-result nil
)))
494 (signal 'inexact-result nil
))
495 (math-sqrt-float a
)))
496 (and (eq (car-safe a
) 'cplx
)
497 (math-with-extra-prec 2
498 (let* ((d (math-abs a
))
499 (imag (math-sqrt (math-mul (math-sub d
(nth 1 a
))
502 (math-sqrt (math-mul (math-add d
(nth 1 a
)) '(float 5 -
1)))
503 (if (math-negp (nth 2 a
)) (math-neg imag
) imag
)))))
504 (and (eq (car-safe a
) 'polar
)
506 (math-sqrt (nth 1 a
))
507 (math-mul (nth 2 a
) '(float 5 -
1))))
508 (and (eq (car-safe a
) 'sdev
)
509 (let ((sqrt (math-sqrt (nth 1 a
))))
511 (math-div (nth 2 a
) (math-mul sqrt
2)))))
512 (and (eq (car-safe a
) 'intv
)
513 (not (math-negp (nth 2 a
)))
514 (math-make-intv (nth 1 a
) (math-sqrt (nth 2 a
)) (math-sqrt (nth 3 a
))))
515 (and (eq (car-safe a
) '*)
516 (or (math-known-nonnegp (nth 1 a
))
517 (math-known-nonnegp (nth 2 a
)))
518 (math-mul (math-sqrt (nth 1 a
)) (math-sqrt (nth 2 a
))))
519 (and (eq (car-safe a
) '/)
520 (or (and (math-known-nonnegp (nth 2 a
))
521 (math-div (math-sqrt (nth 1 a
)) (math-sqrt (nth 2 a
))))
522 (and (math-known-nonnegp (nth 1 a
))
523 (not (math-equal-int (nth 1 a
) 1))
524 (math-mul (math-sqrt (nth 1 a
))
525 (math-sqrt (math-div 1 (nth 2 a
)))))))
526 (and (eq (car-safe a
) '^
)
527 (math-known-evenp (nth 2 a
))
528 (math-known-realp (nth 1 a
))
529 (math-abs (math-pow (nth 1 a
) (math-div (nth 2 a
) 2))))
530 (let ((inf (math-infinitep a
)))
532 (math-mul (math-sqrt (math-infinite-dir a inf
)) inf
)))
534 (calc-record-why 'numberp a
)
535 (list 'calcFunc-sqrt a
))))
536 (defalias 'calcFunc-sqrt
'math-sqrt
)
538 (defun math-infinite-dir (a &optional inf
)
539 (or inf
(setq inf
(math-infinitep a
)))
540 (math-normalize (math-expr-subst a inf
1)))
542 (defun math-sqrt-float (a &optional guess
) ; [F F F]
543 (if calc-symbolic-mode
544 (signal 'inexact-result nil
)
545 (math-with-extra-prec 1 (math-sqrt-raw a guess
))))
547 (defun math-sqrt-raw (a &optional guess
) ; [F F F]
548 (if (not (Math-posp a
))
551 ((math-use-emacs-fn 'sqrt a
))
554 (let ((ldiff (- (math-numdigs (nth 1 a
)) 6)))
555 (or (= (%
(+ (nth 2 a
) ldiff
) 2) 0) (setq ldiff
(1+ ldiff
)))
556 (setq guess
(math-make-float (math-isqrt-small
557 (math-scale-int (nth 1 a
) (- ldiff
)))
558 (/ (+ (nth 2 a
) ldiff
) 2)))))
559 (math-sqrt-float-iter a guess
)))))
561 (defun math-sqrt-float-iter (a guess
) ; [F F F]
562 (math-working "sqrt" guess
)
563 (let ((g2 (math-mul-float (math-add-float guess
(math-div-float a guess
))
565 (if (math-nearly-equal-float g2 guess
)
567 (math-sqrt-float-iter a g2
))))
569 ;;; True if A and B differ only in the last digit of precision. [P F F]
570 (defun math-nearly-equal-float (a b
)
571 (let ((ediff (- (nth 2 a
) (nth 2 b
))))
572 (cond ((= ediff
0) ;; Expanded out for speed
573 (setq ediff
(math-add (Math-integer-neg (nth 1 a
)) (nth 1 b
)))
575 (and (not (consp ediff
))
578 (= (math-numdigs (nth 1 a
)) calc-internal-prec
))))
580 (setq ediff
(math-add (Math-integer-neg (nth 1 b
))
581 (math-scale-int (nth 1 a
) 1)))
582 (and (not (consp ediff
))
585 (= (math-numdigs (nth 1 b
)) calc-internal-prec
)))
587 (setq ediff
(math-add (Math-integer-neg (nth 1 a
))
588 (math-scale-int (nth 1 b
) 1)))
589 (and (not (consp ediff
))
592 (= (math-numdigs (nth 1 a
)) calc-internal-prec
))))))
594 (defun math-nearly-equal (a b
) ; [P N N] [Public]
595 (setq a
(math-float a
))
596 (setq b
(math-float b
))
597 (if (eq (car a
) 'polar
) (setq a
(math-complex a
)))
598 (if (eq (car b
) 'polar
) (setq b
(math-complex b
)))
599 (if (eq (car a
) 'cplx
)
600 (if (eq (car b
) 'cplx
)
601 (and (or (math-nearly-equal-float (nth 1 a
) (nth 1 b
))
602 (and (math-nearly-zerop-float (nth 1 a
) (nth 2 a
))
603 (math-nearly-zerop-float (nth 1 b
) (nth 2 b
))))
604 (or (math-nearly-equal-float (nth 2 a
) (nth 2 b
))
605 (and (math-nearly-zerop-float (nth 2 a
) (nth 1 a
))
606 (math-nearly-zerop-float (nth 2 b
) (nth 1 b
)))))
607 (and (math-nearly-equal-float (nth 1 a
) b
)
608 (math-nearly-zerop-float (nth 2 a
) b
)))
609 (if (eq (car b
) 'cplx
)
610 (and (math-nearly-equal-float a
(nth 1 b
))
611 (math-nearly-zerop-float a
(nth 2 b
)))
612 (math-nearly-equal-float a b
))))
614 ;;; True if A is nearly zero compared to B. [P F F]
615 (defun math-nearly-zerop-float (a b
)
617 (<= (+ (math-numdigs (nth 1 a
)) (nth 2 a
))
618 (1+ (- (+ (math-numdigs (nth 1 b
)) (nth 2 b
)) calc-internal-prec
)))))
620 (defun math-nearly-zerop (a b
) ; [P N R] [Public]
621 (setq a
(math-float a
))
622 (setq b
(math-float b
))
623 (if (eq (car a
) 'cplx
)
624 (and (math-nearly-zerop-float (nth 1 a
) b
)
625 (math-nearly-zerop-float (nth 2 a
) b
))
626 (if (eq (car a
) 'polar
)
627 (math-nearly-zerop-float (nth 1 a
) b
)
628 (math-nearly-zerop-float a b
))))
630 ;;; This implementation could be improved, accuracy-wise.
631 (defun math-hypot (a b
)
632 (cond ((Math-zerop a
) (math-abs b
))
633 ((Math-zerop b
) (math-abs a
))
634 ((not (Math-scalarp a
))
635 (if (math-infinitep a
)
636 (if (math-infinitep b
)
641 (calc-record-why 'scalarp a
)
642 (list 'calcFunc-hypot a b
)))
643 ((not (Math-scalarp b
))
644 (if (math-infinitep b
)
646 (calc-record-why 'scalarp b
)
647 (list 'calcFunc-hypot a b
)))
648 ((and (Math-numberp a
) (Math-numberp b
))
649 (math-with-extra-prec 1
650 (math-sqrt (math-add (calcFunc-abssqr a
) (calcFunc-abssqr b
)))))
651 ((eq (car-safe a
) 'hms
)
652 (if (eq (car-safe b
) 'hms
) ; this helps sdev's of hms forms
653 (math-to-hms (math-hypot (math-from-hms a
'deg
)
654 (math-from-hms b
'deg
)))
655 (math-to-hms (math-hypot (math-from-hms a
'deg
) b
))))
656 ((eq (car-safe b
) 'hms
)
657 (math-to-hms (math-hypot a
(math-from-hms b
'deg
))))
659 (defalias 'calcFunc-hypot
'math-hypot
)
661 (defun calcFunc-sqr (x)
666 (defun math-nth-root (a n
)
667 (cond ((= n
2) (math-sqrt a
))
671 (let ((root (math-nth-root-integer a n
)))
674 (and (not calc-symbolic-mode
)
675 (math-nth-root-float (math-float a
) n
676 (math-float (cdr root
)))))))
677 ((eq (car-safe a
) 'frac
)
678 (let* ((num-root (math-nth-root-integer (nth 1 a
) n
))
679 (den-root (math-nth-root-integer (nth 2 a
) n
)))
680 (if (and (car num-root
) (car den-root
))
681 (list 'frac
(cdr num-root
) (cdr den-root
))
682 (and (not calc-symbolic-mode
)
685 (math-div-float (math-float (cdr num-root
))
686 (math-float (cdr den-root
))))))))
687 ((eq (car-safe a
) 'float
)
688 (and (not calc-symbolic-mode
)
689 (math-nth-root-float a n
)))
690 ((eq (car-safe a
) 'polar
)
691 (let ((root (math-nth-root (nth 1 a
) n
)))
692 (and root
(list 'polar root
(math-div (nth 2 a
) n
)))))
695 ;; The variables math-nrf-n, math-nrf-nf and math-nrf-nfm1 are local
696 ;; to math-nth-root-float, but are used by math-nth-root-float-iter,
697 ;; which is called by math-nth-root-float.
700 (defvar math-nrf-nfm1
)
702 (defun math-nth-root-float (a math-nrf-n
&optional guess
)
703 (math-inexact-result)
704 (math-with-extra-prec 1
705 (let ((math-nrf-nf (math-float math-nrf-n
))
706 (math-nrf-nfm1 (math-float (1- math-nrf-n
))))
707 (math-nth-root-float-iter a
(or guess
709 1 (/ (+ (math-numdigs (nth 1 a
))
714 (defun math-nth-root-float-iter (a guess
)
715 (math-working "root" guess
)
716 (let ((g2 (math-div-float (math-add-float (math-mul math-nrf-nfm1 guess
)
718 a
(math-ipow guess
(1- math-nrf-n
))))
720 (if (math-nearly-equal-float g2 guess
)
722 (math-nth-root-float-iter a g2
))))
724 ;; The variable math-nri-n is local to math-nth-root-integer, but
725 ;; is used by math-nth-root-int-iter, which is called by
726 ;; math-nth-root-int.
729 (defun math-nth-root-integer (a math-nri-n
&optional guess
) ; [I I S]
730 (math-nth-root-int-iter a
(or guess
731 (math-scale-int 1 (/ (+ (math-numdigs a
)
735 (defun math-nth-root-int-iter (a guess
)
736 (math-working "root" guess
)
737 (let* ((q (math-idivmod a
(math-ipow guess
(1- math-nri-n
))))
738 (s (math-add (car q
) (math-mul (1- math-nri-n
) guess
)))
739 (g2 (math-idivmod s math-nri-n
)))
740 (if (Math-natnum-lessp (car g2
) guess
)
741 (math-nth-root-int-iter a
(car g2
))
742 (cons (and (equal (car g2
) guess
)
747 (defun calcFunc-nroot (x n
)
748 (calcFunc-pow x
(if (integerp n
)
755 ;;;; Transcendental functions.
757 ;;; All of these functions are defined on the complex plane.
758 ;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
760 ;;; Most functions increase calc-internal-prec by 2 digits, then round
761 ;;; down afterward. "-raw" functions use the current precision, require
762 ;;; their arguments to be in float (or complex float) format, and always
763 ;;; work in radians (where applicable).
765 (defun math-to-radians (a) ; [N N]
766 (cond ((eq (car-safe a
) 'hms
)
767 (math-from-hms a
'rad
))
768 ((memq calc-angle-mode
'(deg hms
))
769 (math-mul a
(math-pi-over-180)))
772 (defun math-from-radians (a) ; [N N]
773 (cond ((eq calc-angle-mode
'deg
)
775 (math-div a
(math-pi-over-180))
776 (list 'calcFunc-deg a
)))
777 ((eq calc-angle-mode
'hms
)
778 (math-to-hms a
'rad
))
781 (defun math-to-radians-2 (a) ; [N N]
782 (cond ((eq (car-safe a
) 'hms
)
783 (math-from-hms a
'rad
))
784 ((memq calc-angle-mode
'(deg hms
))
785 (if calc-symbolic-mode
786 (math-div (math-mul a
'(var pi var-pi
)) 180)
787 (math-mul a
(math-pi-over-180))))
790 (defun math-from-radians-2 (a) ; [N N]
791 (cond ((memq calc-angle-mode
'(deg hms
))
792 (if calc-symbolic-mode
793 (math-div (math-mul 180 a
) '(var pi var-pi
))
794 (math-div a
(math-pi-over-180))))
799 ;;; Sine, cosine, and tangent.
801 (defun calcFunc-sin (x) ; [N N] [Public]
802 (cond ((and (integerp x
)
803 (if (eq calc-angle-mode
'deg
)
806 (aref [0 1 0 -
1] (math-mod (/ x
90) 4)))
808 (math-with-extra-prec 2
809 (math-sin-raw (math-to-radians (math-float x
)))))
812 (math-with-extra-prec 2
813 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
814 (xs (math-to-radians (math-float (nth 2 x
))))
815 (sc (math-sin-cos-raw xx
)))
816 (math-make-sdev (car sc
) (math-mul xs
(cdr sc
)))))
817 (math-make-sdev (calcFunc-sin (nth 1 x
))
818 (math-mul (nth 2 x
) (calcFunc-cos (nth 1 x
))))))
819 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
820 (calcFunc-cos (math-sub x
(math-quarter-circle nil
))))
821 ((equal x
'(var nan var-nan
))
823 (t (calc-record-why 'scalarp x
)
824 (list 'calcFunc-sin x
))))
826 (defun calcFunc-cos (x) ; [N N] [Public]
827 (cond ((and (integerp x
)
828 (if (eq calc-angle-mode
'deg
)
831 (aref [1 0 -
1 0] (math-mod (/ x
90) 4)))
833 (math-with-extra-prec 2
834 (math-cos-raw (math-to-radians (math-float x
)))))
837 (math-with-extra-prec 2
838 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
839 (xs (math-to-radians (math-float (nth 2 x
))))
840 (sc (math-sin-cos-raw xx
)))
841 (math-make-sdev (cdr sc
) (math-mul xs
(car sc
)))))
842 (math-make-sdev (calcFunc-cos (nth 1 x
))
843 (math-mul (nth 2 x
) (calcFunc-sin (nth 1 x
))))))
844 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
845 (math-with-extra-prec 2
846 (let* ((xx (math-to-radians (math-float x
)))
847 (na (math-floor (math-div (nth 2 xx
) (math-pi))))
848 (nb (math-floor (math-div (nth 3 xx
) (math-pi))))
849 (span (math-sub nb na
)))
850 (if (memq span
'(0 1))
851 (let ((int (math-sort-intv (nth 1 x
)
852 (math-cos-raw (nth 2 xx
))
853 (math-cos-raw (nth 3 xx
)))))
856 (math-make-intv (logior (nth 1 x
) 2)
859 (math-make-intv (logior (nth 1 x
) 1)
863 (list 'intv
3 -
1 1)))))
864 ((equal x
'(var nan var-nan
))
866 (t (calc-record-why 'scalarp x
)
867 (list 'calcFunc-cos x
))))
869 (defun calcFunc-sincos (x) ; [V N] [Public]
871 (math-with-extra-prec 2
872 (let ((sc (math-sin-cos-raw (math-to-radians (math-float x
)))))
873 (list 'vec
(cdr sc
) (car sc
)))) ; the vector [cos, sin]
874 (list 'vec
(calcFunc-sin x
) (calcFunc-cos x
))))
876 (defun calcFunc-tan (x) ; [N N] [Public]
877 (cond ((and (integerp x
)
878 (if (eq calc-angle-mode
'deg
)
883 (math-with-extra-prec 2
884 (math-tan-raw (math-to-radians (math-float x
)))))
887 (math-with-extra-prec 2
888 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
889 (xs (math-to-radians (math-float (nth 2 x
))))
890 (sc (math-sin-cos-raw xx
)))
891 (if (and (math-zerop (cdr sc
)) (not calc-infinite-mode
))
893 (calc-record-why "*Division by zero")
894 (list 'calcFunc-tan x
))
895 (math-make-sdev (math-div-float (car sc
) (cdr sc
))
896 (math-div-float xs
(math-sqr (cdr sc
)))))))
897 (math-make-sdev (calcFunc-tan (nth 1 x
))
899 (math-sqr (calcFunc-cos (nth 1 x
)))))))
900 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
901 (or (math-with-extra-prec 2
902 (let* ((xx (math-to-radians (math-float x
)))
903 (na (math-floor (math-div (math-sub (nth 2 xx
)
906 (nb (math-floor (math-div (math-sub (nth 3 xx
)
910 (math-sort-intv (nth 1 x
)
911 (math-tan-raw (nth 2 xx
))
912 (math-tan-raw (nth 3 xx
))))))
913 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))))
914 ((equal x
'(var nan var-nan
))
916 (t (calc-record-why 'scalarp x
)
917 (list 'calcFunc-tan x
))))
919 (defun calcFunc-sec (x)
920 (cond ((and (integerp x
)
921 (eq calc-angle-mode
'deg
)
927 (eq calc-angle-mode
'rad
)
931 (math-with-extra-prec 2
932 (math-sec-raw (math-to-radians (math-float x
)))))
935 (math-with-extra-prec 2
936 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
937 (xs (math-to-radians (math-float (nth 2 x
))))
938 (sc (math-sin-cos-raw xx
)))
939 (if (and (math-zerop (cdr sc
))
940 (not calc-infinite-mode
))
942 (calc-record-why "*Division by zero")
943 (list 'calcFunc-sec x
))
944 (math-make-sdev (math-div-float '(float 1 0) (cdr sc
))
946 (math-mul xs
(car sc
))
947 (math-sqr (cdr sc
)))))))
948 (math-make-sdev (calcFunc-sec (nth 1 x
))
951 (calcFunc-sin (nth 1 x
)))
952 (math-sqr (calcFunc-cos (nth 1 x
)))))))
953 ((and (eq (car x
) 'intv
)
954 (math-intv-constp x
))
955 (math-with-extra-prec 2
956 (let* ((xx (math-to-radians (math-float x
)))
957 (na (math-floor (math-div (math-sub (nth 2 xx
)
960 (nb (math-floor (math-div (math-sub (nth 3 xx
)
963 (naa (math-floor (math-div (nth 2 xx
) (math-pi-over-2))))
964 (nbb (math-floor (math-div (nth 3 xx
) (math-pi-over-2))))
965 (span (math-sub nbb naa
)))
966 (if (not (equal na nb
))
967 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
968 (let ((int (math-sort-intv (nth 1 x
)
969 (math-sec-raw (nth 2 xx
))
970 (math-sec-raw (nth 3 xx
)))))
972 (if (math-evenp (math-div (math-add naa
1) 2))
973 (math-make-intv (logior (nth 1 int
) 2)
976 (math-make-intv (logior (nth 1 int
) 1)
980 ((equal x
'(var nan var-nan
))
982 (t (calc-record-why 'scalarp x
)
983 (list 'calcFunc-sec x
))))
985 (defun calcFunc-csc (x)
986 (cond ((and (integerp x
)
987 (eq calc-angle-mode
'deg
)
988 (= (%
(- x
90) 180) 0))
989 (if (= (%
(- x
90) 360) 0)
993 (math-with-extra-prec 2
994 (math-csc-raw (math-to-radians (math-float x
)))))
997 (math-with-extra-prec 2
998 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
999 (xs (math-to-radians (math-float (nth 2 x
))))
1000 (sc (math-sin-cos-raw xx
)))
1001 (if (and (math-zerop (car sc
))
1002 (not calc-infinite-mode
))
1004 (calc-record-why "*Division by zero")
1005 (list 'calcFunc-csc x
))
1006 (math-make-sdev (math-div-float '(float 1 0) (car sc
))
1008 (math-mul xs
(cdr sc
))
1009 (math-sqr (car sc
)))))))
1010 (math-make-sdev (calcFunc-csc (nth 1 x
))
1013 (calcFunc-cos (nth 1 x
)))
1014 (math-sqr (calcFunc-sin (nth 1 x
)))))))
1015 ((and (eq (car x
) 'intv
)
1016 (math-intv-constp x
))
1017 (math-with-extra-prec 2
1018 (let* ((xx (math-to-radians (math-float x
)))
1019 (na (math-floor (math-div (nth 2 xx
) (math-pi))))
1020 (nb (math-floor (math-div (nth 3 xx
) (math-pi))))
1021 (naa (math-floor (math-div (nth 2 xx
) (math-pi-over-2))))
1022 (nbb (math-floor (math-div (nth 3 xx
) (math-pi-over-2))))
1023 (span (math-sub nbb naa
)))
1024 (if (not (equal na nb
))
1025 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
1026 (let ((int (math-sort-intv (nth 1 x
)
1027 (math-csc-raw (nth 2 xx
))
1028 (math-csc-raw (nth 3 xx
)))))
1030 (if (math-evenp (math-div naa
2))
1031 (math-make-intv (logior (nth 1 int
) 2)
1034 (math-make-intv (logior (nth 1 int
) 1)
1038 ((equal x
'(var nan var-nan
))
1040 (t (calc-record-why 'scalarp x
)
1041 (list 'calcFunc-csc x
))))
1043 (defun calcFunc-cot (x) ; [N N] [Public]
1044 (cond ((and (integerp x
)
1045 (if (eq calc-angle-mode
'deg
)
1046 (= (%
(- x
90) 180) 0)
1050 (math-with-extra-prec 2
1051 (math-cot-raw (math-to-radians (math-float x
)))))
1054 (math-with-extra-prec 2
1055 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
1056 (xs (math-to-radians (math-float (nth 2 x
))))
1057 (sc (math-sin-cos-raw xx
)))
1058 (if (and (math-zerop (car sc
)) (not calc-infinite-mode
))
1060 (calc-record-why "*Division by zero")
1061 (list 'calcFunc-cot x
))
1062 (math-make-sdev (math-div-float (cdr sc
) (car sc
))
1063 (math-div-float xs
(math-sqr (car sc
)))))))
1064 (math-make-sdev (calcFunc-cot (nth 1 x
))
1066 (math-sqr (calcFunc-sin (nth 1 x
)))))))
1067 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1068 (or (math-with-extra-prec 2
1069 (let* ((xx (math-to-radians (math-float x
)))
1070 (na (math-floor (math-div (nth 2 xx
) (math-pi))))
1071 (nb (math-floor (math-div (nth 3 xx
) (math-pi)))))
1073 (math-sort-intv (nth 1 x
)
1074 (math-cot-raw (nth 2 xx
))
1075 (math-cot-raw (nth 3 xx
))))))
1076 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))))
1077 ((equal x
'(var nan var-nan
))
1079 (t (calc-record-why 'scalarp x
)
1080 (list 'calcFunc-cot x
))))
1082 (defun math-sin-raw (x &optional orgx
) ; [N N]
1083 (cond ((eq (car x
) 'cplx
)
1084 (let* ((expx (math-exp-raw (nth 2 x
)))
1085 (expmx (math-div-float '(float 1 0) expx
))
1086 (sc (math-sin-cos-raw (nth 1 x
))))
1088 (math-mul-float (car sc
)
1089 (math-mul-float (math-add-float expx expmx
)
1091 (math-mul-float (cdr sc
)
1092 (math-mul-float (math-sub-float expx expmx
)
1094 ((eq (car x
) 'polar
)
1095 (math-polar (math-sin-raw (math-complex x
))))
1096 ((Math-integer-negp (nth 1 x
))
1097 (math-neg-float (math-sin-raw (math-neg-float x
) (if orgx orgx x
))))
1098 ((math-lessp-float '(float 7 0) x
) ; avoid inf loops due to roundoff
1099 (math-sin-raw (math-mod x
(math-two-pi)) (if orgx orgx x
)))
1100 (t (math-sin-raw-2 x
(if orgx orgx x
)))))
1102 (defun math-cos-raw (x) ; [N N]
1103 (if (eq (car-safe x
) 'polar
)
1104 (math-polar (math-cos-raw (math-complex x
)))
1105 (math-sin-raw (math-sub (math-pi-over-2) x
) x
)))
1107 (defun math-sec-raw (x) ; [N N]
1108 (cond ((eq (car x
) 'cplx
)
1109 (let* ((x (math-mul x
'(float 1 0)))
1110 (expx (math-exp-raw (nth 2 x
)))
1111 (expmx (math-div-float '(float 1 0) expx
))
1112 (sh (math-mul-float (math-sub-float expx expmx
) '(float 5 -
1)))
1113 (ch (math-mul-float (math-add-float expx expmx
) '(float 5 -
1)))
1114 (sc (math-sin-cos-raw (nth 1 x
)))
1116 (math-mul-float (math-sqr (car sc
))
1118 (math-mul-float (math-sqr (cdr sc
))
1120 (and (not (eq (nth 1 d
) 0))
1122 (math-div-float (math-mul-float (cdr sc
) ch
) d
)
1123 (math-div-float (math-mul-float (car sc
) sh
) d
)))))
1124 ((eq (car x
) 'polar
)
1125 (math-polar (math-sec-raw (math-complex x
))))
1127 (let ((cs (math-cos-raw x
)))
1130 (math-div-float '(float 1 0) cs
))))))
1132 (defun math-csc-raw (x) ; [N N]
1133 (cond ((eq (car x
) 'cplx
)
1134 (let* ((x (math-mul x
'(float 1 0)))
1135 (expx (math-exp-raw (nth 2 x
)))
1136 (expmx (math-div-float '(float 1 0) expx
))
1137 (sh (math-mul-float (math-sub-float expx expmx
) '(float 5 -
1)))
1138 (ch (math-mul-float (math-add-float expx expmx
) '(float 5 -
1)))
1139 (sc (math-sin-cos-raw (nth 1 x
)))
1141 (math-mul-float (math-sqr (car sc
))
1143 (math-mul-float (math-sqr (cdr sc
))
1145 (and (not (eq (nth 1 d
) 0))
1147 (math-div-float (math-mul-float (car sc
) ch
) d
)
1148 (math-div-float (math-mul-float (cdr sc
) sh
) d
)))))
1149 ((eq (car x
) 'polar
)
1150 (math-polar (math-csc-raw (math-complex x
))))
1152 (let ((sn (math-sin-raw x
)))
1155 (math-div-float '(float 1 0) sn
))))))
1157 (defun math-cot-raw (x) ; [N N]
1158 (cond ((eq (car x
) 'cplx
)
1159 (let* ((x (math-mul x
'(float 1 0)))
1160 (expx (math-exp-raw (nth 2 x
)))
1161 (expmx (math-div-float '(float 1 0) expx
))
1162 (sh (math-mul-float (math-sub-float expx expmx
) '(float 5 -
1)))
1163 (ch (math-mul-float (math-add-float expx expmx
) '(float 5 -
1)))
1164 (sc (math-sin-cos-raw (nth 1 x
)))
1168 (and (not (eq (nth 1 d
) 0))
1171 (math-mul-float (car sc
) (cdr sc
))
1175 (math-mul-float sh ch
)
1177 ((eq (car x
) 'polar
)
1178 (math-polar (math-cot-raw (math-complex x
))))
1180 (let ((sc (math-sin-cos-raw x
)))
1181 (if (eq (nth 1 (car sc
)) 0)
1182 (math-div (cdr sc
) 0)
1183 (math-div-float (cdr sc
) (car sc
)))))))
1186 ;;; This could use a smarter method: Reduce x as in math-sin-raw, then
1187 ;;; compute either sin(x) or cos(x), whichever is smaller, and compute
1188 ;;; the other using the identity sin(x)^2 + cos(x)^2 = 1.
1189 (defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x))
1190 (cons (math-sin-raw x
) (math-cos-raw x
)))
1192 (defun math-tan-raw (x) ; [N N]
1193 (cond ((eq (car x
) 'cplx
)
1194 (let* ((x (math-mul x
'(float 2 0)))
1195 (expx (math-exp-raw (nth 2 x
)))
1196 (expmx (math-div-float '(float 1 0) expx
))
1197 (sc (math-sin-cos-raw (nth 1 x
)))
1198 (d (math-add-float (cdr sc
)
1199 (math-mul-float (math-add-float expx expmx
)
1201 (and (not (eq (nth 1 d
) 0))
1203 (math-div-float (car sc
) d
)
1204 (math-div-float (math-mul-float (math-sub-float expx
1206 '(float 5 -
1)) d
)))))
1207 ((eq (car x
) 'polar
)
1208 (math-polar (math-tan-raw (math-complex x
))))
1210 (let ((sc (math-sin-cos-raw x
)))
1211 (if (eq (nth 1 (cdr sc
)) 0)
1212 (math-div (car sc
) 0)
1213 (math-div-float (car sc
) (cdr sc
)))))))
1215 (defun math-sin-raw-2 (x orgx
) ; This avoids poss of inf recursion. [F F]
1216 (let ((xmpo2 (math-sub-float (math-pi-over-2) x
)))
1217 (cond ((Math-integer-negp (nth 1 xmpo2
))
1218 (math-neg-float (math-sin-raw-2 (math-sub-float x
(math-pi))
1220 ((math-lessp-float (math-pi-over-4) x
)
1221 (math-cos-raw-2 xmpo2 orgx
))
1222 ((math-lessp-float x
(math-neg (math-pi-over-4)))
1223 (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x
) orgx
)))
1224 ((math-with-extra-prec -
1 (math-nearly-zerop-float x orgx
))
1226 ((math-use-emacs-fn 'sin x
))
1227 (calc-symbolic-mode (signal 'inexact-result nil
))
1228 (t (math-sin-series x
6 4 x
(math-neg-float (math-sqr-float x
)))))))
1230 (defun math-cos-raw-2 (x orgx
) ; [F F]
1231 (cond ((math-with-extra-prec -
1 (math-nearly-zerop-float x orgx
))
1233 ((math-use-emacs-fn 'cos x
))
1234 (calc-symbolic-mode (signal 'inexact-result nil
))
1235 (t (let ((xnegsqr (math-neg-float (math-sqr-float x
))))
1237 (math-add-float '(float 1 0)
1238 (math-mul-float xnegsqr
'(float 5 -
1)))
1239 24 5 xnegsqr xnegsqr
)))))
1241 (defun math-sin-series (sum nfac n x xnegsqr
)
1242 (math-working "sin" sum
)
1243 (let* ((nextx (math-mul-float x xnegsqr
))
1244 (nextsum (math-add-float sum
(math-div-float nextx
1245 (math-float nfac
)))))
1246 (if (math-nearly-equal-float sum nextsum
)
1248 (math-sin-series nextsum
(math-mul nfac
(* n
(1+ n
)))
1249 (+ n
2) nextx xnegsqr
))))
1252 ;;; Inverse sine, cosine, tangent.
1254 (defun calcFunc-arcsin (x) ; [N N] [Public]
1256 ((and (eq x
1) (eq calc-angle-mode
'deg
)) 90)
1257 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) -
90)
1258 (calc-symbolic-mode (signal 'inexact-result nil
))
1260 (math-with-extra-prec 2
1261 (math-from-radians (math-arcsin-raw (math-float x
)))))
1263 (math-make-sdev (calcFunc-arcsin (nth 1 x
))
1267 (math-sub 1 (math-sqr (nth 1 x
))))))))
1269 (math-sort-intv (nth 1 x
)
1270 (calcFunc-arcsin (nth 2 x
))
1271 (calcFunc-arcsin (nth 3 x
))))
1272 ((equal x
'(var nan var-nan
))
1274 (t (calc-record-why 'numberp x
)
1275 (list 'calcFunc-arcsin x
))))
1277 (defun calcFunc-arccos (x) ; [N N] [Public]
1279 ((and (eq x
0) (eq calc-angle-mode
'deg
)) 90)
1280 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) 180)
1281 (calc-symbolic-mode (signal 'inexact-result nil
))
1283 (math-with-extra-prec 2
1284 (math-from-radians (math-arccos-raw (math-float x
)))))
1286 (math-make-sdev (calcFunc-arccos (nth 1 x
))
1290 (math-sub 1 (math-sqr (nth 1 x
))))))))
1292 (math-sort-intv (nth 1 x
)
1293 (calcFunc-arccos (nth 2 x
))
1294 (calcFunc-arccos (nth 3 x
))))
1295 ((equal x
'(var nan var-nan
))
1297 (t (calc-record-why 'numberp x
)
1298 (list 'calcFunc-arccos x
))))
1300 (defun calcFunc-arctan (x) ; [N N] [Public]
1302 ((and (eq x
1) (eq calc-angle-mode
'deg
)) 45)
1303 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) -
45)
1305 (math-with-extra-prec 2
1306 (math-from-radians (math-arctan-raw (math-float x
)))))
1308 (math-make-sdev (calcFunc-arctan (nth 1 x
))
1311 (math-add 1 (math-sqr (nth 1 x
)))))))
1313 (math-sort-intv (nth 1 x
)
1314 (calcFunc-arctan (nth 2 x
))
1315 (calcFunc-arctan (nth 3 x
))))
1316 ((equal x
'(var inf var-inf
))
1317 (math-quarter-circle t
))
1318 ((equal x
'(neg (var inf var-inf
)))
1319 (math-neg (math-quarter-circle t
)))
1320 ((equal x
'(var nan var-nan
))
1322 (t (calc-record-why 'numberp x
)
1323 (list 'calcFunc-arctan x
))))
1325 (defun math-arcsin-raw (x) ; [N N]
1326 (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x
)))))
1327 (if (or (memq (car x
) '(cplx polar
))
1328 (memq (car a
) '(cplx polar
)))
1329 (math-with-extra-prec 2 ; use extra precision for difficult case
1330 (math-mul '(cplx 0 -
1)
1331 (math-ln-raw (math-add (math-mul '(cplx 0 1) x
) a
))))
1332 (math-arctan2-raw x a
))))
1334 (defun math-arccos-raw (x) ; [N N]
1335 (math-sub (math-pi-over-2) (math-arcsin-raw x
)))
1337 (defun math-arctan-raw (x) ; [N N]
1338 (cond ((memq (car x
) '(cplx polar
))
1339 (math-with-extra-prec 2 ; extra-extra
1341 (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x
)))
1342 (math-ln-raw (math-add 1 (math-mul '(cplx 0 -
1) x
))))
1344 ((Math-integer-negp (nth 1 x
))
1345 (math-neg-float (math-arctan-raw (math-neg-float x
))))
1347 ((math-use-emacs-fn 'atan x
))
1348 (calc-symbolic-mode (signal 'inexact-result nil
))
1349 ((math-equal-int x
1) (math-pi-over-4))
1350 ((math-equal-int x -
1) (math-neg (math-pi-over-4)))
1351 ((math-lessp-float '(float 414214 -
6) x
) ; if x > sqrt(2) - 1, reduce
1352 (if (math-lessp-float '(float 1 0) x
)
1353 (math-sub-float (math-mul-float (math-pi) '(float 5 -
1))
1354 (math-arctan-raw (math-div-float '(float 1 0) x
)))
1355 (math-sub-float (math-mul-float (math-pi) '(float 25 -
2))
1356 (math-arctan-raw (math-div-float
1357 (math-sub-float '(float 1 0) x
)
1358 (math-add-float '(float 1 0)
1360 (t (math-arctan-series x
3 x
(math-neg-float (math-sqr-float x
))))))
1362 (defun math-arctan-series (sum n x xnegsqr
)
1363 (math-working "arctan" sum
)
1364 (let* ((nextx (math-mul-float x xnegsqr
))
1365 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1366 (if (math-nearly-equal-float sum nextsum
)
1368 (math-arctan-series nextsum
(+ n
2) nextx xnegsqr
))))
1370 (defun calcFunc-arctan2 (y x
) ; [F R R] [Public]
1373 (math-with-extra-prec 2
1374 (math-from-radians (math-arctan2-raw (math-float y
)
1376 (calc-record-why 'anglep x
)
1377 (list 'calcFunc-arctan2 y x
))
1378 (if (and (or (math-infinitep x
) (math-anglep x
))
1379 (or (math-infinitep y
) (math-anglep y
)))
1394 (calcFunc-arctan2 y x
)
1395 '(var nan var-nan
)))
1396 (calc-record-why 'anglep y
)
1397 (list 'calcFunc-arctan2 y x
))))
1399 (defun math-arctan2-raw (y x
) ; [F R R]
1400 (cond ((math-zerop y
)
1401 (if (math-negp x
) (math-pi)
1402 (if (or (math-floatp x
) (math-floatp y
)) '(float 0 0) 0)))
1406 (math-neg (math-pi-over-2))))
1408 (math-arctan-raw (math-div-float y x
)))
1410 (math-add-float (math-arctan-raw (math-div-float y x
))
1413 (math-sub-float (math-arctan-raw (math-div-float y x
))
1416 (defun calcFunc-arcsincos (x) ; [V N] [Public]
1417 (if (and (Math-vectorp x
)
1419 (calcFunc-arctan2 (nth 2 x
) (nth 1 x
))
1420 (math-reject-arg x
"*Two-element vector expected")))
1424 ;;; Exponential function.
1426 (defun calcFunc-exp (x) ; [N N] [Public]
1428 ((and (memq x
'(1 -
1)) calc-symbolic-mode
)
1429 (if (eq x
1) '(var e var-e
) (math-div 1 '(var e var-e
))))
1431 (math-with-extra-prec 2 (math-exp-raw (math-float x
))))
1432 ((eq (car-safe x
) 'sdev
)
1433 (let ((ex (calcFunc-exp (nth 1 x
))))
1434 (math-make-sdev ex
(math-mul (nth 2 x
) ex
))))
1435 ((eq (car-safe x
) 'intv
)
1436 (math-make-intv (nth 1 x
) (calcFunc-exp (nth 2 x
))
1437 (calcFunc-exp (nth 3 x
))))
1438 ((equal x
'(var inf var-inf
))
1440 ((equal x
'(neg (var inf var-inf
)))
1442 ((equal x
'(var nan var-nan
))
1444 (t (calc-record-why 'numberp x
)
1445 (list 'calcFunc-exp x
))))
1447 (defun calcFunc-expm1 (x) ; [N N] [Public]
1449 ((math-zerop x
) '(float 0 0))
1450 (calc-symbolic-mode (signal 'inexact-result nil
))
1452 (math-with-extra-prec 2
1453 (let ((x (math-float x
)))
1454 (if (and (eq (car x
) 'float
)
1455 (math-lessp-float x
'(float 1 0))
1456 (math-lessp-float '(float -
1 0) x
))
1457 (math-exp-minus-1-raw x
)
1458 (math-add (math-exp-raw x
) -
1)))))
1459 ((eq (car-safe x
) 'sdev
)
1461 (let ((ex (calcFunc-expm1 (nth 1 x
))))
1462 (math-make-sdev ex
(math-mul (nth 2 x
) (math-add ex
1))))
1463 (math-make-sdev (calcFunc-expm1 (nth 1 x
))
1464 (math-mul (nth 2 x
) (calcFunc-exp (nth 1 x
))))))
1465 ((eq (car-safe x
) 'intv
)
1466 (math-make-intv (nth 1 x
)
1467 (calcFunc-expm1 (nth 2 x
))
1468 (calcFunc-expm1 (nth 3 x
))))
1469 ((equal x
'(var inf var-inf
))
1471 ((equal x
'(neg (var inf var-inf
)))
1473 ((equal x
'(var nan var-nan
))
1475 (t (calc-record-why 'numberp x
)
1476 (list 'calcFunc-expm1 x
))))
1478 (defun calcFunc-exp10 (x) ; [N N] [Public]
1481 (math-pow '(float 1 1) x
)))
1483 (defun math-exp-raw (x) ; [N N]
1484 (cond ((math-zerop x
) '(float 1 0))
1485 (calc-symbolic-mode (signal 'inexact-result nil
))
1487 (let ((expx (math-exp-raw (nth 1 x
)))
1488 (sc (math-sin-cos-raw (nth 2 x
))))
1490 (math-mul-float expx
(cdr sc
))
1491 (math-mul-float expx
(car sc
)))))
1492 ((eq (car x
) 'polar
)
1493 (let ((xc (math-complex x
)))
1495 (math-exp-raw (nth 1 xc
))
1496 (math-from-radians (nth 2 xc
)))))
1497 ((math-use-emacs-fn 'exp x
))
1498 ((or (math-lessp-float '(float 5 -
1) x
)
1499 (math-lessp-float x
'(float -
5 -
1)))
1500 (if (math-lessp-float '(float 921035 1) x
)
1502 (if (math-lessp-float x
'(float -
921035 1))
1504 (let* ((two-x (math-mul-float x
'(float 2 0)))
1505 (hint (math-scale-int (nth 1 two-x
) (nth 2 two-x
)))
1506 (hfrac (math-sub-float x
(math-mul-float (math-float hint
)
1508 (math-mul-float (math-ipow (math-sqrt-e) hint
)
1509 (math-add-float '(float 1 0)
1510 (math-exp-minus-1-raw hfrac
)))))
1511 (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x
)))))
1513 (defun math-exp-minus-1-raw (x) ; [F F]
1514 (math-exp-series x
2 3 x x
))
1516 (defun math-exp-series (sum nfac n xpow x
)
1517 (math-working "exp" sum
)
1518 (let* ((nextx (math-mul-float xpow x
))
1519 (nextsum (math-add-float sum
(math-div-float nextx
1520 (math-float nfac
)))))
1521 (if (math-nearly-equal-float sum nextsum
)
1523 (math-exp-series nextsum
(math-mul nfac n
) (1+ n
) nextx x
))))
1529 (defun calcFunc-ln (x) ; [N N] [Public]
1530 (cond ((math-zerop x
)
1531 (if calc-infinite-mode
1532 '(neg (var inf var-inf
))
1533 (math-reject-arg x
"*Logarithm of zero")))
1536 (math-with-extra-prec 2 (math-ln-raw (math-float x
))))
1537 ((eq (car-safe x
) 'sdev
)
1538 (math-make-sdev (calcFunc-ln (nth 1 x
))
1539 (math-div (nth 2 x
) (nth 1 x
))))
1540 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1541 (Math-zerop (nth 2 x
))
1542 (not (math-intv-constp x
))))
1543 (let ((calc-infinite-mode t
))
1544 (math-make-intv (nth 1 x
) (calcFunc-ln (nth 2 x
))
1545 (calcFunc-ln (nth 3 x
)))))
1546 ((equal x
'(var e var-e
))
1548 ((and (eq (car-safe x
) '^
)
1549 (equal (nth 1 x
) '(var e var-e
))
1550 (math-known-realp (nth 2 x
)))
1553 (if (equal x
'(var nan var-nan
))
1555 '(var inf var-inf
)))
1556 (t (calc-record-why 'numberp x
)
1557 (list 'calcFunc-ln x
))))
1559 (defun calcFunc-log10 (x) ; [N N] [Public]
1560 (cond ((math-equal-int x
1)
1561 (if (math-floatp x
) '(float 0 0) 0))
1562 ((and (Math-integerp x
)
1564 (let ((res (math-integer-log x
10)))
1566 (setq x
(cdr res
)))))
1568 ((and (eq (car-safe x
) 'frac
)
1570 (let ((res (math-integer-log (nth 2 x
) 10)))
1572 (setq x
(- (cdr res
))))))
1575 (if calc-infinite-mode
1576 '(neg (var inf var-inf
))
1577 (math-reject-arg x
"*Logarithm of zero")))
1578 (calc-symbolic-mode (signal 'inexact-result nil
))
1580 (math-with-extra-prec 2
1581 (let ((xf (math-float x
)))
1582 (if (eq (nth 1 xf
) 0)
1583 (math-reject-arg x
"*Logarithm of zero"))
1584 (if (Math-integer-posp (nth 1 xf
))
1585 (if (eq (nth 1 xf
) 1) ; log10(1*10^n) = n
1586 (math-float (nth 2 xf
))
1587 (let ((xdigs (1- (math-numdigs (nth 1 xf
)))))
1589 (math-div-float (math-ln-raw-2
1590 (list 'float
(nth 1 xf
) (- xdigs
)))
1592 (math-float (+ (nth 2 xf
) xdigs
)))))
1593 (math-div (calcFunc-ln xf
) (math-ln-10))))))
1594 ((eq (car-safe x
) 'sdev
)
1595 (math-make-sdev (calcFunc-log10 (nth 1 x
))
1597 (math-mul (nth 1 x
) (math-ln-10)))))
1598 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1599 (not (math-intv-constp x
))))
1600 (math-make-intv (nth 1 x
)
1601 (calcFunc-log10 (nth 2 x
))
1602 (calcFunc-log10 (nth 3 x
))))
1604 (if (equal x
'(var nan var-nan
))
1606 '(var inf var-inf
)))
1607 (t (calc-record-why 'numberp x
)
1608 (list 'calcFunc-log10 x
))))
1610 (defun calcFunc-log (x &optional b
) ; [N N N] [Public]
1611 (cond ((or (null b
) (equal b
'(var e var-e
)))
1612 (math-normalize (list 'calcFunc-ln x
)))
1613 ((or (eq b
10) (equal b
'(float 1 1)))
1614 (math-normalize (list 'calcFunc-log10 x
)))
1616 (if calc-infinite-mode
1617 (math-div (calcFunc-ln x
) (calcFunc-ln b
))
1618 (math-reject-arg x
"*Logarithm of zero")))
1620 (if calc-infinite-mode
1621 (math-div (calcFunc-ln x
) (calcFunc-ln b
))
1622 (math-reject-arg b
"*Logarithm of zero")))
1623 ((math-equal-int b
1)
1624 (if calc-infinite-mode
1625 (math-div (calcFunc-ln x
) 0)
1626 (math-reject-arg b
"*Logarithm base one")))
1627 ((math-equal-int x
1)
1628 (if (math-floatp b
) '(float 0 0) 0))
1629 ((and (Math-ratp x
) (Math-ratp b
)
1630 (math-posp x
) (math-posp b
)
1631 (let* ((sign 1) (inv nil
)
1632 (xx (if (Math-lessp 1 x
)
1636 (bb (if (Math-lessp 1 b
)
1638 (setq sign
(- sign
))
1640 (res (if (Math-lessp xx bb
)
1641 (setq inv
(math-integer-log bb xx
))
1642 (math-integer-log xx bb
))))
1645 (math-div 1 (* sign
(cdr res
)))
1646 (* sign
(cdr res
)))))))
1648 (calc-symbolic-mode (signal 'inexact-result nil
))
1649 ((and (Math-numberp x
) (Math-numberp b
))
1650 (math-with-extra-prec 2
1651 (math-div (math-ln-raw (math-float x
))
1652 (math-log-base-raw b
))))
1653 ((and (eq (car-safe x
) 'sdev
)
1655 (math-make-sdev (calcFunc-log (nth 1 x
) b
)
1658 (math-log-base-raw b
)))))
1659 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1660 (not (math-intv-constp x
)))
1662 (math-make-intv (nth 1 x
)
1663 (calcFunc-log (nth 2 x
) b
)
1664 (calcFunc-log (nth 3 x
) b
)))
1665 ((or (eq (car-safe x
) 'intv
) (eq (car-safe b
) 'intv
))
1666 (math-div (calcFunc-ln x
) (calcFunc-ln b
)))
1667 ((or (math-infinitep x
)
1669 (math-div (calcFunc-ln x
) (calcFunc-ln b
)))
1670 (t (if (Math-numberp b
)
1671 (calc-record-why 'numberp x
)
1672 (calc-record-why 'numberp b
))
1673 (list 'calcFunc-log x b
))))
1675 (defun calcFunc-alog (x &optional b
)
1676 (cond ((or (null b
) (equal b
'(var e var-e
)))
1677 (math-normalize (list 'calcFunc-exp x
)))
1678 (t (math-pow b x
))))
1680 (defun calcFunc-ilog (x b
)
1681 (if (and (math-natnump x
) (not (eq x
0))
1682 (math-natnump b
) (not (eq b
0)))
1684 (math-reject-arg x
"*Logarithm base one")
1685 (if (Math-natnum-lessp x b
)
1687 (cdr (math-integer-log x b
))))
1688 (math-floor (calcFunc-log x b
))))
1690 (defun math-integer-log (x b
)
1691 (let ((pows (list b
))
1695 (while (not (Math-lessp x pow
))
1696 (setq pows
(cons pow pows
)
1697 pow
(math-sqr pow
)))
1698 (setq n
(lsh 1 (1- (length pows
)))
1701 (while (and (setq pows
(cdr pows
))
1704 next
(math-mul pow
(car pows
)))
1705 (or (Math-lessp x next
)
1708 (cons (equal pow x
) sum
)))
1711 (defvar math-log-base-cache nil
)
1712 (defun math-log-base-raw (b) ; [N N]
1713 (if (not (and (equal (car math-log-base-cache
) b
)
1714 (eq (nth 1 math-log-base-cache
) calc-internal-prec
)))
1715 (setq math-log-base-cache
(list b calc-internal-prec
1716 (math-ln-raw (math-float b
)))))
1717 (nth 2 math-log-base-cache
))
1719 (defun calcFunc-lnp1 (x) ; [N N] [Public]
1720 (cond ((Math-equal-int x -
1)
1721 (if calc-infinite-mode
1722 '(neg (var inf var-inf
))
1723 (math-reject-arg x
"*Logarithm of zero")))
1725 ((math-zerop x
) '(float 0 0))
1726 (calc-symbolic-mode (signal 'inexact-result nil
))
1728 (math-with-extra-prec 2
1729 (let ((x (math-float x
)))
1730 (if (and (eq (car x
) 'float
)
1731 (math-lessp-float x
'(float 5 -
1))
1732 (math-lessp-float '(float -
5 -
1) x
))
1733 (math-ln-plus-1-raw x
)
1734 (math-ln-raw (math-add-float x
'(float 1 0)))))))
1735 ((eq (car-safe x
) 'sdev
)
1736 (math-make-sdev (calcFunc-lnp1 (nth 1 x
))
1737 (math-div (nth 2 x
) (math-add (nth 1 x
) 1))))
1738 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1739 (not (math-intv-constp x
))))
1740 (math-make-intv (nth 1 x
)
1741 (calcFunc-lnp1 (nth 2 x
))
1742 (calcFunc-lnp1 (nth 3 x
))))
1744 (if (equal x
'(var nan var-nan
))
1746 '(var inf var-inf
)))
1747 (t (calc-record-why 'numberp x
)
1748 (list 'calcFunc-lnp1 x
))))
1750 (defun math-ln-raw (x) ; [N N] --- must be float format!
1751 (cond ((eq (car-safe x
) 'cplx
)
1753 (math-mul-float (math-ln-raw
1754 (math-add-float (math-sqr-float (nth 1 x
))
1755 (math-sqr-float (nth 2 x
))))
1757 (math-arctan2-raw (nth 2 x
) (nth 1 x
))))
1758 ((eq (car x
) 'polar
)
1759 (math-polar (list 'cplx
1760 (math-ln-raw (nth 1 x
))
1761 (math-to-radians (nth 2 x
)))))
1762 ((Math-equal-int x
1)
1764 (calc-symbolic-mode (signal 'inexact-result nil
))
1765 ((math-posp (nth 1 x
)) ; positive and real
1767 ((math-use-emacs-fn 'log x
))
1769 (let ((xdigs (1- (math-numdigs (nth 1 x
)))))
1770 (math-add-float (math-ln-raw-2 (list 'float
(nth 1 x
) (- xdigs
)))
1771 (math-mul-float (math-float (+ (nth 2 x
) xdigs
))
1774 (math-reject-arg x
"*Logarithm of zero"))
1775 ((eq calc-complex-mode
'polar
) ; negative and real
1777 (list 'cplx
; negative and real
1778 (math-ln-raw (math-neg-float x
))
1780 (t (list 'cplx
; negative and real
1781 (math-ln-raw (math-neg-float x
))
1784 (defun math-ln-raw-2 (x) ; [F F]
1785 (cond ((math-lessp-float '(float 14 -
1) x
)
1786 (math-add-float (math-ln-raw-2 (math-mul-float x
'(float 5 -
1)))
1788 (t ; now .7 < x <= 1.4
1789 (math-ln-raw-3 (math-div-float (math-sub-float x
'(float 1 0))
1790 (math-add-float x
'(float 1 0)))))))
1792 (defun math-ln-raw-3 (x) ; [F F]
1793 (math-mul-float (math-ln-raw-series x
3 x
(math-sqr-float x
))
1796 ;;; Compute ln((1+x)/(1-x))
1797 (defun math-ln-raw-series (sum n x xsqr
)
1798 (math-working "log" sum
)
1799 (let* ((nextx (math-mul-float x xsqr
))
1800 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1801 (if (math-nearly-equal-float sum nextsum
)
1803 (math-ln-raw-series nextsum
(+ n
2) nextx xsqr
))))
1805 (defun math-ln-plus-1-raw (x)
1806 (math-lnp1-series x
2 x
(math-neg x
)))
1808 (defun math-lnp1-series (sum n xpow x
)
1809 (math-working "lnp1" sum
)
1810 (let* ((nextx (math-mul-float xpow x
))
1811 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1812 (if (math-nearly-equal-float sum nextsum
)
1814 (math-lnp1-series nextsum
(1+ n
) nextx x
))))
1816 (defconst math-approx-ln-10
1817 (math-read-number-simple "2.302585092994045684018")
1818 "An approximation for ln(10).")
1820 (math-defcache math-ln-10 math-approx-ln-10
1821 (math-ln-raw-2 '(float 1 1)))
1823 (defconst math-approx-ln-2
1824 (math-read-number-simple "0.693147180559945309417")
1825 "An approximation for ln(2).")
1827 (math-defcache math-ln-2 math-approx-ln-2
1828 (math-ln-raw-3 (math-float '(frac 1 3))))
1832 ;;; Hyperbolic functions.
1834 (defun calcFunc-sinh (x) ; [N N] [Public]
1836 (math-expand-formulas
1838 (list '/ (list '-
(list 'calcFunc-exp x
)
1839 (list 'calcFunc-exp
(list 'neg x
))) 2)))
1841 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1842 (math-with-extra-prec 2
1843 (let ((expx (math-exp-raw (math-float x
))))
1844 (math-mul (math-add expx
(math-div -
1 expx
)) '(float 5 -
1)))))
1845 ((eq (car-safe x
) 'sdev
)
1846 (math-make-sdev (calcFunc-sinh (nth 1 x
))
1847 (math-mul (nth 2 x
) (calcFunc-cosh (nth 1 x
)))))
1849 (math-sort-intv (nth 1 x
)
1850 (calcFunc-sinh (nth 2 x
))
1851 (calcFunc-sinh (nth 3 x
))))
1852 ((or (equal x
'(var inf var-inf
))
1853 (equal x
'(neg (var inf var-inf
)))
1854 (equal x
'(var nan var-nan
)))
1856 (t (calc-record-why 'numberp x
)
1857 (list 'calcFunc-sinh x
))))
1858 (put 'calcFunc-sinh
'math-expandable t
)
1860 (defun calcFunc-cosh (x) ; [N N] [Public]
1862 (math-expand-formulas
1864 (list '/ (list '+ (list 'calcFunc-exp x
)
1865 (list 'calcFunc-exp
(list 'neg x
))) 2)))
1867 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1868 (math-with-extra-prec 2
1869 (let ((expx (math-exp-raw (math-float x
))))
1870 (math-mul (math-add expx
(math-div 1 expx
)) '(float 5 -
1)))))
1871 ((eq (car-safe x
) 'sdev
)
1872 (math-make-sdev (calcFunc-cosh (nth 1 x
))
1874 (calcFunc-sinh (nth 1 x
)))))
1875 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1876 (setq x
(math-abs x
))
1877 (math-sort-intv (nth 1 x
)
1878 (calcFunc-cosh (nth 2 x
))
1879 (calcFunc-cosh (nth 3 x
))))
1880 ((or (equal x
'(var inf var-inf
))
1881 (equal x
'(neg (var inf var-inf
)))
1882 (equal x
'(var nan var-nan
)))
1884 (t (calc-record-why 'numberp x
)
1885 (list 'calcFunc-cosh x
))))
1886 (put 'calcFunc-cosh
'math-expandable t
)
1888 (defun calcFunc-tanh (x) ; [N N] [Public]
1890 (math-expand-formulas
1892 (let ((expx (list 'calcFunc-exp x
))
1893 (expmx (list 'calcFunc-exp
(list 'neg x
))))
1895 (list '/ (list '- expx expmx
) (list '+ expx expmx
))))))
1897 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1898 (math-with-extra-prec 2
1899 (let* ((expx (calcFunc-exp (math-float x
)))
1900 (expmx (math-div 1 expx
)))
1901 (math-div (math-sub expx expmx
)
1902 (math-add expx expmx
)))))
1903 ((eq (car-safe x
) 'sdev
)
1904 (math-make-sdev (calcFunc-tanh (nth 1 x
))
1906 (math-sqr (calcFunc-cosh (nth 1 x
))))))
1908 (math-sort-intv (nth 1 x
)
1909 (calcFunc-tanh (nth 2 x
))
1910 (calcFunc-tanh (nth 3 x
))))
1911 ((equal x
'(var inf var-inf
))
1913 ((equal x
'(neg (var inf var-inf
)))
1915 ((equal x
'(var nan var-nan
))
1917 (t (calc-record-why 'numberp x
)
1918 (list 'calcFunc-tanh x
))))
1919 (put 'calcFunc-tanh
'math-expandable t
)
1921 (defun calcFunc-sech (x) ; [N N] [Public]
1923 (math-expand-formulas
1925 (list '/ 2 (list '+ (list 'calcFunc-exp x
)
1926 (list 'calcFunc-exp
(list 'neg x
))))))
1928 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1929 (math-with-extra-prec 2
1930 (let ((expx (math-exp-raw (math-float x
))))
1931 (math-div '(float 2 0) (math-add expx
(math-div 1 expx
))))))
1932 ((eq (car-safe x
) 'sdev
)
1933 (math-make-sdev (calcFunc-sech (nth 1 x
))
1935 (math-mul (calcFunc-sech (nth 1 x
))
1936 (calcFunc-tanh (nth 1 x
))))))
1937 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1938 (setq x
(math-abs x
))
1939 (math-sort-intv (nth 1 x
)
1940 (calcFunc-sech (nth 2 x
))
1941 (calcFunc-sech (nth 3 x
))))
1942 ((or (equal x
'(var inf var-inf
))
1943 (equal x
'(neg (var inf var-inf
))))
1945 ((equal x
'(var nan var-nan
))
1947 (t (calc-record-why 'numberp x
)
1948 (list 'calcFunc-sech x
))))
1949 (put 'calcFunc-sech
'math-expandable t
)
1951 (defun calcFunc-csch (x) ; [N N] [Public]
1952 (cond ((eq x
0) (math-div 1 0))
1953 (math-expand-formulas
1955 (list '/ 2 (list '-
(list 'calcFunc-exp x
)
1956 (list 'calcFunc-exp
(list 'neg x
))))))
1958 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1959 (math-with-extra-prec 2
1960 (let ((expx (math-exp-raw (math-float x
))))
1961 (math-div '(float 2 0) (math-add expx
(math-div -
1 expx
))))))
1962 ((eq (car-safe x
) 'sdev
)
1963 (math-make-sdev (calcFunc-csch (nth 1 x
))
1965 (math-mul (calcFunc-csch (nth 1 x
))
1966 (calcFunc-coth (nth 1 x
))))))
1968 (if (and (Math-negp (nth 2 x
))
1969 (Math-posp (nth 3 x
)))
1970 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
1971 (math-sort-intv (nth 1 x
)
1972 (calcFunc-csch (nth 2 x
))
1973 (calcFunc-csch (nth 3 x
)))))
1974 ((or (equal x
'(var inf var-inf
))
1975 (equal x
'(neg (var inf var-inf
))))
1977 ((equal x
'(var nan var-nan
))
1979 (t (calc-record-why 'numberp x
)
1980 (list 'calcFunc-csch x
))))
1981 (put 'calcFunc-csch
'math-expandable t
)
1983 (defun calcFunc-coth (x) ; [N N] [Public]
1984 (cond ((eq x
0) (math-div 1 0))
1985 (math-expand-formulas
1987 (let ((expx (list 'calcFunc-exp x
))
1988 (expmx (list 'calcFunc-exp
(list 'neg x
))))
1990 (list '/ (list '+ expx expmx
) (list '- expx expmx
))))))
1992 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1993 (math-with-extra-prec 2
1994 (let* ((expx (calcFunc-exp (math-float x
)))
1995 (expmx (math-div 1 expx
)))
1996 (math-div (math-add expx expmx
)
1997 (math-sub expx expmx
)))))
1998 ((eq (car-safe x
) 'sdev
)
1999 (math-make-sdev (calcFunc-coth (nth 1 x
))
2001 (math-sqr (calcFunc-sinh (nth 1 x
))))))
2003 (if (and (Math-negp (nth 2 x
))
2004 (Math-posp (nth 3 x
)))
2005 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
2006 (math-sort-intv (nth 1 x
)
2007 (calcFunc-coth (nth 2 x
))
2008 (calcFunc-coth (nth 3 x
)))))
2009 ((equal x
'(var inf var-inf
))
2011 ((equal x
'(neg (var inf var-inf
)))
2013 ((equal x
'(var nan var-nan
))
2015 (t (calc-record-why 'numberp x
)
2016 (list 'calcFunc-coth x
))))
2017 (put 'calcFunc-coth
'math-expandable t
)
2019 (defun calcFunc-arcsinh (x) ; [N N] [Public]
2021 (math-expand-formulas
2023 (list 'calcFunc-ln
(list '+ x
(list 'calcFunc-sqrt
2024 (list '+ (list '^ x
2) 1))))))
2026 (if calc-symbolic-mode
(signal 'inexact-result nil
))
2027 (math-with-extra-prec 2
2028 (math-ln-raw (math-add x
(math-sqrt-raw (math-add (math-sqr x
)
2030 ((eq (car-safe x
) 'sdev
)
2031 (math-make-sdev (calcFunc-arcsinh (nth 1 x
))
2034 (math-add (math-sqr (nth 1 x
)) 1)))))
2036 (math-sort-intv (nth 1 x
)
2037 (calcFunc-arcsinh (nth 2 x
))
2038 (calcFunc-arcsinh (nth 3 x
))))
2039 ((or (equal x
'(var inf var-inf
))
2040 (equal x
'(neg (var inf var-inf
)))
2041 (equal x
'(var nan var-nan
)))
2043 (t (calc-record-why 'numberp x
)
2044 (list 'calcFunc-arcsinh x
))))
2045 (put 'calcFunc-arcsinh
'math-expandable t
)
2047 (defun calcFunc-arccosh (x) ; [N N] [Public]
2049 ((and (eq x -
1) calc-symbolic-mode
)
2051 ((and (eq x
0) calc-symbolic-mode
)
2052 (math-div (math-mul '(var pi var-pi
) '(var i var-i
)) 2))
2053 (math-expand-formulas
2055 (list 'calcFunc-ln
(list '+ x
(list 'calcFunc-sqrt
2056 (list '-
(list '^ x
2) 1))))))
2058 (if calc-symbolic-mode
(signal 'inexact-result nil
))
2059 (if (Math-equal-int x -
1)
2060 (math-imaginary (math-pi))
2061 (math-with-extra-prec 2
2062 (if (or t
; need to do this even in the real case!
2063 (memq (car-safe x
) '(cplx polar
)))
2064 (let ((xp1 (math-add 1 x
))) ; this gets the branch cuts right
2066 (math-add x
(math-mul xp1
2073 (math-add x
(math-sqrt-raw (math-add (math-sqr x
)
2074 '(float -
1 0)))))))))
2075 ((eq (car-safe x
) 'sdev
)
2076 (math-make-sdev (calcFunc-arccosh (nth 1 x
))
2079 (math-add (math-sqr (nth 1 x
)) -
1)))))
2081 (math-sort-intv (nth 1 x
)
2082 (calcFunc-arccosh (nth 2 x
))
2083 (calcFunc-arccosh (nth 3 x
))))
2084 ((or (equal x
'(var inf var-inf
))
2085 (equal x
'(neg (var inf var-inf
)))
2086 (equal x
'(var nan var-nan
)))
2088 (t (calc-record-why 'numberp x
)
2089 (list 'calcFunc-arccosh x
))))
2090 (put 'calcFunc-arccosh
'math-expandable t
)
2092 (defun calcFunc-arctanh (x) ; [N N] [Public]
2094 ((and (Math-equal-int x
1) calc-infinite-mode
)
2096 ((and (Math-equal-int x -
1) calc-infinite-mode
)
2097 '(neg (var inf var-inf
)))
2098 (math-expand-formulas
2100 (list 'calcFunc-ln
(list '+ 1 x
))
2101 (list 'calcFunc-ln
(list '-
1 x
))) 2))
2103 (if calc-symbolic-mode
(signal 'inexact-result nil
))
2104 (math-with-extra-prec 2
2105 (if (or (memq (car-safe x
) '(cplx polar
))
2107 (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x
))
2108 (math-ln-raw (math-sub '(float 1 0) x
)))
2110 (if (and (math-equal-int x
1) calc-infinite-mode
)
2112 (if (and (math-equal-int x -
1) calc-infinite-mode
)
2113 '(neg (var inf var-inf
))
2114 (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x
)
2117 ((eq (car-safe x
) 'sdev
)
2118 (math-make-sdev (calcFunc-arctanh (nth 1 x
))
2120 (math-sub 1 (math-sqr (nth 1 x
))))))
2122 (math-sort-intv (nth 1 x
)
2123 (calcFunc-arctanh (nth 2 x
))
2124 (calcFunc-arctanh (nth 3 x
))))
2125 ((equal x
'(var nan var-nan
))
2127 (t (calc-record-why 'numberp x
)
2128 (list 'calcFunc-arctanh x
))))
2129 (put 'calcFunc-arctanh
'math-expandable t
)
2132 ;;; Convert A from HMS or degrees to radians.
2133 (defun calcFunc-rad (a) ; [R R] [Public]
2134 (cond ((or (Math-numberp a
)
2136 (math-with-extra-prec 2
2137 (math-mul a
(math-pi-over-180))))
2139 (math-from-hms a
'rad
))
2141 (math-make-sdev (calcFunc-rad (nth 1 a
))
2142 (calcFunc-rad (nth 2 a
))))
2143 (math-expand-formulas
2144 (math-div (math-mul a
'(var pi var-pi
)) 180))
2145 ((math-infinitep a
) a
)
2146 (t (list 'calcFunc-rad a
))))
2147 (put 'calcFunc-rad
'math-expandable t
)
2149 ;;; Convert A from HMS or radians to degrees.
2150 (defun calcFunc-deg (a) ; [R R] [Public]
2151 (cond ((or (Math-numberp a
)
2153 (math-with-extra-prec 2
2154 (math-div a
(math-pi-over-180))))
2156 (math-from-hms a
'deg
))
2158 (math-make-sdev (calcFunc-deg (nth 1 a
))
2159 (calcFunc-deg (nth 2 a
))))
2160 (math-expand-formulas
2161 (math-div (math-mul 180 a
) '(var pi var-pi
)))
2162 ((math-infinitep a
) a
)
2163 (t (list 'calcFunc-deg a
))))
2164 (put 'calcFunc-deg
'math-expandable t
)
2166 (provide 'calc-math
)
2168 ;; arch-tag: c7367e8e-d0b8-4f70-8577-2fb3f31dbb4c
2169 ;;; calc-math.el ends here