1 ;;; calc-math.el --- mathematical functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
4 ;; 2005, 2006, 2007, 2008 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
)
108 (string-to-number (math-format-number (math-float x
))))))
111 (defun calc-sqrt (arg)
114 (if (calc-is-inverse)
115 (calc-unary-op "^2" 'calcFunc-sqr arg
)
116 (calc-unary-op "sqrt" 'calcFunc-sqrt arg
))))
118 (defun calc-isqrt (arg)
121 (if (calc-is-inverse)
122 (calc-unary-op "^2" 'calcFunc-sqr arg
)
123 (calc-unary-op "isqt" 'calcFunc-isqrt arg
))))
126 (defun calc-hypot (arg)
129 (calc-binary-op "hypt" 'calcFunc-hypot arg
)))
136 (defun calc-log10 (arg)
138 (calc-hyperbolic-func)
141 (defun calc-log (arg)
144 (if (calc-is-inverse)
145 (calc-binary-op "alog" 'calcFunc-alog arg
)
146 (calc-binary-op "log" 'calcFunc-log arg
))))
148 (defun calc-ilog (arg)
151 (if (calc-is-inverse)
152 (calc-binary-op "alog" 'calcFunc-alog arg
)
153 (calc-binary-op "ilog" 'calcFunc-ilog arg
))))
155 (defun calc-lnp1 (arg)
160 (defun calc-exp (arg)
163 (if (calc-is-hyperbolic)
164 (if (calc-is-inverse)
165 (calc-unary-op "lg10" 'calcFunc-log10 arg
)
166 (calc-unary-op "10^" 'calcFunc-exp10 arg
))
167 (if (calc-is-inverse)
168 (calc-unary-op "ln" 'calcFunc-ln arg
)
169 (calc-unary-op "exp" 'calcFunc-exp arg
)))))
171 (defun calc-expm1 (arg)
174 (if (calc-is-inverse)
175 (calc-unary-op "ln+1" 'calcFunc-lnp1 arg
)
176 (calc-unary-op "ex-1" 'calcFunc-expm1 arg
))))
181 (if (calc-is-inverse)
182 (if (calc-is-hyperbolic)
183 (if calc-symbolic-mode
184 (calc-pop-push-record 0 "phi" '(var phi var-phi
))
185 (calc-pop-push-record 0 "phi" (math-phi)))
186 (if calc-symbolic-mode
187 (calc-pop-push-record 0 "gmma" '(var gamma var-gamma
))
188 (calc-pop-push-record 0 "gmma" (math-gamma-const))))
189 (if (calc-is-hyperbolic)
190 (if calc-symbolic-mode
191 (calc-pop-push-record 0 "e" '(var e var-e
))
192 (calc-pop-push-record 0 "e" (math-e)))
193 (if calc-symbolic-mode
194 (calc-pop-push-record 0 "pi" '(var pi var-pi
))
195 (calc-pop-push-record 0 "pi" (math-pi)))))))
197 (defun calc-sin (arg)
200 (if (calc-is-hyperbolic)
201 (if (calc-is-inverse)
202 (calc-unary-op "asnh" 'calcFunc-arcsinh arg
)
203 (calc-unary-op "sinh" 'calcFunc-sinh arg
))
204 (if (calc-is-inverse)
205 (calc-unary-op "asin" 'calcFunc-arcsin arg
)
206 (calc-unary-op "sin" 'calcFunc-sin arg
)))))
208 (defun calc-arcsin (arg)
213 (defun calc-sinh (arg)
215 (calc-hyperbolic-func)
218 (defun calc-arcsinh (arg)
221 (calc-hyperbolic-func)
224 (defun calc-sec (arg)
227 (if (calc-is-hyperbolic)
228 (calc-unary-op "sech" 'calcFunc-sech arg
)
229 (calc-unary-op "sec" 'calcFunc-sec arg
))))
231 (defun calc-sech (arg)
233 (calc-hyperbolic-func)
236 (defun calc-cos (arg)
239 (if (calc-is-hyperbolic)
240 (if (calc-is-inverse)
241 (calc-unary-op "acsh" 'calcFunc-arccosh arg
)
242 (calc-unary-op "cosh" 'calcFunc-cosh arg
))
243 (if (calc-is-inverse)
244 (calc-unary-op "acos" 'calcFunc-arccos arg
)
245 (calc-unary-op "cos" 'calcFunc-cos arg
)))))
247 (defun calc-arccos (arg)
252 (defun calc-cosh (arg)
254 (calc-hyperbolic-func)
257 (defun calc-arccosh (arg)
260 (calc-hyperbolic-func)
263 (defun calc-csc (arg)
266 (if (calc-is-hyperbolic)
267 (calc-unary-op "csch" 'calcFunc-csch arg
)
268 (calc-unary-op "csc" 'calcFunc-csc arg
))))
270 (defun calc-csch (arg)
272 (calc-hyperbolic-func)
275 (defun calc-sincos ()
278 (if (calc-is-inverse)
279 (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos
(calc-top-n 1)))
280 (calc-enter-result 1 "sncs" (list 'calcFunc-sincos
(calc-top-n 1))))))
282 (defun calc-tan (arg)
285 (if (calc-is-hyperbolic)
286 (if (calc-is-inverse)
287 (calc-unary-op "atnh" 'calcFunc-arctanh arg
)
288 (calc-unary-op "tanh" 'calcFunc-tanh arg
))
289 (if (calc-is-inverse)
290 (calc-unary-op "atan" 'calcFunc-arctan arg
)
291 (calc-unary-op "tan" 'calcFunc-tan arg
)))))
293 (defun calc-arctan (arg)
298 (defun calc-tanh (arg)
300 (calc-hyperbolic-func)
303 (defun calc-arctanh (arg)
306 (calc-hyperbolic-func)
309 (defun calc-cot (arg)
312 (if (calc-is-hyperbolic)
313 (calc-unary-op "coth" 'calcFunc-coth arg
)
314 (calc-unary-op "cot" 'calcFunc-cot arg
))))
316 (defun calc-coth (arg)
318 (calc-hyperbolic-func)
321 (defun calc-arctan2 ()
324 (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2
(calc-top-list-n 2)))))
326 (defun calc-conj (arg)
329 (calc-unary-op "conj" 'calcFunc-conj arg
)))
331 (defun calc-imaginary ()
334 (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))))
336 (defun calc-to-degrees (arg)
339 (calc-unary-op ">deg" 'calcFunc-deg arg
)))
341 (defun calc-to-radians (arg)
344 (calc-unary-op ">rad" 'calcFunc-rad arg
)))
347 (defun calc-degrees-mode (arg)
351 (calc-change-mode 'calc-angle-mode
'deg
)
352 (message "Angles measured in degrees")))
353 ((= arg
2) (calc-radians-mode))
354 ((= arg
3) (calc-hms-mode))
355 (t (error "Prefix argument out of range"))))
357 (defun calc-radians-mode ()
360 (calc-change-mode 'calc-angle-mode
'rad
)
361 (message "Angles measured in radians")))
364 ;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public]
365 ;;; This method takes advantage of the fact that Newton's method starting
366 ;;; with an overestimate always works, even using truncating integer division!
367 (defun math-isqrt (a)
368 (cond ((Math-zerop a
) a
)
369 ((not (math-natnump a
))
370 (math-reject-arg a
'natnump
))
372 (math-isqrt-small a
))
374 (math-normalize (cons 'bigpos
(cdr (math-isqrt-bignum (cdr a
))))))))
376 (defun calcFunc-isqrt (a)
378 (math-isqrt (math-floor a
))
379 (math-floor (math-sqrt a
))))
382 ;;; This returns (flag . result) where the flag is t if A is a perfect square.
383 (defun math-isqrt-bignum (a) ; [P.l L]
384 (let ((len (length a
)))
386 (let* ((top (nthcdr (- len
2) a
)))
387 (math-isqrt-bignum-iter
389 (math-scale-bignum-digit-size
391 (1+ (math-isqrt-small
392 (+ (* (nth 1 top
) math-bignum-digit-size
) (car top
)))))
394 (let* ((top (nth (1- len
) a
)))
395 (math-isqrt-bignum-iter
397 (math-scale-bignum-digit-size
398 (list (1+ (math-isqrt-small top
)))
401 (defun math-isqrt-bignum-iter (a guess
) ; [l L l]
402 (math-working "isqrt" (cons 'bigpos guess
))
403 (let* ((q (math-div-bignum a guess
))
404 (s (math-add-bignum (car q
) guess
))
405 (g2 (math-div2-bignum s
))
406 (comp (math-compare-bignum g2 guess
)))
408 (math-isqrt-bignum-iter a g2
)
409 (cons (and (= comp
0)
410 (math-zerop-bignum (cdr q
))
414 (defun math-zerop-bignum (a)
417 (while (eq (car (setq a
(cdr a
))) 0))
420 (defun math-scale-bignum-digit-size (a n
) ; [L L S]
426 (defun math-isqrt-small (a) ; A > 0. [S S]
427 (let ((g (cond ((>= a
1000000) 10000)
432 (while (< (setq g2
(/ (+ g
(/ a g
)) 2)) g
)
439 ;;; Compute the square root of a number.
440 ;;; [T N] if possible, else [F N] if possible, else [C N]. [Public]
443 (and (Math-zerop a
) a
)
444 (and (math-known-nonposp a
)
445 (math-imaginary (math-sqrt (math-neg a
))))
447 (let ((sqrt (math-isqrt-small a
)))
448 (if (= (* sqrt sqrt
) a
)
450 (if calc-symbolic-mode
451 (list 'calcFunc-sqrt a
)
452 (math-sqrt-float (math-float a
) (math-float sqrt
))))))
453 (and (eq (car-safe a
) 'bigpos
)
454 (let* ((res (math-isqrt-bignum (cdr a
)))
455 (sqrt (math-normalize (cons 'bigpos
(cdr res
)))))
458 (if calc-symbolic-mode
459 (list 'calcFunc-sqrt a
)
460 (math-sqrt-float (math-float a
) (math-float sqrt
))))))
461 (and (eq (car-safe a
) 'frac
)
462 (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a
)))))
463 (num-sqrt (math-normalize (cons 'bigpos
(cdr num-res
))))
464 (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a
)))))
465 (den-sqrt (math-normalize (cons 'bigpos
(cdr den-res
)))))
466 (if (and (car num-res
) (car den-res
))
467 (list 'frac num-sqrt den-sqrt
)
468 (if calc-symbolic-mode
469 (if (or (car num-res
) (car den-res
))
470 (math-div (if (car num-res
)
471 num-sqrt
(list 'calcFunc-sqrt
(nth 1 a
)))
473 den-sqrt
(list 'calcFunc-sqrt
(nth 2 a
))))
474 (list 'calcFunc-sqrt a
))
475 (math-sqrt-float (math-float a
)
476 (math-div (math-float num-sqrt
) den-sqrt
))))))
477 (and (eq (car-safe a
) 'float
)
478 (if calc-symbolic-mode
479 (if (= (%
(nth 2 a
) 2) 0)
480 (let ((res (math-isqrt-bignum
481 (cdr (Math-bignum-test (nth 1 a
))))))
483 (math-make-float (math-normalize
484 (cons 'bigpos
(cdr res
)))
486 (signal 'inexact-result nil
)))
487 (signal 'inexact-result nil
))
488 (math-sqrt-float a
)))
489 (and (eq (car-safe a
) 'cplx
)
490 (math-with-extra-prec 2
491 (let* ((d (math-abs a
))
492 (imag (math-sqrt (math-mul (math-sub d
(nth 1 a
))
495 (math-sqrt (math-mul (math-add d
(nth 1 a
)) '(float 5 -
1)))
496 (if (math-negp (nth 2 a
)) (math-neg imag
) imag
)))))
497 (and (eq (car-safe a
) 'polar
)
499 (math-sqrt (nth 1 a
))
500 (math-mul (nth 2 a
) '(float 5 -
1))))
501 (and (eq (car-safe a
) 'sdev
)
502 (let ((sqrt (math-sqrt (nth 1 a
))))
504 (math-div (nth 2 a
) (math-mul sqrt
2)))))
505 (and (eq (car-safe a
) 'intv
)
506 (not (math-negp (nth 2 a
)))
507 (math-make-intv (nth 1 a
) (math-sqrt (nth 2 a
)) (math-sqrt (nth 3 a
))))
508 (and (eq (car-safe a
) '*)
509 (or (math-known-nonnegp (nth 1 a
))
510 (math-known-nonnegp (nth 2 a
)))
511 (math-mul (math-sqrt (nth 1 a
)) (math-sqrt (nth 2 a
))))
512 (and (eq (car-safe a
) '/)
513 (or (and (math-known-nonnegp (nth 2 a
))
514 (math-div (math-sqrt (nth 1 a
)) (math-sqrt (nth 2 a
))))
515 (and (math-known-nonnegp (nth 1 a
))
516 (not (math-equal-int (nth 1 a
) 1))
517 (math-mul (math-sqrt (nth 1 a
))
518 (math-sqrt (math-div 1 (nth 2 a
)))))))
519 (and (eq (car-safe a
) '^
)
520 (math-known-evenp (nth 2 a
))
521 (math-known-realp (nth 1 a
))
522 (math-abs (math-pow (nth 1 a
) (math-div (nth 2 a
) 2))))
523 (let ((inf (math-infinitep a
)))
525 (math-mul (math-sqrt (math-infinite-dir a inf
)) inf
)))
527 (calc-record-why 'numberp a
)
528 (list 'calcFunc-sqrt a
))))
529 (defalias 'calcFunc-sqrt
'math-sqrt
)
531 (defun math-infinite-dir (a &optional inf
)
532 (or inf
(setq inf
(math-infinitep a
)))
533 (math-normalize (math-expr-subst a inf
1)))
535 (defun math-sqrt-float (a &optional guess
) ; [F F F]
536 (if calc-symbolic-mode
537 (signal 'inexact-result nil
)
538 (math-with-extra-prec 1 (math-sqrt-raw a guess
))))
540 (defun math-sqrt-raw (a &optional guess
) ; [F F F]
541 (if (not (Math-posp a
))
544 ((math-use-emacs-fn 'sqrt a
))
547 (let ((ldiff (- (math-numdigs (nth 1 a
)) 6)))
548 (or (= (%
(+ (nth 2 a
) ldiff
) 2) 0) (setq ldiff
(1+ ldiff
)))
549 (setq guess
(math-make-float (math-isqrt-small
550 (math-scale-int (nth 1 a
) (- ldiff
)))
551 (/ (+ (nth 2 a
) ldiff
) 2)))))
552 (math-sqrt-float-iter a guess
)))))
554 (defun math-sqrt-float-iter (a guess
) ; [F F F]
555 (math-working "sqrt" guess
)
556 (let ((g2 (math-mul-float (math-add-float guess
(math-div-float a guess
))
558 (if (math-nearly-equal-float g2 guess
)
560 (math-sqrt-float-iter a g2
))))
562 ;;; True if A and B differ only in the last digit of precision. [P F F]
563 (defun math-nearly-equal-float (a b
)
564 (let ((ediff (- (nth 2 a
) (nth 2 b
))))
565 (cond ((= ediff
0) ;; Expanded out for speed
566 (setq ediff
(math-add (Math-integer-neg (nth 1 a
)) (nth 1 b
)))
568 (and (not (consp ediff
))
571 (= (math-numdigs (nth 1 a
)) calc-internal-prec
))))
573 (setq ediff
(math-add (Math-integer-neg (nth 1 b
))
574 (math-scale-int (nth 1 a
) 1)))
575 (and (not (consp ediff
))
578 (= (math-numdigs (nth 1 b
)) calc-internal-prec
)))
580 (setq ediff
(math-add (Math-integer-neg (nth 1 a
))
581 (math-scale-int (nth 1 b
) 1)))
582 (and (not (consp ediff
))
585 (= (math-numdigs (nth 1 a
)) calc-internal-prec
))))))
587 (defun math-nearly-equal (a b
) ; [P N N] [Public]
588 (setq a
(math-float a
))
589 (setq b
(math-float b
))
590 (if (eq (car a
) 'polar
) (setq a
(math-complex a
)))
591 (if (eq (car b
) 'polar
) (setq b
(math-complex b
)))
592 (if (eq (car a
) 'cplx
)
593 (if (eq (car b
) 'cplx
)
594 (and (or (math-nearly-equal-float (nth 1 a
) (nth 1 b
))
595 (and (math-nearly-zerop-float (nth 1 a
) (nth 2 a
))
596 (math-nearly-zerop-float (nth 1 b
) (nth 2 b
))))
597 (or (math-nearly-equal-float (nth 2 a
) (nth 2 b
))
598 (and (math-nearly-zerop-float (nth 2 a
) (nth 1 a
))
599 (math-nearly-zerop-float (nth 2 b
) (nth 1 b
)))))
600 (and (math-nearly-equal-float (nth 1 a
) b
)
601 (math-nearly-zerop-float (nth 2 a
) b
)))
602 (if (eq (car b
) 'cplx
)
603 (and (math-nearly-equal-float a
(nth 1 b
))
604 (math-nearly-zerop-float a
(nth 2 b
)))
605 (math-nearly-equal-float a b
))))
607 ;;; True if A is nearly zero compared to B. [P F F]
608 (defun math-nearly-zerop-float (a b
)
610 (<= (+ (math-numdigs (nth 1 a
)) (nth 2 a
))
611 (1+ (- (+ (math-numdigs (nth 1 b
)) (nth 2 b
)) calc-internal-prec
)))))
613 (defun math-nearly-zerop (a b
) ; [P N R] [Public]
614 (setq a
(math-float a
))
615 (setq b
(math-float b
))
616 (if (eq (car a
) 'cplx
)
617 (and (math-nearly-zerop-float (nth 1 a
) b
)
618 (math-nearly-zerop-float (nth 2 a
) b
))
619 (if (eq (car a
) 'polar
)
620 (math-nearly-zerop-float (nth 1 a
) b
)
621 (math-nearly-zerop-float a b
))))
623 ;;; This implementation could be improved, accuracy-wise.
624 (defun math-hypot (a b
)
625 (cond ((Math-zerop a
) (math-abs b
))
626 ((Math-zerop b
) (math-abs a
))
627 ((not (Math-scalarp a
))
628 (if (math-infinitep a
)
629 (if (math-infinitep b
)
634 (calc-record-why 'scalarp a
)
635 (list 'calcFunc-hypot a b
)))
636 ((not (Math-scalarp b
))
637 (if (math-infinitep b
)
639 (calc-record-why 'scalarp b
)
640 (list 'calcFunc-hypot a b
)))
641 ((and (Math-numberp a
) (Math-numberp b
))
642 (math-with-extra-prec 1
643 (math-sqrt (math-add (calcFunc-abssqr a
) (calcFunc-abssqr b
)))))
644 ((eq (car-safe a
) 'hms
)
645 (if (eq (car-safe b
) 'hms
) ; this helps sdev's of hms forms
646 (math-to-hms (math-hypot (math-from-hms a
'deg
)
647 (math-from-hms b
'deg
)))
648 (math-to-hms (math-hypot (math-from-hms a
'deg
) b
))))
649 ((eq (car-safe b
) 'hms
)
650 (math-to-hms (math-hypot a
(math-from-hms b
'deg
))))
652 (defalias 'calcFunc-hypot
'math-hypot
)
654 (defun calcFunc-sqr (x)
659 (defun math-nth-root (a n
)
660 (cond ((= n
2) (math-sqrt a
))
664 (let ((root (math-nth-root-integer a n
)))
667 (and (not calc-symbolic-mode
)
668 (math-nth-root-float (math-float a
) n
669 (math-float (cdr root
)))))))
670 ((eq (car-safe a
) 'frac
)
671 (let* ((num-root (math-nth-root-integer (nth 1 a
) n
))
672 (den-root (math-nth-root-integer (nth 2 a
) n
)))
673 (if (and (car num-root
) (car den-root
))
674 (list 'frac
(cdr num-root
) (cdr den-root
))
675 (and (not calc-symbolic-mode
)
678 (math-div-float (math-float (cdr num-root
))
679 (math-float (cdr den-root
))))))))
680 ((eq (car-safe a
) 'float
)
681 (and (not calc-symbolic-mode
)
682 (math-nth-root-float a n
)))
683 ((eq (car-safe a
) 'polar
)
684 (let ((root (math-nth-root (nth 1 a
) n
)))
685 (and root
(list 'polar root
(math-div (nth 2 a
) n
)))))
688 ;; The variables math-nrf-n, math-nrf-nf and math-nrf-nfm1 are local
689 ;; to math-nth-root-float, but are used by math-nth-root-float-iter,
690 ;; which is called by math-nth-root-float.
693 (defvar math-nrf-nfm1
)
695 (defun math-nth-root-float (a math-nrf-n
&optional guess
)
696 (math-inexact-result)
697 (math-with-extra-prec 1
698 (let ((math-nrf-nf (math-float math-nrf-n
))
699 (math-nrf-nfm1 (math-float (1- math-nrf-n
))))
700 (math-nth-root-float-iter a
(or guess
702 1 (/ (+ (math-numdigs (nth 1 a
))
707 (defun math-nth-root-float-iter (a guess
)
708 (math-working "root" guess
)
709 (let ((g2 (math-div-float (math-add-float (math-mul math-nrf-nfm1 guess
)
711 a
(math-ipow guess
(1- math-nrf-n
))))
713 (if (math-nearly-equal-float g2 guess
)
715 (math-nth-root-float-iter a g2
))))
717 ;; The variable math-nri-n is local to math-nth-root-integer, but
718 ;; is used by math-nth-root-int-iter, which is called by
719 ;; math-nth-root-int.
722 (defun math-nth-root-integer (a math-nri-n
&optional guess
) ; [I I S]
723 (math-nth-root-int-iter a
(or guess
724 (math-scale-int 1 (/ (+ (math-numdigs a
)
728 (defun math-nth-root-int-iter (a guess
)
729 (math-working "root" guess
)
730 (let* ((q (math-idivmod a
(math-ipow guess
(1- math-nri-n
))))
731 (s (math-add (car q
) (math-mul (1- math-nri-n
) guess
)))
732 (g2 (math-idivmod s math-nri-n
)))
733 (if (Math-natnum-lessp (car g2
) guess
)
734 (math-nth-root-int-iter a
(car g2
))
735 (cons (and (equal (car g2
) guess
)
740 (defun calcFunc-nroot (x n
)
741 (calcFunc-pow x
(if (integerp n
)
748 ;;;; Transcendental functions.
750 ;;; All of these functions are defined on the complex plane.
751 ;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
753 ;;; Most functions increase calc-internal-prec by 2 digits, then round
754 ;;; down afterward. "-raw" functions use the current precision, require
755 ;;; their arguments to be in float (or complex float) format, and always
756 ;;; work in radians (where applicable).
758 (defun math-to-radians (a) ; [N N]
759 (cond ((eq (car-safe a
) 'hms
)
760 (math-from-hms a
'rad
))
761 ((memq calc-angle-mode
'(deg hms
))
762 (math-mul a
(math-pi-over-180)))
765 (defun math-from-radians (a) ; [N N]
766 (cond ((eq calc-angle-mode
'deg
)
768 (math-div a
(math-pi-over-180))
769 (list 'calcFunc-deg a
)))
770 ((eq calc-angle-mode
'hms
)
771 (math-to-hms a
'rad
))
774 (defun math-to-radians-2 (a) ; [N N]
775 (cond ((eq (car-safe a
) 'hms
)
776 (math-from-hms a
'rad
))
777 ((memq calc-angle-mode
'(deg hms
))
778 (if calc-symbolic-mode
779 (math-div (math-mul a
'(var pi var-pi
)) 180)
780 (math-mul a
(math-pi-over-180))))
783 (defun math-from-radians-2 (a) ; [N N]
784 (cond ((memq calc-angle-mode
'(deg hms
))
785 (if calc-symbolic-mode
786 (math-div (math-mul 180 a
) '(var pi var-pi
))
787 (math-div a
(math-pi-over-180))))
792 ;;; Sine, cosine, and tangent.
794 (defun calcFunc-sin (x) ; [N N] [Public]
795 (cond ((and (integerp x
)
796 (if (eq calc-angle-mode
'deg
)
799 (aref [0 1 0 -
1] (math-mod (/ x
90) 4)))
801 (math-with-extra-prec 2
802 (math-sin-raw (math-to-radians (math-float x
)))))
805 (math-with-extra-prec 2
806 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
807 (xs (math-to-radians (math-float (nth 2 x
))))
808 (sc (math-sin-cos-raw xx
)))
809 (math-make-sdev (car sc
) (math-mul xs
(cdr sc
)))))
810 (math-make-sdev (calcFunc-sin (nth 1 x
))
811 (math-mul (nth 2 x
) (calcFunc-cos (nth 1 x
))))))
812 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
813 (calcFunc-cos (math-sub x
(math-quarter-circle nil
))))
814 ((equal x
'(var nan var-nan
))
816 (t (calc-record-why 'scalarp x
)
817 (list 'calcFunc-sin x
))))
819 (defun calcFunc-cos (x) ; [N N] [Public]
820 (cond ((and (integerp x
)
821 (if (eq calc-angle-mode
'deg
)
824 (aref [1 0 -
1 0] (math-mod (/ x
90) 4)))
826 (math-with-extra-prec 2
827 (math-cos-raw (math-to-radians (math-float x
)))))
830 (math-with-extra-prec 2
831 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
832 (xs (math-to-radians (math-float (nth 2 x
))))
833 (sc (math-sin-cos-raw xx
)))
834 (math-make-sdev (cdr sc
) (math-mul xs
(car sc
)))))
835 (math-make-sdev (calcFunc-cos (nth 1 x
))
836 (math-mul (nth 2 x
) (calcFunc-sin (nth 1 x
))))))
837 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
838 (math-with-extra-prec 2
839 (let* ((xx (math-to-radians (math-float x
)))
840 (na (math-floor (math-div (nth 2 xx
) (math-pi))))
841 (nb (math-floor (math-div (nth 3 xx
) (math-pi))))
842 (span (math-sub nb na
)))
843 (if (memq span
'(0 1))
844 (let ((int (math-sort-intv (nth 1 x
)
845 (math-cos-raw (nth 2 xx
))
846 (math-cos-raw (nth 3 xx
)))))
849 (math-make-intv (logior (nth 1 x
) 2)
852 (math-make-intv (logior (nth 1 x
) 1)
856 (list 'intv
3 -
1 1)))))
857 ((equal x
'(var nan var-nan
))
859 (t (calc-record-why 'scalarp x
)
860 (list 'calcFunc-cos x
))))
862 (defun calcFunc-sincos (x) ; [V N] [Public]
864 (math-with-extra-prec 2
865 (let ((sc (math-sin-cos-raw (math-to-radians (math-float x
)))))
866 (list 'vec
(cdr sc
) (car sc
)))) ; the vector [cos, sin]
867 (list 'vec
(calcFunc-sin x
) (calcFunc-cos x
))))
869 (defun calcFunc-tan (x) ; [N N] [Public]
870 (cond ((and (integerp x
)
871 (if (eq calc-angle-mode
'deg
)
876 (math-with-extra-prec 2
877 (math-tan-raw (math-to-radians (math-float x
)))))
880 (math-with-extra-prec 2
881 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
882 (xs (math-to-radians (math-float (nth 2 x
))))
883 (sc (math-sin-cos-raw xx
)))
884 (if (and (math-zerop (cdr sc
)) (not calc-infinite-mode
))
886 (calc-record-why "*Division by zero")
887 (list 'calcFunc-tan x
))
888 (math-make-sdev (math-div-float (car sc
) (cdr sc
))
889 (math-div-float xs
(math-sqr (cdr sc
)))))))
890 (math-make-sdev (calcFunc-tan (nth 1 x
))
892 (math-sqr (calcFunc-cos (nth 1 x
)))))))
893 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
894 (or (math-with-extra-prec 2
895 (let* ((xx (math-to-radians (math-float x
)))
896 (na (math-floor (math-div (math-sub (nth 2 xx
)
899 (nb (math-floor (math-div (math-sub (nth 3 xx
)
903 (math-sort-intv (nth 1 x
)
904 (math-tan-raw (nth 2 xx
))
905 (math-tan-raw (nth 3 xx
))))))
906 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))))
907 ((equal x
'(var nan var-nan
))
909 (t (calc-record-why 'scalarp x
)
910 (list 'calcFunc-tan x
))))
912 (defun calcFunc-sec (x)
913 (cond ((and (integerp x
)
914 (eq calc-angle-mode
'deg
)
920 (eq calc-angle-mode
'rad
)
924 (math-with-extra-prec 2
925 (math-sec-raw (math-to-radians (math-float x
)))))
928 (math-with-extra-prec 2
929 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
930 (xs (math-to-radians (math-float (nth 2 x
))))
931 (sc (math-sin-cos-raw xx
)))
932 (if (and (math-zerop (cdr sc
))
933 (not calc-infinite-mode
))
935 (calc-record-why "*Division by zero")
936 (list 'calcFunc-sec x
))
937 (math-make-sdev (math-div-float '(float 1 0) (cdr sc
))
939 (math-mul xs
(car sc
))
940 (math-sqr (cdr sc
)))))))
941 (math-make-sdev (calcFunc-sec (nth 1 x
))
944 (calcFunc-sin (nth 1 x
)))
945 (math-sqr (calcFunc-cos (nth 1 x
)))))))
946 ((and (eq (car x
) 'intv
)
947 (math-intv-constp x
))
948 (math-with-extra-prec 2
949 (let* ((xx (math-to-radians (math-float x
)))
950 (na (math-floor (math-div (math-sub (nth 2 xx
)
953 (nb (math-floor (math-div (math-sub (nth 3 xx
)
956 (naa (math-floor (math-div (nth 2 xx
) (math-pi-over-2))))
957 (nbb (math-floor (math-div (nth 3 xx
) (math-pi-over-2))))
958 (span (math-sub nbb naa
)))
959 (if (not (equal na nb
))
960 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
961 (let ((int (math-sort-intv (nth 1 x
)
962 (math-sec-raw (nth 2 xx
))
963 (math-sec-raw (nth 3 xx
)))))
965 (if (math-evenp (math-div (math-add naa
1) 2))
966 (math-make-intv (logior (nth 1 int
) 2)
969 (math-make-intv (logior (nth 1 int
) 1)
973 ((equal x
'(var nan var-nan
))
975 (t (calc-record-why 'scalarp x
)
976 (list 'calcFunc-sec x
))))
978 (defun calcFunc-csc (x)
979 (cond ((and (integerp x
)
980 (eq calc-angle-mode
'deg
)
981 (= (%
(- x
90) 180) 0))
982 (if (= (%
(- x
90) 360) 0)
986 (math-with-extra-prec 2
987 (math-csc-raw (math-to-radians (math-float x
)))))
990 (math-with-extra-prec 2
991 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
992 (xs (math-to-radians (math-float (nth 2 x
))))
993 (sc (math-sin-cos-raw xx
)))
994 (if (and (math-zerop (car sc
))
995 (not calc-infinite-mode
))
997 (calc-record-why "*Division by zero")
998 (list 'calcFunc-csc x
))
999 (math-make-sdev (math-div-float '(float 1 0) (car sc
))
1001 (math-mul xs
(cdr sc
))
1002 (math-sqr (car sc
)))))))
1003 (math-make-sdev (calcFunc-csc (nth 1 x
))
1006 (calcFunc-cos (nth 1 x
)))
1007 (math-sqr (calcFunc-sin (nth 1 x
)))))))
1008 ((and (eq (car x
) 'intv
)
1009 (math-intv-constp x
))
1010 (math-with-extra-prec 2
1011 (let* ((xx (math-to-radians (math-float x
)))
1012 (na (math-floor (math-div (nth 2 xx
) (math-pi))))
1013 (nb (math-floor (math-div (nth 3 xx
) (math-pi))))
1014 (naa (math-floor (math-div (nth 2 xx
) (math-pi-over-2))))
1015 (nbb (math-floor (math-div (nth 3 xx
) (math-pi-over-2))))
1016 (span (math-sub nbb naa
)))
1017 (if (not (equal na nb
))
1018 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
1019 (let ((int (math-sort-intv (nth 1 x
)
1020 (math-csc-raw (nth 2 xx
))
1021 (math-csc-raw (nth 3 xx
)))))
1023 (if (math-evenp (math-div naa
2))
1024 (math-make-intv (logior (nth 1 int
) 2)
1027 (math-make-intv (logior (nth 1 int
) 1)
1031 ((equal x
'(var nan var-nan
))
1033 (t (calc-record-why 'scalarp x
)
1034 (list 'calcFunc-csc x
))))
1036 (defun calcFunc-cot (x) ; [N N] [Public]
1037 (cond ((and (integerp x
)
1038 (if (eq calc-angle-mode
'deg
)
1039 (= (%
(- x
90) 180) 0)
1043 (math-with-extra-prec 2
1044 (math-cot-raw (math-to-radians (math-float x
)))))
1047 (math-with-extra-prec 2
1048 (let* ((xx (math-to-radians (math-float (nth 1 x
))))
1049 (xs (math-to-radians (math-float (nth 2 x
))))
1050 (sc (math-sin-cos-raw xx
)))
1051 (if (and (math-zerop (car sc
)) (not calc-infinite-mode
))
1053 (calc-record-why "*Division by zero")
1054 (list 'calcFunc-cot x
))
1055 (math-make-sdev (math-div-float (cdr sc
) (car sc
))
1056 (math-div-float xs
(math-sqr (car sc
)))))))
1057 (math-make-sdev (calcFunc-cot (nth 1 x
))
1059 (math-sqr (calcFunc-sin (nth 1 x
)))))))
1060 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1061 (or (math-with-extra-prec 2
1062 (let* ((xx (math-to-radians (math-float x
)))
1063 (na (math-floor (math-div (nth 2 xx
) (math-pi))))
1064 (nb (math-floor (math-div (nth 3 xx
) (math-pi)))))
1066 (math-sort-intv (nth 1 x
)
1067 (math-cot-raw (nth 2 xx
))
1068 (math-cot-raw (nth 3 xx
))))))
1069 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))))
1070 ((equal x
'(var nan var-nan
))
1072 (t (calc-record-why 'scalarp x
)
1073 (list 'calcFunc-cot x
))))
1075 (defun math-sin-raw (x &optional orgx
) ; [N N]
1076 (cond ((eq (car x
) 'cplx
)
1077 (let* ((expx (math-exp-raw (nth 2 x
)))
1078 (expmx (math-div-float '(float 1 0) expx
))
1079 (sc (math-sin-cos-raw (nth 1 x
))))
1081 (math-mul-float (car sc
)
1082 (math-mul-float (math-add-float expx expmx
)
1084 (math-mul-float (cdr sc
)
1085 (math-mul-float (math-sub-float expx expmx
)
1087 ((eq (car x
) 'polar
)
1088 (math-polar (math-sin-raw (math-complex x
))))
1089 ((Math-integer-negp (nth 1 x
))
1090 (math-neg-float (math-sin-raw (math-neg-float x
) (if orgx orgx x
))))
1091 ((math-lessp-float '(float 7 0) x
) ; avoid inf loops due to roundoff
1092 (math-sin-raw (math-mod x
(math-two-pi)) (if orgx orgx x
)))
1093 (t (math-sin-raw-2 x
(if orgx orgx x
)))))
1095 (defun math-cos-raw (x) ; [N N]
1096 (if (eq (car-safe x
) 'polar
)
1097 (math-polar (math-cos-raw (math-complex x
)))
1098 (math-sin-raw (math-sub (math-pi-over-2) x
) x
)))
1100 (defun math-sec-raw (x) ; [N N]
1101 (cond ((eq (car x
) 'cplx
)
1102 (let* ((x (math-mul x
'(float 1 0)))
1103 (expx (math-exp-raw (nth 2 x
)))
1104 (expmx (math-div-float '(float 1 0) expx
))
1105 (sh (math-mul-float (math-sub-float expx expmx
) '(float 5 -
1)))
1106 (ch (math-mul-float (math-add-float expx expmx
) '(float 5 -
1)))
1107 (sc (math-sin-cos-raw (nth 1 x
)))
1109 (math-mul-float (math-sqr (car sc
))
1111 (math-mul-float (math-sqr (cdr sc
))
1113 (and (not (eq (nth 1 d
) 0))
1115 (math-div-float (math-mul-float (cdr sc
) ch
) d
)
1116 (math-div-float (math-mul-float (car sc
) sh
) d
)))))
1117 ((eq (car x
) 'polar
)
1118 (math-polar (math-sec-raw (math-complex x
))))
1120 (let ((cs (math-cos-raw x
)))
1123 (math-div-float '(float 1 0) cs
))))))
1125 (defun math-csc-raw (x) ; [N N]
1126 (cond ((eq (car x
) 'cplx
)
1127 (let* ((x (math-mul x
'(float 1 0)))
1128 (expx (math-exp-raw (nth 2 x
)))
1129 (expmx (math-div-float '(float 1 0) expx
))
1130 (sh (math-mul-float (math-sub-float expx expmx
) '(float 5 -
1)))
1131 (ch (math-mul-float (math-add-float expx expmx
) '(float 5 -
1)))
1132 (sc (math-sin-cos-raw (nth 1 x
)))
1134 (math-mul-float (math-sqr (car sc
))
1136 (math-mul-float (math-sqr (cdr sc
))
1138 (and (not (eq (nth 1 d
) 0))
1140 (math-div-float (math-mul-float (car sc
) ch
) d
)
1141 (math-div-float (math-mul-float (cdr sc
) sh
) d
)))))
1142 ((eq (car x
) 'polar
)
1143 (math-polar (math-csc-raw (math-complex x
))))
1145 (let ((sn (math-sin-raw x
)))
1148 (math-div-float '(float 1 0) sn
))))))
1150 (defun math-cot-raw (x) ; [N N]
1151 (cond ((eq (car x
) 'cplx
)
1152 (let* ((x (math-mul x
'(float 1 0)))
1153 (expx (math-exp-raw (nth 2 x
)))
1154 (expmx (math-div-float '(float 1 0) expx
))
1155 (sh (math-mul-float (math-sub-float expx expmx
) '(float 5 -
1)))
1156 (ch (math-mul-float (math-add-float expx expmx
) '(float 5 -
1)))
1157 (sc (math-sin-cos-raw (nth 1 x
)))
1161 (and (not (eq (nth 1 d
) 0))
1164 (math-mul-float (car sc
) (cdr sc
))
1168 (math-mul-float sh ch
)
1170 ((eq (car x
) 'polar
)
1171 (math-polar (math-cot-raw (math-complex x
))))
1173 (let ((sc (math-sin-cos-raw x
)))
1174 (if (eq (nth 1 (car sc
)) 0)
1175 (math-div (cdr sc
) 0)
1176 (math-div-float (cdr sc
) (car sc
)))))))
1179 ;;; This could use a smarter method: Reduce x as in math-sin-raw, then
1180 ;;; compute either sin(x) or cos(x), whichever is smaller, and compute
1181 ;;; the other using the identity sin(x)^2 + cos(x)^2 = 1.
1182 (defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x))
1183 (cons (math-sin-raw x
) (math-cos-raw x
)))
1185 (defun math-tan-raw (x) ; [N N]
1186 (cond ((eq (car x
) 'cplx
)
1187 (let* ((x (math-mul x
'(float 2 0)))
1188 (expx (math-exp-raw (nth 2 x
)))
1189 (expmx (math-div-float '(float 1 0) expx
))
1190 (sc (math-sin-cos-raw (nth 1 x
)))
1191 (d (math-add-float (cdr sc
)
1192 (math-mul-float (math-add-float expx expmx
)
1194 (and (not (eq (nth 1 d
) 0))
1196 (math-div-float (car sc
) d
)
1197 (math-div-float (math-mul-float (math-sub-float expx
1199 '(float 5 -
1)) d
)))))
1200 ((eq (car x
) 'polar
)
1201 (math-polar (math-tan-raw (math-complex x
))))
1203 (let ((sc (math-sin-cos-raw x
)))
1204 (if (eq (nth 1 (cdr sc
)) 0)
1205 (math-div (car sc
) 0)
1206 (math-div-float (car sc
) (cdr sc
)))))))
1208 (defun math-sin-raw-2 (x orgx
) ; This avoids poss of inf recursion. [F F]
1209 (let ((xmpo2 (math-sub-float (math-pi-over-2) x
)))
1210 (cond ((Math-integer-negp (nth 1 xmpo2
))
1211 (math-neg-float (math-sin-raw-2 (math-sub-float x
(math-pi))
1213 ((math-lessp-float (math-pi-over-4) x
)
1214 (math-cos-raw-2 xmpo2 orgx
))
1215 ((math-lessp-float x
(math-neg (math-pi-over-4)))
1216 (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x
) orgx
)))
1217 ((math-with-extra-prec -
1 (math-nearly-zerop-float x orgx
))
1219 ((math-use-emacs-fn 'sin x
))
1220 (calc-symbolic-mode (signal 'inexact-result nil
))
1221 (t (math-sin-series x
6 4 x
(math-neg-float (math-sqr-float x
)))))))
1223 (defun math-cos-raw-2 (x orgx
) ; [F F]
1224 (cond ((math-with-extra-prec -
1 (math-nearly-zerop-float x orgx
))
1226 ((math-use-emacs-fn 'cos x
))
1227 (calc-symbolic-mode (signal 'inexact-result nil
))
1228 (t (let ((xnegsqr (math-neg-float (math-sqr-float x
))))
1230 (math-add-float '(float 1 0)
1231 (math-mul-float xnegsqr
'(float 5 -
1)))
1232 24 5 xnegsqr xnegsqr
)))))
1234 (defun math-sin-series (sum nfac n x xnegsqr
)
1235 (math-working "sin" sum
)
1236 (let* ((nextx (math-mul-float x xnegsqr
))
1237 (nextsum (math-add-float sum
(math-div-float nextx
1238 (math-float nfac
)))))
1239 (if (math-nearly-equal-float sum nextsum
)
1241 (math-sin-series nextsum
(math-mul nfac
(* n
(1+ n
)))
1242 (+ n
2) nextx xnegsqr
))))
1245 ;;; Inverse sine, cosine, tangent.
1247 (defun calcFunc-arcsin (x) ; [N N] [Public]
1249 ((and (eq x
1) (eq calc-angle-mode
'deg
)) 90)
1250 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) -
90)
1251 (calc-symbolic-mode (signal 'inexact-result nil
))
1253 (math-with-extra-prec 2
1254 (math-from-radians (math-arcsin-raw (math-float x
)))))
1256 (math-make-sdev (calcFunc-arcsin (nth 1 x
))
1260 (math-sub 1 (math-sqr (nth 1 x
))))))))
1262 (math-sort-intv (nth 1 x
)
1263 (calcFunc-arcsin (nth 2 x
))
1264 (calcFunc-arcsin (nth 3 x
))))
1265 ((equal x
'(var nan var-nan
))
1267 (t (calc-record-why 'numberp x
)
1268 (list 'calcFunc-arcsin x
))))
1270 (defun calcFunc-arccos (x) ; [N N] [Public]
1272 ((and (eq x
0) (eq calc-angle-mode
'deg
)) 90)
1273 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) 180)
1274 (calc-symbolic-mode (signal 'inexact-result nil
))
1276 (math-with-extra-prec 2
1277 (math-from-radians (math-arccos-raw (math-float x
)))))
1279 (math-make-sdev (calcFunc-arccos (nth 1 x
))
1283 (math-sub 1 (math-sqr (nth 1 x
))))))))
1285 (math-sort-intv (nth 1 x
)
1286 (calcFunc-arccos (nth 2 x
))
1287 (calcFunc-arccos (nth 3 x
))))
1288 ((equal x
'(var nan var-nan
))
1290 (t (calc-record-why 'numberp x
)
1291 (list 'calcFunc-arccos x
))))
1293 (defun calcFunc-arctan (x) ; [N N] [Public]
1295 ((and (eq x
1) (eq calc-angle-mode
'deg
)) 45)
1296 ((and (eq x -
1) (eq calc-angle-mode
'deg
)) -
45)
1298 (math-with-extra-prec 2
1299 (math-from-radians (math-arctan-raw (math-float x
)))))
1301 (math-make-sdev (calcFunc-arctan (nth 1 x
))
1304 (math-add 1 (math-sqr (nth 1 x
)))))))
1306 (math-sort-intv (nth 1 x
)
1307 (calcFunc-arctan (nth 2 x
))
1308 (calcFunc-arctan (nth 3 x
))))
1309 ((equal x
'(var inf var-inf
))
1310 (math-quarter-circle t
))
1311 ((equal x
'(neg (var inf var-inf
)))
1312 (math-neg (math-quarter-circle t
)))
1313 ((equal x
'(var nan var-nan
))
1315 (t (calc-record-why 'numberp x
)
1316 (list 'calcFunc-arctan x
))))
1318 (defun math-arcsin-raw (x) ; [N N]
1319 (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x
)))))
1320 (if (or (memq (car x
) '(cplx polar
))
1321 (memq (car a
) '(cplx polar
)))
1322 (math-with-extra-prec 2 ; use extra precision for difficult case
1323 (math-mul '(cplx 0 -
1)
1324 (math-ln-raw (math-add (math-mul '(cplx 0 1) x
) a
))))
1325 (math-arctan2-raw x a
))))
1327 (defun math-arccos-raw (x) ; [N N]
1328 (math-sub (math-pi-over-2) (math-arcsin-raw x
)))
1330 (defun math-arctan-raw (x) ; [N N]
1331 (cond ((memq (car x
) '(cplx polar
))
1332 (math-with-extra-prec 2 ; extra-extra
1334 (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x
)))
1335 (math-ln-raw (math-add 1 (math-mul '(cplx 0 -
1) x
))))
1337 ((Math-integer-negp (nth 1 x
))
1338 (math-neg-float (math-arctan-raw (math-neg-float x
))))
1340 ((math-use-emacs-fn 'atan x
))
1341 (calc-symbolic-mode (signal 'inexact-result nil
))
1342 ((math-equal-int x
1) (math-pi-over-4))
1343 ((math-equal-int x -
1) (math-neg (math-pi-over-4)))
1344 ((math-lessp-float '(float 414214 -
6) x
) ; if x > sqrt(2) - 1, reduce
1345 (if (math-lessp-float '(float 1 0) x
)
1346 (math-sub-float (math-mul-float (math-pi) '(float 5 -
1))
1347 (math-arctan-raw (math-div-float '(float 1 0) x
)))
1348 (math-sub-float (math-mul-float (math-pi) '(float 25 -
2))
1349 (math-arctan-raw (math-div-float
1350 (math-sub-float '(float 1 0) x
)
1351 (math-add-float '(float 1 0)
1353 (t (math-arctan-series x
3 x
(math-neg-float (math-sqr-float x
))))))
1355 (defun math-arctan-series (sum n x xnegsqr
)
1356 (math-working "arctan" sum
)
1357 (let* ((nextx (math-mul-float x xnegsqr
))
1358 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1359 (if (math-nearly-equal-float sum nextsum
)
1361 (math-arctan-series nextsum
(+ n
2) nextx xnegsqr
))))
1363 (defun calcFunc-arctan2 (y x
) ; [F R R] [Public]
1366 (math-with-extra-prec 2
1367 (math-from-radians (math-arctan2-raw (math-float y
)
1369 (calc-record-why 'anglep x
)
1370 (list 'calcFunc-arctan2 y x
))
1371 (if (and (or (math-infinitep x
) (math-anglep x
))
1372 (or (math-infinitep y
) (math-anglep y
)))
1387 (calcFunc-arctan2 y x
)
1388 '(var nan var-nan
)))
1389 (calc-record-why 'anglep y
)
1390 (list 'calcFunc-arctan2 y x
))))
1392 (defun math-arctan2-raw (y x
) ; [F R R]
1393 (cond ((math-zerop y
)
1394 (if (math-negp x
) (math-pi)
1395 (if (or (math-floatp x
) (math-floatp y
)) '(float 0 0) 0)))
1399 (math-neg (math-pi-over-2))))
1401 (math-arctan-raw (math-div-float y x
)))
1403 (math-add-float (math-arctan-raw (math-div-float y x
))
1406 (math-sub-float (math-arctan-raw (math-div-float y x
))
1409 (defun calcFunc-arcsincos (x) ; [V N] [Public]
1410 (if (and (Math-vectorp x
)
1412 (calcFunc-arctan2 (nth 2 x
) (nth 1 x
))
1413 (math-reject-arg x
"*Two-element vector expected")))
1417 ;;; Exponential function.
1419 (defun calcFunc-exp (x) ; [N N] [Public]
1421 ((and (memq x
'(1 -
1)) calc-symbolic-mode
)
1422 (if (eq x
1) '(var e var-e
) (math-div 1 '(var e var-e
))))
1424 (math-with-extra-prec 2 (math-exp-raw (math-float x
))))
1425 ((eq (car-safe x
) 'sdev
)
1426 (let ((ex (calcFunc-exp (nth 1 x
))))
1427 (math-make-sdev ex
(math-mul (nth 2 x
) ex
))))
1428 ((eq (car-safe x
) 'intv
)
1429 (math-make-intv (nth 1 x
) (calcFunc-exp (nth 2 x
))
1430 (calcFunc-exp (nth 3 x
))))
1431 ((equal x
'(var inf var-inf
))
1433 ((equal x
'(neg (var inf var-inf
)))
1435 ((equal x
'(var nan var-nan
))
1437 (t (calc-record-why 'numberp x
)
1438 (list 'calcFunc-exp x
))))
1440 (defun calcFunc-expm1 (x) ; [N N] [Public]
1442 ((math-zerop x
) '(float 0 0))
1443 (calc-symbolic-mode (signal 'inexact-result nil
))
1445 (math-with-extra-prec 2
1446 (let ((x (math-float x
)))
1447 (if (and (eq (car x
) 'float
)
1448 (math-lessp-float x
'(float 1 0))
1449 (math-lessp-float '(float -
1 0) x
))
1450 (math-exp-minus-1-raw x
)
1451 (math-add (math-exp-raw x
) -
1)))))
1452 ((eq (car-safe x
) 'sdev
)
1454 (let ((ex (calcFunc-expm1 (nth 1 x
))))
1455 (math-make-sdev ex
(math-mul (nth 2 x
) (math-add ex
1))))
1456 (math-make-sdev (calcFunc-expm1 (nth 1 x
))
1457 (math-mul (nth 2 x
) (calcFunc-exp (nth 1 x
))))))
1458 ((eq (car-safe x
) 'intv
)
1459 (math-make-intv (nth 1 x
)
1460 (calcFunc-expm1 (nth 2 x
))
1461 (calcFunc-expm1 (nth 3 x
))))
1462 ((equal x
'(var inf var-inf
))
1464 ((equal x
'(neg (var inf var-inf
)))
1466 ((equal x
'(var nan var-nan
))
1468 (t (calc-record-why 'numberp x
)
1469 (list 'calcFunc-expm1 x
))))
1471 (defun calcFunc-exp10 (x) ; [N N] [Public]
1474 (math-pow '(float 1 1) x
)))
1476 (defun math-exp-raw (x) ; [N N]
1477 (cond ((math-zerop x
) '(float 1 0))
1478 (calc-symbolic-mode (signal 'inexact-result nil
))
1480 (let ((expx (math-exp-raw (nth 1 x
)))
1481 (sc (math-sin-cos-raw (nth 2 x
))))
1483 (math-mul-float expx
(cdr sc
))
1484 (math-mul-float expx
(car sc
)))))
1485 ((eq (car x
) 'polar
)
1486 (let ((xc (math-complex x
)))
1488 (math-exp-raw (nth 1 xc
))
1489 (math-from-radians (nth 2 xc
)))))
1490 ((math-use-emacs-fn 'exp x
))
1491 ((or (math-lessp-float '(float 5 -
1) x
)
1492 (math-lessp-float x
'(float -
5 -
1)))
1493 (if (math-lessp-float '(float 921035 1) x
)
1495 (if (math-lessp-float x
'(float -
921035 1))
1497 (let* ((two-x (math-mul-float x
'(float 2 0)))
1498 (hint (math-scale-int (nth 1 two-x
) (nth 2 two-x
)))
1499 (hfrac (math-sub-float x
(math-mul-float (math-float hint
)
1501 (math-mul-float (math-ipow (math-sqrt-e) hint
)
1502 (math-add-float '(float 1 0)
1503 (math-exp-minus-1-raw hfrac
)))))
1504 (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x
)))))
1506 (defun math-exp-minus-1-raw (x) ; [F F]
1507 (math-exp-series x
2 3 x x
))
1509 (defun math-exp-series (sum nfac n xpow x
)
1510 (math-working "exp" sum
)
1511 (let* ((nextx (math-mul-float xpow x
))
1512 (nextsum (math-add-float sum
(math-div-float nextx
1513 (math-float nfac
)))))
1514 (if (math-nearly-equal-float sum nextsum
)
1516 (math-exp-series nextsum
(math-mul nfac n
) (1+ n
) nextx x
))))
1522 (defun calcFunc-ln (x) ; [N N] [Public]
1523 (cond ((math-zerop x
)
1524 (if calc-infinite-mode
1525 '(neg (var inf var-inf
))
1526 (math-reject-arg x
"*Logarithm of zero")))
1529 (math-with-extra-prec 2 (math-ln-raw (math-float x
))))
1530 ((eq (car-safe x
) 'sdev
)
1531 (math-make-sdev (calcFunc-ln (nth 1 x
))
1532 (math-div (nth 2 x
) (nth 1 x
))))
1533 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1534 (Math-zerop (nth 2 x
))
1535 (not (math-intv-constp x
))))
1536 (let ((calc-infinite-mode t
))
1537 (math-make-intv (nth 1 x
) (calcFunc-ln (nth 2 x
))
1538 (calcFunc-ln (nth 3 x
)))))
1539 ((equal x
'(var e var-e
))
1541 ((and (eq (car-safe x
) '^
)
1542 (equal (nth 1 x
) '(var e var-e
))
1543 (math-known-realp (nth 2 x
)))
1546 (if (equal x
'(var nan var-nan
))
1548 '(var inf var-inf
)))
1549 (t (calc-record-why 'numberp x
)
1550 (list 'calcFunc-ln x
))))
1552 (defun calcFunc-log10 (x) ; [N N] [Public]
1553 (cond ((math-equal-int x
1)
1554 (if (math-floatp x
) '(float 0 0) 0))
1555 ((and (Math-integerp x
)
1557 (let ((res (math-integer-log x
10)))
1559 (setq x
(cdr res
)))))
1561 ((and (eq (car-safe x
) 'frac
)
1563 (let ((res (math-integer-log (nth 2 x
) 10)))
1565 (setq x
(- (cdr res
))))))
1568 (if calc-infinite-mode
1569 '(neg (var inf var-inf
))
1570 (math-reject-arg x
"*Logarithm of zero")))
1571 (calc-symbolic-mode (signal 'inexact-result nil
))
1573 (math-with-extra-prec 2
1574 (let ((xf (math-float x
)))
1575 (if (eq (nth 1 xf
) 0)
1576 (math-reject-arg x
"*Logarithm of zero"))
1577 (if (Math-integer-posp (nth 1 xf
))
1578 (if (eq (nth 1 xf
) 1) ; log10(1*10^n) = n
1579 (math-float (nth 2 xf
))
1580 (let ((xdigs (1- (math-numdigs (nth 1 xf
)))))
1582 (math-div-float (math-ln-raw-2
1583 (list 'float
(nth 1 xf
) (- xdigs
)))
1585 (math-float (+ (nth 2 xf
) xdigs
)))))
1586 (math-div (calcFunc-ln xf
) (math-ln-10))))))
1587 ((eq (car-safe x
) 'sdev
)
1588 (math-make-sdev (calcFunc-log10 (nth 1 x
))
1590 (math-mul (nth 1 x
) (math-ln-10)))))
1591 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1592 (not (math-intv-constp x
))))
1593 (math-make-intv (nth 1 x
)
1594 (calcFunc-log10 (nth 2 x
))
1595 (calcFunc-log10 (nth 3 x
))))
1597 (if (equal x
'(var nan var-nan
))
1599 '(var inf var-inf
)))
1600 (t (calc-record-why 'numberp x
)
1601 (list 'calcFunc-log10 x
))))
1603 (defun calcFunc-log (x &optional b
) ; [N N N] [Public]
1604 (cond ((or (null b
) (equal b
'(var e var-e
)))
1605 (math-normalize (list 'calcFunc-ln x
)))
1606 ((or (eq b
10) (equal b
'(float 1 1)))
1607 (math-normalize (list 'calcFunc-log10 x
)))
1609 (if calc-infinite-mode
1610 (math-div (calcFunc-ln x
) (calcFunc-ln b
))
1611 (math-reject-arg x
"*Logarithm of zero")))
1613 (if calc-infinite-mode
1614 (math-div (calcFunc-ln x
) (calcFunc-ln b
))
1615 (math-reject-arg b
"*Logarithm of zero")))
1616 ((math-equal-int b
1)
1617 (if calc-infinite-mode
1618 (math-div (calcFunc-ln x
) 0)
1619 (math-reject-arg b
"*Logarithm base one")))
1620 ((math-equal-int x
1)
1621 (if (math-floatp b
) '(float 0 0) 0))
1622 ((and (Math-ratp x
) (Math-ratp b
)
1623 (math-posp x
) (math-posp b
)
1624 (let* ((sign 1) (inv nil
)
1625 (xx (if (Math-lessp 1 x
)
1629 (bb (if (Math-lessp 1 b
)
1631 (setq sign
(- sign
))
1633 (res (if (Math-lessp xx bb
)
1634 (setq inv
(math-integer-log bb xx
))
1635 (math-integer-log xx bb
))))
1638 (math-div 1 (* sign
(cdr res
)))
1639 (* sign
(cdr res
)))))))
1641 (calc-symbolic-mode (signal 'inexact-result nil
))
1642 ((and (Math-numberp x
) (Math-numberp b
))
1643 (math-with-extra-prec 2
1644 (math-div (math-ln-raw (math-float x
))
1645 (math-log-base-raw b
))))
1646 ((and (eq (car-safe x
) 'sdev
)
1648 (math-make-sdev (calcFunc-log (nth 1 x
) b
)
1651 (math-log-base-raw b
)))))
1652 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1653 (not (math-intv-constp x
)))
1655 (math-make-intv (nth 1 x
)
1656 (calcFunc-log (nth 2 x
) b
)
1657 (calcFunc-log (nth 3 x
) b
)))
1658 ((or (eq (car-safe x
) 'intv
) (eq (car-safe b
) 'intv
))
1659 (math-div (calcFunc-ln x
) (calcFunc-ln b
)))
1660 ((or (math-infinitep x
)
1662 (math-div (calcFunc-ln x
) (calcFunc-ln b
)))
1663 (t (if (Math-numberp b
)
1664 (calc-record-why 'numberp x
)
1665 (calc-record-why 'numberp b
))
1666 (list 'calcFunc-log x b
))))
1668 (defun calcFunc-alog (x &optional b
)
1669 (cond ((or (null b
) (equal b
'(var e var-e
)))
1670 (math-normalize (list 'calcFunc-exp x
)))
1671 (t (math-pow b x
))))
1673 (defun calcFunc-ilog (x b
)
1674 (if (and (math-natnump x
) (not (eq x
0))
1675 (math-natnump b
) (not (eq b
0)))
1677 (math-reject-arg x
"*Logarithm base one")
1678 (if (Math-natnum-lessp x b
)
1680 (cdr (math-integer-log x b
))))
1681 (math-floor (calcFunc-log x b
))))
1683 (defun math-integer-log (x b
)
1684 (let ((pows (list b
))
1688 (while (not (Math-lessp x pow
))
1689 (setq pows
(cons pow pows
)
1690 pow
(math-sqr pow
)))
1691 (setq n
(lsh 1 (1- (length pows
)))
1694 (while (and (setq pows
(cdr pows
))
1697 next
(math-mul pow
(car pows
)))
1698 (or (Math-lessp x next
)
1701 (cons (equal pow x
) sum
)))
1704 (defvar math-log-base-cache nil
)
1705 (defun math-log-base-raw (b) ; [N N]
1706 (if (not (and (equal (car math-log-base-cache
) b
)
1707 (eq (nth 1 math-log-base-cache
) calc-internal-prec
)))
1708 (setq math-log-base-cache
(list b calc-internal-prec
1709 (math-ln-raw (math-float b
)))))
1710 (nth 2 math-log-base-cache
))
1712 (defun calcFunc-lnp1 (x) ; [N N] [Public]
1713 (cond ((Math-equal-int x -
1)
1714 (if calc-infinite-mode
1715 '(neg (var inf var-inf
))
1716 (math-reject-arg x
"*Logarithm of zero")))
1718 ((math-zerop x
) '(float 0 0))
1719 (calc-symbolic-mode (signal 'inexact-result nil
))
1721 (math-with-extra-prec 2
1722 (let ((x (math-float x
)))
1723 (if (and (eq (car x
) 'float
)
1724 (math-lessp-float x
'(float 5 -
1))
1725 (math-lessp-float '(float -
5 -
1) x
))
1726 (math-ln-plus-1-raw x
)
1727 (math-ln-raw (math-add-float x
'(float 1 0)))))))
1728 ((eq (car-safe x
) 'sdev
)
1729 (math-make-sdev (calcFunc-lnp1 (nth 1 x
))
1730 (math-div (nth 2 x
) (math-add (nth 1 x
) 1))))
1731 ((and (eq (car-safe x
) 'intv
) (or (Math-posp (nth 2 x
))
1732 (not (math-intv-constp x
))))
1733 (math-make-intv (nth 1 x
)
1734 (calcFunc-lnp1 (nth 2 x
))
1735 (calcFunc-lnp1 (nth 3 x
))))
1737 (if (equal x
'(var nan var-nan
))
1739 '(var inf var-inf
)))
1740 (t (calc-record-why 'numberp x
)
1741 (list 'calcFunc-lnp1 x
))))
1743 (defun math-ln-raw (x) ; [N N] --- must be float format!
1744 (cond ((eq (car-safe x
) 'cplx
)
1746 (math-mul-float (math-ln-raw
1747 (math-add-float (math-sqr-float (nth 1 x
))
1748 (math-sqr-float (nth 2 x
))))
1750 (math-arctan2-raw (nth 2 x
) (nth 1 x
))))
1751 ((eq (car x
) 'polar
)
1752 (math-polar (list 'cplx
1753 (math-ln-raw (nth 1 x
))
1754 (math-to-radians (nth 2 x
)))))
1755 ((Math-equal-int x
1)
1757 (calc-symbolic-mode (signal 'inexact-result nil
))
1758 ((math-posp (nth 1 x
)) ; positive and real
1760 ((math-use-emacs-fn 'log x
))
1762 (let ((xdigs (1- (math-numdigs (nth 1 x
)))))
1763 (math-add-float (math-ln-raw-2 (list 'float
(nth 1 x
) (- xdigs
)))
1764 (math-mul-float (math-float (+ (nth 2 x
) xdigs
))
1767 (math-reject-arg x
"*Logarithm of zero"))
1768 ((eq calc-complex-mode
'polar
) ; negative and real
1770 (list 'cplx
; negative and real
1771 (math-ln-raw (math-neg-float x
))
1773 (t (list 'cplx
; negative and real
1774 (math-ln-raw (math-neg-float x
))
1777 (defun math-ln-raw-2 (x) ; [F F]
1778 (cond ((math-lessp-float '(float 14 -
1) x
)
1779 (math-add-float (math-ln-raw-2 (math-mul-float x
'(float 5 -
1)))
1781 (t ; now .7 < x <= 1.4
1782 (math-ln-raw-3 (math-div-float (math-sub-float x
'(float 1 0))
1783 (math-add-float x
'(float 1 0)))))))
1785 (defun math-ln-raw-3 (x) ; [F F]
1786 (math-mul-float (math-ln-raw-series x
3 x
(math-sqr-float x
))
1789 ;;; Compute ln((1+x)/(1-x))
1790 (defun math-ln-raw-series (sum n x xsqr
)
1791 (math-working "log" sum
)
1792 (let* ((nextx (math-mul-float x xsqr
))
1793 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1794 (if (math-nearly-equal-float sum nextsum
)
1796 (math-ln-raw-series nextsum
(+ n
2) nextx xsqr
))))
1798 (defun math-ln-plus-1-raw (x)
1799 (math-lnp1-series x
2 x
(math-neg x
)))
1801 (defun math-lnp1-series (sum n xpow x
)
1802 (math-working "lnp1" sum
)
1803 (let* ((nextx (math-mul-float xpow x
))
1804 (nextsum (math-add-float sum
(math-div-float nextx
(math-float n
)))))
1805 (if (math-nearly-equal-float sum nextsum
)
1807 (math-lnp1-series nextsum
(1+ n
) nextx x
))))
1809 (defconst math-approx-ln-10
1810 (math-read-number-simple "2.302585092994045684018")
1811 "An approximation for ln(10).")
1813 (math-defcache math-ln-10 math-approx-ln-10
1814 (math-ln-raw-2 '(float 1 1)))
1816 (defconst math-approx-ln-2
1817 (math-read-number-simple "0.693147180559945309417")
1818 "An approximation for ln(2).")
1820 (math-defcache math-ln-2 math-approx-ln-2
1821 (math-ln-raw-3 (math-float '(frac 1 3))))
1825 ;;; Hyperbolic functions.
1827 (defun calcFunc-sinh (x) ; [N N] [Public]
1829 (math-expand-formulas
1831 (list '/ (list '-
(list 'calcFunc-exp x
)
1832 (list 'calcFunc-exp
(list 'neg x
))) 2)))
1834 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1835 (math-with-extra-prec 2
1836 (let ((expx (math-exp-raw (math-float x
))))
1837 (math-mul (math-add expx
(math-div -
1 expx
)) '(float 5 -
1)))))
1838 ((eq (car-safe x
) 'sdev
)
1839 (math-make-sdev (calcFunc-sinh (nth 1 x
))
1840 (math-mul (nth 2 x
) (calcFunc-cosh (nth 1 x
)))))
1842 (math-sort-intv (nth 1 x
)
1843 (calcFunc-sinh (nth 2 x
))
1844 (calcFunc-sinh (nth 3 x
))))
1845 ((or (equal x
'(var inf var-inf
))
1846 (equal x
'(neg (var inf var-inf
)))
1847 (equal x
'(var nan var-nan
)))
1849 (t (calc-record-why 'numberp x
)
1850 (list 'calcFunc-sinh x
))))
1851 (put 'calcFunc-sinh
'math-expandable t
)
1853 (defun calcFunc-cosh (x) ; [N N] [Public]
1855 (math-expand-formulas
1857 (list '/ (list '+ (list 'calcFunc-exp x
)
1858 (list 'calcFunc-exp
(list 'neg x
))) 2)))
1860 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1861 (math-with-extra-prec 2
1862 (let ((expx (math-exp-raw (math-float x
))))
1863 (math-mul (math-add expx
(math-div 1 expx
)) '(float 5 -
1)))))
1864 ((eq (car-safe x
) 'sdev
)
1865 (math-make-sdev (calcFunc-cosh (nth 1 x
))
1867 (calcFunc-sinh (nth 1 x
)))))
1868 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1869 (setq x
(math-abs x
))
1870 (math-sort-intv (nth 1 x
)
1871 (calcFunc-cosh (nth 2 x
))
1872 (calcFunc-cosh (nth 3 x
))))
1873 ((or (equal x
'(var inf var-inf
))
1874 (equal x
'(neg (var inf var-inf
)))
1875 (equal x
'(var nan var-nan
)))
1877 (t (calc-record-why 'numberp x
)
1878 (list 'calcFunc-cosh x
))))
1879 (put 'calcFunc-cosh
'math-expandable t
)
1881 (defun calcFunc-tanh (x) ; [N N] [Public]
1883 (math-expand-formulas
1885 (let ((expx (list 'calcFunc-exp x
))
1886 (expmx (list 'calcFunc-exp
(list 'neg x
))))
1888 (list '/ (list '- expx expmx
) (list '+ expx expmx
))))))
1890 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1891 (math-with-extra-prec 2
1892 (let* ((expx (calcFunc-exp (math-float x
)))
1893 (expmx (math-div 1 expx
)))
1894 (math-div (math-sub expx expmx
)
1895 (math-add expx expmx
)))))
1896 ((eq (car-safe x
) 'sdev
)
1897 (math-make-sdev (calcFunc-tanh (nth 1 x
))
1899 (math-sqr (calcFunc-cosh (nth 1 x
))))))
1901 (math-sort-intv (nth 1 x
)
1902 (calcFunc-tanh (nth 2 x
))
1903 (calcFunc-tanh (nth 3 x
))))
1904 ((equal x
'(var inf var-inf
))
1906 ((equal x
'(neg (var inf var-inf
)))
1908 ((equal x
'(var nan var-nan
))
1910 (t (calc-record-why 'numberp x
)
1911 (list 'calcFunc-tanh x
))))
1912 (put 'calcFunc-tanh
'math-expandable t
)
1914 (defun calcFunc-sech (x) ; [N N] [Public]
1916 (math-expand-formulas
1918 (list '/ 2 (list '+ (list 'calcFunc-exp x
)
1919 (list 'calcFunc-exp
(list 'neg x
))))))
1921 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1922 (math-with-extra-prec 2
1923 (let ((expx (math-exp-raw (math-float x
))))
1924 (math-div '(float 2 0) (math-add expx
(math-div 1 expx
))))))
1925 ((eq (car-safe x
) 'sdev
)
1926 (math-make-sdev (calcFunc-sech (nth 1 x
))
1928 (math-mul (calcFunc-sech (nth 1 x
))
1929 (calcFunc-tanh (nth 1 x
))))))
1930 ((and (eq (car x
) 'intv
) (math-intv-constp x
))
1931 (setq x
(math-abs x
))
1932 (math-sort-intv (nth 1 x
)
1933 (calcFunc-sech (nth 2 x
))
1934 (calcFunc-sech (nth 3 x
))))
1935 ((or (equal x
'(var inf var-inf
))
1936 (equal x
'(neg (var inf var-inf
))))
1938 ((equal x
'(var nan var-nan
))
1940 (t (calc-record-why 'numberp x
)
1941 (list 'calcFunc-sech x
))))
1942 (put 'calcFunc-sech
'math-expandable t
)
1944 (defun calcFunc-csch (x) ; [N N] [Public]
1945 (cond ((eq x
0) (math-div 1 0))
1946 (math-expand-formulas
1948 (list '/ 2 (list '-
(list 'calcFunc-exp x
)
1949 (list 'calcFunc-exp
(list 'neg x
))))))
1951 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1952 (math-with-extra-prec 2
1953 (let ((expx (math-exp-raw (math-float x
))))
1954 (math-div '(float 2 0) (math-add expx
(math-div -
1 expx
))))))
1955 ((eq (car-safe x
) 'sdev
)
1956 (math-make-sdev (calcFunc-csch (nth 1 x
))
1958 (math-mul (calcFunc-csch (nth 1 x
))
1959 (calcFunc-coth (nth 1 x
))))))
1961 (if (and (Math-negp (nth 2 x
))
1962 (Math-posp (nth 3 x
)))
1963 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
1964 (math-sort-intv (nth 1 x
)
1965 (calcFunc-csch (nth 2 x
))
1966 (calcFunc-csch (nth 3 x
)))))
1967 ((or (equal x
'(var inf var-inf
))
1968 (equal x
'(neg (var inf var-inf
))))
1970 ((equal x
'(var nan var-nan
))
1972 (t (calc-record-why 'numberp x
)
1973 (list 'calcFunc-csch x
))))
1974 (put 'calcFunc-csch
'math-expandable t
)
1976 (defun calcFunc-coth (x) ; [N N] [Public]
1977 (cond ((eq x
0) (math-div 1 0))
1978 (math-expand-formulas
1980 (let ((expx (list 'calcFunc-exp x
))
1981 (expmx (list 'calcFunc-exp
(list 'neg x
))))
1983 (list '/ (list '+ expx expmx
) (list '- expx expmx
))))))
1985 (if calc-symbolic-mode
(signal 'inexact-result nil
))
1986 (math-with-extra-prec 2
1987 (let* ((expx (calcFunc-exp (math-float x
)))
1988 (expmx (math-div 1 expx
)))
1989 (math-div (math-add expx expmx
)
1990 (math-sub expx expmx
)))))
1991 ((eq (car-safe x
) 'sdev
)
1992 (math-make-sdev (calcFunc-coth (nth 1 x
))
1994 (math-sqr (calcFunc-sinh (nth 1 x
))))))
1996 (if (and (Math-negp (nth 2 x
))
1997 (Math-posp (nth 3 x
)))
1998 '(intv 3 (neg (var inf var-inf
)) (var inf var-inf
))
1999 (math-sort-intv (nth 1 x
)
2000 (calcFunc-coth (nth 2 x
))
2001 (calcFunc-coth (nth 3 x
)))))
2002 ((equal x
'(var inf var-inf
))
2004 ((equal x
'(neg (var inf var-inf
)))
2006 ((equal x
'(var nan var-nan
))
2008 (t (calc-record-why 'numberp x
)
2009 (list 'calcFunc-coth x
))))
2010 (put 'calcFunc-coth
'math-expandable t
)
2012 (defun calcFunc-arcsinh (x) ; [N N] [Public]
2014 (math-expand-formulas
2016 (list 'calcFunc-ln
(list '+ x
(list 'calcFunc-sqrt
2017 (list '+ (list '^ x
2) 1))))))
2019 (if calc-symbolic-mode
(signal 'inexact-result nil
))
2020 (math-with-extra-prec 2
2021 (math-ln-raw (math-add x
(math-sqrt-raw (math-add (math-sqr x
)
2023 ((eq (car-safe x
) 'sdev
)
2024 (math-make-sdev (calcFunc-arcsinh (nth 1 x
))
2027 (math-add (math-sqr (nth 1 x
)) 1)))))
2029 (math-sort-intv (nth 1 x
)
2030 (calcFunc-arcsinh (nth 2 x
))
2031 (calcFunc-arcsinh (nth 3 x
))))
2032 ((or (equal x
'(var inf var-inf
))
2033 (equal x
'(neg (var inf var-inf
)))
2034 (equal x
'(var nan var-nan
)))
2036 (t (calc-record-why 'numberp x
)
2037 (list 'calcFunc-arcsinh x
))))
2038 (put 'calcFunc-arcsinh
'math-expandable t
)
2040 (defun calcFunc-arccosh (x) ; [N N] [Public]
2042 ((and (eq x -
1) calc-symbolic-mode
)
2044 ((and (eq x
0) calc-symbolic-mode
)
2045 (math-div (math-mul '(var pi var-pi
) '(var i var-i
)) 2))
2046 (math-expand-formulas
2048 (list 'calcFunc-ln
(list '+ x
(list 'calcFunc-sqrt
2049 (list '-
(list '^ x
2) 1))))))
2051 (if calc-symbolic-mode
(signal 'inexact-result nil
))
2052 (if (Math-equal-int x -
1)
2053 (math-imaginary (math-pi))
2054 (math-with-extra-prec 2
2055 (if (or t
; need to do this even in the real case!
2056 (memq (car-safe x
) '(cplx polar
)))
2057 (let ((xp1 (math-add 1 x
))) ; this gets the branch cuts right
2059 (math-add x
(math-mul xp1
2066 (math-add x
(math-sqrt-raw (math-add (math-sqr x
)
2067 '(float -
1 0)))))))))
2068 ((eq (car-safe x
) 'sdev
)
2069 (math-make-sdev (calcFunc-arccosh (nth 1 x
))
2072 (math-add (math-sqr (nth 1 x
)) -
1)))))
2074 (math-sort-intv (nth 1 x
)
2075 (calcFunc-arccosh (nth 2 x
))
2076 (calcFunc-arccosh (nth 3 x
))))
2077 ((or (equal x
'(var inf var-inf
))
2078 (equal x
'(neg (var inf var-inf
)))
2079 (equal x
'(var nan var-nan
)))
2081 (t (calc-record-why 'numberp x
)
2082 (list 'calcFunc-arccosh x
))))
2083 (put 'calcFunc-arccosh
'math-expandable t
)
2085 (defun calcFunc-arctanh (x) ; [N N] [Public]
2087 ((and (Math-equal-int x
1) calc-infinite-mode
)
2089 ((and (Math-equal-int x -
1) calc-infinite-mode
)
2090 '(neg (var inf var-inf
)))
2091 (math-expand-formulas
2093 (list 'calcFunc-ln
(list '+ 1 x
))
2094 (list 'calcFunc-ln
(list '-
1 x
))) 2))
2096 (if calc-symbolic-mode
(signal 'inexact-result nil
))
2097 (math-with-extra-prec 2
2098 (if (or (memq (car-safe x
) '(cplx polar
))
2100 (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x
))
2101 (math-ln-raw (math-sub '(float 1 0) x
)))
2103 (if (and (math-equal-int x
1) calc-infinite-mode
)
2105 (if (and (math-equal-int x -
1) calc-infinite-mode
)
2106 '(neg (var inf var-inf
))
2107 (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x
)
2110 ((eq (car-safe x
) 'sdev
)
2111 (math-make-sdev (calcFunc-arctanh (nth 1 x
))
2113 (math-sub 1 (math-sqr (nth 1 x
))))))
2115 (math-sort-intv (nth 1 x
)
2116 (calcFunc-arctanh (nth 2 x
))
2117 (calcFunc-arctanh (nth 3 x
))))
2118 ((equal x
'(var nan var-nan
))
2120 (t (calc-record-why 'numberp x
)
2121 (list 'calcFunc-arctanh x
))))
2122 (put 'calcFunc-arctanh
'math-expandable t
)
2125 ;;; Convert A from HMS or degrees to radians.
2126 (defun calcFunc-rad (a) ; [R R] [Public]
2127 (cond ((or (Math-numberp a
)
2129 (math-with-extra-prec 2
2130 (math-mul a
(math-pi-over-180))))
2132 (math-from-hms a
'rad
))
2134 (math-make-sdev (calcFunc-rad (nth 1 a
))
2135 (calcFunc-rad (nth 2 a
))))
2136 (math-expand-formulas
2137 (math-div (math-mul a
'(var pi var-pi
)) 180))
2138 ((math-infinitep a
) a
)
2139 (t (list 'calcFunc-rad a
))))
2140 (put 'calcFunc-rad
'math-expandable t
)
2142 ;;; Convert A from HMS or radians to degrees.
2143 (defun calcFunc-deg (a) ; [R R] [Public]
2144 (cond ((or (Math-numberp a
)
2146 (math-with-extra-prec 2
2147 (math-div a
(math-pi-over-180))))
2149 (math-from-hms a
'deg
))
2151 (math-make-sdev (calcFunc-deg (nth 1 a
))
2152 (calcFunc-deg (nth 2 a
))))
2153 (math-expand-formulas
2154 (math-div (math-mul 180 a
) '(var pi var-pi
)))
2155 ((math-infinitep a
) a
)
2156 (t (list 'calcFunc-deg a
))))
2157 (put 'calcFunc-deg
'math-expandable t
)
2159 (provide 'calc-math
)
2161 ;; arch-tag: c7367e8e-d0b8-4f70-8577-2fb3f31dbb4c
2162 ;;; calc-math.el ends here