Newsgroups: comp.sources.misc
From: daveg@synaptics.com (David Gillespie)
Subject:  v24i069:  gnucalc - GNU Emacs Calculator, v2.00, Part21/56
Message-ID: <1991Oct31.072644.17973@sparky.imd.sterling.com>
X-Md4-Signature: 935912531147058885be56ee6ad9de74
Date: Thu, 31 Oct 1991 07:26:44 GMT
Approved: kent@sparky.imd.sterling.com

Submitted-by: daveg@synaptics.com (David Gillespie)
Posting-number: Volume 24, Issue 69
Archive-name: gnucalc/part21
Environment: Emacs
Supersedes: gmcalc: Volume 13, Issue 27-45

---- Cut Here and unpack ----
#!/bin/sh
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc-mat.el continued
#
if test ! -r _shar_seq_.tmp; then
	echo 'Please unpack part 1 first!'
	exit 1
fi
(read Scheck
 if test "$Scheck" != 21; then
	echo Please unpack part "$Scheck" next!
	exit 1
 else
	exit 0
 fi
) < _shar_seq_.tmp || exit 1
if test ! -f _shar_wnt_.tmp; then
	echo 'x - still skipping calc-mat.el'
else
echo 'x - continuing file calc-mat.el'
sed 's/^X//' << 'SHAR_EOF' >> 'calc-mat.el' &&
X					    (nth j (nth k lu))))
X		k (1+ k)))
X	(setcar (nthcdr j (nth i lu)) sum)
X	(let ((dum (math-abs-approx sum)))
X	  (if (Math-lessp big dum)
X	      (setq big dum
X		    imax i)))
X	(setq i (1+ i)))
X      (if (> imax j)
X	  (setq lu (math-swap-rows lu j imax)
X		d (- d)))
X      (setq index (cons imax index))
X      (let ((pivot (nth j (nth j lu))))
X	(if (math-zerop pivot)
X	    (throw 'singular nil)
X	  (setq i j)
X	  (while (<= (setq i (1+ i)) n)
X	    (setcar (nthcdr j (nth i lu))
X		    (math-div (nth j (nth i lu)) pivot)))))
X      (setq j (1+ j)))
X    (list lu (nreverse index) d))
)
X
(defun math-swap-rows (m r1 r2)
X  (or (= r1 r2)
X      (let* ((r1prev (nthcdr (1- r1) m))
X	     (row1 (cdr r1prev))
X	     (r2prev (nthcdr (1- r2) m))
X	     (row2 (cdr r2prev))
X	     (r2next (cdr row2)))
X	(setcdr r2prev row1)
X	(setcdr r1prev row2)
X	(setcdr row2 (cdr row1))
X	(setcdr row1 r2next)))
X  m
)
X
X
(defun math-lud-solve (lud b &optional need)
X  (if lud
X      (let* ((x (math-copy-matrix b))
X	     (n (1- (length x)))
X	     (m (1- (length (nth 1 x))))
X	     (lu (car lud))
X	     (col 1)
X	     i j ip ii index sum)
X	(while (<= col m)
X	  (math-working "LUD solver step" col)
X	  (setq i 1
X		ii nil
X		index (nth 1 lud))
X	  (while (<= i n)
X	    (setq ip (car index)
X		  index (cdr index)
X		  sum (nth col (nth ip x)))
X	    (setcar (nthcdr col (nth ip x)) (nth col (nth i x)))
X	    (if (null ii)
X		(or (math-zerop sum)
X		    (setq ii i))
X	      (setq j ii)
X	      (while (< j i)
X		(setq sum (math-sub sum (math-mul (nth j (nth i lu))
X						  (nth col (nth j x))))
X		      j (1+ j))))
X	    (setcar (nthcdr col (nth i x)) sum)
X	    (setq i (1+ i)))
X	  (while (>= (setq i (1- i)) 1)
X	    (setq sum (nth col (nth i x))
X		  j i)
X	    (while (<= (setq j (1+ j)) n)
X	      (setq sum (math-sub sum (math-mul (nth j (nth i lu))
X						(nth col (nth j x))))))
X	    (setcar (nthcdr col (nth i x))
X		    (math-div sum (nth i (nth i lu)))))
X	  (setq col (1+ col)))
X	x)
X    (and need
X	 (math-reject-arg need "*Singular matrix")))
)
X
(defun calcFunc-lud (m)
X  (if (math-square-matrixp m)
X      (or (math-with-extra-prec 2
X	    (let ((lud (math-matrix-lud m)))
X	      (and lud
X		   (let* ((lmat (math-copy-matrix (car lud)))
X			  (umat (math-copy-matrix (car lud)))
X			  (n (1- (length (car lud))))
X			  (perm (calcFunc-idn 1 n))
X			  i (j 1))
X		     (while (<= j n)
X		       (setq i 1)
X		       (while (< i j)
X			 (setcar (nthcdr j (nth i lmat)) 0)
X			 (setq i (1+ i)))
X		       (setcar (nthcdr j (nth j lmat)) 1)
X		       (while (<= (setq i (1+ i)) n)
X			 (setcar (nthcdr j (nth i umat)) 0))
X		       (setq j (1+ j)))
X		     (while (>= (setq j (1- j)) 1)
X		       (let ((pos (nth (1- j) (nth 1 lud))))
X			 (or (= pos j)
X			     (setq perm (math-swap-rows perm j pos)))))
X		     (list 'vec perm lmat umat)))))
X	  (math-reject-arg m "*Singular matrix"))
X    (math-reject-arg m 'square-matrixp))
)
X
SHAR_EOF
echo 'File calc-mat.el is complete' &&
chmod 0644 calc-mat.el ||
echo 'restore of calc-mat.el failed'
Wc_c="`wc -c < 'calc-mat.el'`"
test 10372 -eq "$Wc_c" ||
	echo 'calc-mat.el: original size 10372, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= calc-math.el ==============
if test -f 'calc-math.el' -a X"$1" != X"-c"; then
	echo 'x - skipping calc-math.el (File already exists)'
	rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting calc-math.el (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'calc-math.el' &&
;; Calculator for GNU Emacs, part II [calc-math.el]
;; Copyright (C) 1990, 1991 Free Software Foundation, Inc.
;; Written by Dave Gillespie, daveg@csvax.cs.caltech.edu.
X
;; This file is part of GNU Emacs.
X
;; GNU Emacs is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY.  No author or distributor
;; accepts responsibility to anyone for the consequences of using it
;; or for whether it serves any particular purpose or works at all,
;; unless he says so in writing.  Refer to the GNU Emacs General Public
;; License for full details.
X
;; Everyone is granted permission to copy, modify and redistribute
;; GNU Emacs, but only under the conditions described in the
;; GNU Emacs General Public License.   A copy of this license is
;; supposed to have been given to you along with GNU Emacs so you
;; can know your rights and responsibilities.  It should be in a
;; file named COPYING.  Among other things, the copyright notice
;; and this notice must be preserved on all copies.
X
X
X
;; This file is autoloaded from calc-ext.el.
(require 'calc-ext)
X
(require 'calc-macs)
X
(defun calc-Need-calc-math () nil)
X
X
(defun calc-sqrt (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (if (calc-is-inverse)
X       (calc-unary-op "^2" 'calcFunc-sqr arg)
X     (calc-unary-op "sqrt" 'calcFunc-sqrt arg)))
)
X
(defun calc-isqrt (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (if (calc-is-inverse)
X       (calc-unary-op "^2" 'calcFunc-sqr arg)
X     (calc-unary-op "isqt" 'calcFunc-isqrt arg)))
)
X
X
(defun calc-hypot (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (calc-binary-op "hypt" 'calcFunc-hypot arg))
)
X
(defun calc-ln (arg)
X  (interactive "P")
X  (calc-invert-func)
X  (calc-exp arg)
)
X
(defun calc-log10 (arg)
X  (interactive "P")
X  (calc-hyperbolic-func)
X  (calc-ln arg)
)
X
(defun calc-log (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (if (calc-is-inverse)
X       (calc-binary-op "alog" 'calcFunc-alog arg)
X     (calc-binary-op "log" 'calcFunc-log arg)))
)
X
(defun calc-ilog (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (if (calc-is-inverse)
X       (calc-binary-op "alog" 'calcFunc-alog arg)
X     (calc-binary-op "ilog" 'calcFunc-ilog arg)))
)
X
(defun calc-lnp1 (arg)
X  (interactive "P")
X  (calc-invert-func)
X  (calc-expm1 arg)
)
X
(defun calc-exp (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (if (calc-is-hyperbolic)
X       (if (calc-is-inverse)
X	   (calc-unary-op "lg10" 'calcFunc-log10 arg)
X	 (calc-unary-op "10^" 'calcFunc-exp10 arg))
X     (if (calc-is-inverse)
X	 (calc-unary-op "ln" 'calcFunc-ln arg)
X       (calc-unary-op "exp" 'calcFunc-exp arg))))
)
X
(defun calc-expm1 (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (if (calc-is-inverse)
X       (calc-unary-op "ln+1" 'calcFunc-lnp1 arg)
X     (calc-unary-op "ex-1" 'calcFunc-expm1 arg)))
)
X
(defun calc-pi ()
X  (interactive)
X  (calc-slow-wrapper
X   (if (calc-is-inverse)
X       (if (calc-is-hyperbolic)
X	   (if calc-symbolic-mode
X	       (calc-pop-push-record 0 "phi" '(var phi var-phi))
X	     (calc-pop-push-record 0 "phi" (math-phi)))
X	 (if calc-symbolic-mode
X	     (calc-pop-push-record 0 "gmma" '(var gamma var-gamma))
X	   (calc-pop-push-record 0 "gmma" (math-gamma-const))))
X     (if (calc-is-hyperbolic)
X	 (if calc-symbolic-mode
X	     (calc-pop-push-record 0 "e" '(var e var-e))
X	   (calc-pop-push-record 0 "e" (math-e)))
X       (if calc-symbolic-mode
X	   (calc-pop-push-record 0 "pi" '(var pi var-pi))
X	 (calc-pop-push-record 0 "pi" (math-pi))))))
)
X
(defun calc-sin (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (if (calc-is-hyperbolic)
X       (if (calc-is-inverse)
X	   (calc-unary-op "asnh" 'calcFunc-arcsinh arg)
X	 (calc-unary-op "sinh" 'calcFunc-sinh arg))
X     (if (calc-is-inverse)
X	 (calc-unary-op "asin" 'calcFunc-arcsin arg)
X       (calc-unary-op "sin" 'calcFunc-sin arg))))
)
X
(defun calc-arcsin (arg)
X  (interactive "P")
X  (calc-invert-func)
X  (calc-sin arg)
)
X
(defun calc-sinh (arg)
X  (interactive "P")
X  (calc-hyperbolic-func)
X  (calc-sin arg)
)
X
(defun calc-arcsinh (arg)
X  (interactive "P")
X  (calc-invert-func)
X  (calc-hyperbolic-func)
X  (calc-sin arg)
)
X
(defun calc-cos (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (if (calc-is-hyperbolic)
X       (if (calc-is-inverse)
X	   (calc-unary-op "acsh" 'calcFunc-arccosh arg)
X	 (calc-unary-op "cosh" 'calcFunc-cosh arg))
X     (if (calc-is-inverse)
X	 (calc-unary-op "acos" 'calcFunc-arccos arg)
X       (calc-unary-op "cos" 'calcFunc-cos arg))))
)
X
(defun calc-arccos (arg)
X  (interactive "P")
X  (calc-invert-func)
X  (calc-cos arg)
)
X
(defun calc-cosh (arg)
X  (interactive "P")
X  (calc-hyperbolic-func)
X  (calc-cos arg)
)
X
(defun calc-arccosh (arg)
X  (interactive "P")
X  (calc-invert-func)
X  (calc-hyperbolic-func)
X  (calc-cos arg)
)
X
(defun calc-sincos ()
X  (interactive)
X  (calc-slow-wrapper
X   (if (calc-is-inverse)
X       (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1)))
X     (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1)))))
)
X
(defun calc-tan (arg)
X  (interactive "P")
X  (calc-slow-wrapper
X   (if (calc-is-hyperbolic)
X       (if (calc-is-inverse)
X	   (calc-unary-op "atnh" 'calcFunc-arctanh arg)
X	 (calc-unary-op "tanh" 'calcFunc-tanh arg))
X     (if (calc-is-inverse)
X	 (calc-unary-op "atan" 'calcFunc-arctan arg)
X       (calc-unary-op "tan" 'calcFunc-tan arg))))
)
X
(defun calc-arctan (arg)
X  (interactive "P")
X  (calc-invert-func)
X  (calc-tan arg)
)
X
(defun calc-tanh (arg)
X  (interactive "P")
X  (calc-hyperbolic-func)
X  (calc-tan arg)
)
X
(defun calc-arctanh (arg)
X  (interactive "P")
X  (calc-invert-func)
X  (calc-hyperbolic-func)
X  (calc-tan arg)
)
X
(defun calc-arctan2 ()
X  (interactive)
X  (calc-slow-wrapper
X   (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2))))
)
X
(defun calc-conj (arg)
X  (interactive "P")
X  (calc-wrapper
X   (calc-unary-op "conj" 'calcFunc-conj arg))
)
X
(defun calc-imaginary ()
X  (interactive)
X  (calc-slow-wrapper
X   (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1))))
)
X
X
X
(defun calc-to-degrees (arg)
X  (interactive "P")
X  (calc-wrapper
X   (calc-unary-op ">deg" 'calcFunc-deg arg))
)
X
(defun calc-to-radians (arg)
X  (interactive "P")
X  (calc-wrapper
X   (calc-unary-op ">rad" 'calcFunc-rad arg))
)
X
X
(defun calc-degrees-mode ()
X  (interactive)
X  (calc-wrapper
X   (calc-change-mode 'calc-angle-mode 'deg)
X   (message "Angles measured in degrees."))
)
X
(defun calc-radians-mode ()
X  (interactive)
X  (calc-wrapper
X   (calc-change-mode 'calc-angle-mode 'rad)
X   (message "Angles measured in radians."))
)
X
X
;;; Compute the integer square-root floor(sqrt(A)).  A > 0.  [I I] [Public]
;;; This method takes advantage of the fact that Newton's method starting
;;; with an overestimate always works, even using truncating integer division!
(defun math-isqrt (a)
X  (cond ((Math-zerop a) a)
X	((not (math-natnump a))
X	 (math-reject-arg a 'natnump))
X	((integerp a)
X	 (math-isqrt-small a))
X	(t
X	 (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a)))))))
)
X
(defun calcFunc-isqrt (a)
X  (if (math-realp a)
X      (math-isqrt (math-floor a))
X    (math-floor (math-sqrt a)))
)
X
X
;;; This returns (flag . result) where the flag is T if A is a perfect square.
(defun math-isqrt-bignum (a)   ; [P.l L]
X  (let ((len (length a)))
X    (if (= (% len 2) 0)
X	(let* ((top (nthcdr (- len 2) a)))
X	  (math-isqrt-bignum-iter
X	   a
X	   (math-scale-bignum-3
X	    (math-bignum-big
X	     (1+ (math-isqrt-small
X		  (+ (* (nth 1 top) 1000) (car top)))))
X	    (1- (/ len 2)))))
X      (let* ((top (nth (1- len) a)))
X	(math-isqrt-bignum-iter
X	 a
X	 (math-scale-bignum-3
X	  (list (1+ (math-isqrt-small top)))
X	  (/ len 2))))))
)
X
(defun math-isqrt-bignum-iter (a guess)   ; [l L l]
X  (math-working "isqrt" (cons 'bigpos guess))
X  (let* ((q (math-div-bignum a guess))
X	 (s (math-add-bignum (car q) guess))
X	 (g2 (math-div2-bignum s))
X	 (comp (math-compare-bignum g2 guess)))
X    (if (< comp 0)
X	(math-isqrt-bignum-iter a g2)
X      (cons (and (= comp 0)
X		 (math-zerop-bignum (cdr q))
X		 (= (% (car s) 2) 0))
X	    guess)))
)
X
(defun math-zerop-bignum (a)
X  (and (eq (car a) 0)
X       (progn
X	 (while (eq (car (setq a (cdr a))) 0))
X	 (null a)))
)
X
(defun math-scale-bignum-3 (a n)   ; [L L S]
X  (while (> n 0)
X    (setq a (cons 0 a)
X	  n (1- n)))
X  a
)
X
(defun math-isqrt-small (a)   ; A > 0.  [S S]
X  (let ((g (cond ((>= a 10000) 1000)
X		 ((>= a 100) 100)
X		 (t 10)))
X	g2)
X    (while (< (setq g2 (/ (+ g (/ a g)) 2)) g)
X      (setq g g2))
X    g)
)
X
X
X
X
;;; Compute the square root of a number.
;;; [T N] if possible, else [F N] if possible, else [C N].  [Public]
(defun math-sqrt (a)
X  (or
X   (and (Math-zerop a) a)
X   (and (math-known-nonposp a)
X	(math-imaginary (math-sqrt (math-neg a))))
X   (and (integerp a)
X	(let ((sqrt (math-isqrt-small a)))
X	  (if (= (* sqrt sqrt) a)
X	      sqrt
X	    (if calc-symbolic-mode
X		(list 'calcFunc-sqrt a)
X	      (math-sqrt-float (math-float a) (math-float sqrt))))))
X   (and (eq (car-safe a) 'bigpos)
X	(let* ((res (math-isqrt-bignum (cdr a)))
X	       (sqrt (math-normalize (cons 'bigpos (cdr res)))))
X	  (if (car res)
X	      sqrt
X	    (if calc-symbolic-mode
X		(list 'calcFunc-sqrt a)
X	      (math-sqrt-float (math-float a) (math-float sqrt))))))
X   (and (eq (car-safe a) 'frac)
X	(let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a)))))
X	       (num-sqrt (math-normalize (cons 'bigpos (cdr num-res))))
X	       (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a)))))
X	       (den-sqrt (math-normalize (cons 'bigpos (cdr den-res)))))
X	  (if (and (car num-res) (car den-res))
X	      (list 'frac num-sqrt den-sqrt)
X	    (if calc-symbolic-mode
X		(if (or (car num-res) (car den-res))
X		    (math-div (if (car num-res)
X				  num-sqrt (list 'calcFunc-sqrt (nth 1 a)))
X			      (if (car den-res)
X				  den-sqrt (list 'calcFunc-sqrt (nth 2 a))))
X		  (list 'calcFunc-sqrt a))
X	      (math-sqrt-float (math-float a)
X			       (math-div (math-float num-sqrt) den-sqrt))))))
X   (and (eq (car-safe a) 'float)
X	(if calc-symbolic-mode
X	    (if (= (% (nth 2 a) 2) 0)
X		(let ((res (math-isqrt-bignum
X			    (cdr (Math-bignum-test (nth 1 a))))))
X		  (if (car res)
X		      (math-make-float (math-normalize
X					(cons 'bigpos (cdr res)))
X				       (/ (nth 2 a) 2))
X		    (signal 'inexact-result nil)))
X	      (signal 'inexact-result nil))
X	  (math-sqrt-float a)))
X   (and (eq (car-safe a) 'cplx)
X	(math-with-extra-prec 2
X	  (let* ((d (math-abs a))
X		 (imag (math-sqrt (math-mul (math-sub d (nth 1 a))
X					    '(float 5 -1)))))
X	    (list 'cplx
X		  (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1)))
X		  (if (math-negp (nth 2 a)) (math-neg imag) imag)))))
X   (and (eq (car-safe a) 'polar)
X	(list 'polar
X	      (math-sqrt (nth 1 a))
X	      (math-mul (nth 2 a) '(float 5 -1))))
X   (and (eq (car-safe a) 'sdev)
X	(let ((sqrt (math-sqrt (nth 1 a))))
X	  (math-make-sdev sqrt
X			  (math-div (nth 2 a) (math-mul sqrt 2)))))
X   (and (eq (car-safe a) 'intv)
X	(not (math-negp (nth 2 a)))
X	(math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a))))
X   (and (eq (car-safe a) '*)
X	(or (math-known-nonnegp (nth 1 a))
X	    (math-known-nonnegp (nth 2 a)))
X	(math-mul (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
X   (and (eq (car-safe a) '/)
X	(or (and (math-known-nonnegp (nth 2 a))
X		 (math-div (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
X	    (and (math-known-nonnegp (nth 1 a))
X		 (not (math-equal-int (nth 1 a) 1))
X		 (math-mul (math-sqrt (nth 1 a))
X			   (math-sqrt (math-div 1 (nth 2 a)))))))
X   (and (eq (car-safe a) '^)
X	(math-known-evenp (nth 2 a))
X	(math-known-realp (nth 1 a))
X	(math-abs (math-pow (nth 1 a) (math-div (nth 2 a) 2))))
X   (let ((inf (math-infinitep a)))
X     (and inf
X	  (math-mul (math-sqrt (math-infinite-dir a inf)) inf)))
X   (progn
X     (calc-record-why 'numberp a)
X     (list 'calcFunc-sqrt a)))
)
(fset 'calcFunc-sqrt (symbol-function 'math-sqrt))
X
(defun math-infinite-dir (a &optional inf)
X  (or inf (setq inf (math-infinitep a)))
X  (math-normalize (math-expr-subst a inf 1))
)
X
(defun math-sqrt-float (a &optional guess)   ; [F F F]
X  (if calc-symbolic-mode
X      (signal 'inexact-result nil)
X    (math-with-extra-prec 1 (math-sqrt-raw a guess)))
)
X
(defun math-sqrt-raw (a &optional guess)   ; [F F F]
X  (if (not (Math-posp a))
X      (math-sqrt a)
X    (if (null guess)
X	(let ((ldiff (- (math-numdigs (nth 1 a)) 6)))
X	  (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff)))
X	  (setq guess (math-make-float (math-isqrt-small
X					(math-scale-int (nth 1 a) (- ldiff)))
X				       (/ (+ (nth 2 a) ldiff) 2)))))
X    (math-sqrt-float-iter a guess))
)
X
(defun math-sqrt-float-iter (a guess)   ; [F F F]
X  (math-working "sqrt" guess)
X  (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess))
X			    '(float 5 -1))))
X     (if (math-nearly-equal-float g2 guess)
X	 g2
X       (math-sqrt-float-iter a g2)))
)
X
;;; True if A and B differ only in the last digit of precision.  [P F F]
(defun math-nearly-equal-float (a b)
X  (let ((ediff (- (nth 2 a) (nth 2 b))))
X    (cond ((= ediff 0)   ;; Expanded out for speed
X	   (setq ediff (math-add (Math-integer-neg (nth 1 a)) (nth 1 b)))
X	   (or (eq ediff 0)
X	       (and (not (consp ediff))
X		    (< ediff 10)
X		    (> ediff -10)
X		    (= (math-numdigs (nth 1 a)) calc-internal-prec))))
X	  ((= ediff 1)
X	   (setq ediff (math-add (Math-integer-neg (nth 1 b))
X				 (math-scale-int (nth 1 a) 1)))
X	   (and (not (consp ediff))
X		(< ediff 10)
X		(> ediff -10)
X		(= (math-numdigs (nth 1 b)) calc-internal-prec)))
X	  ((= ediff -1)
X	   (setq ediff (math-add (Math-integer-neg (nth 1 a))
X				 (math-scale-int (nth 1 b) 1)))
X	   (and (not (consp ediff))
X		(< ediff 10)
X		(> ediff -10)
X		(= (math-numdigs (nth 1 a)) calc-internal-prec)))))
)
X
(defun math-nearly-equal (a b)   ;  [P N N] [Public]
X  (setq a (math-float a))
X  (setq b (math-float b))
X  (if (eq (car a) 'polar) (setq a (math-complex a)))
X  (if (eq (car b) 'polar) (setq b (math-complex b)))
X  (if (eq (car a) 'cplx)
X      (if (eq (car b) 'cplx)
X	  (and (or (math-nearly-equal-float (nth 1 a) (nth 1 b))
X		   (and (math-nearly-zerop-float (nth 1 a) (nth 2 a))
X			(math-nearly-zerop-float (nth 1 b) (nth 2 b))))
X	       (or (math-nearly-equal-float (nth 2 a) (nth 2 b))
X		   (and (math-nearly-zerop-float (nth 2 a) (nth 1 a))
X			(math-nearly-zerop-float (nth 2 b) (nth 1 b)))))
X	(and (math-nearly-equal-float (nth 1 a) b)
X	     (math-nearly-zerop-float (nth 2 a) b)))
X      (if (eq (car b) 'cplx)
X	  (and (math-nearly-equal-float a (nth 1 b))
X	       (math-nearly-zerop-float a (nth 2 b)))
X	(math-nearly-equal-float a b)))
)
X
;;; True if A is nearly zero compared to B.  [P F F]
(defun math-nearly-zerop-float (a b)
X  (or (eq (nth 1 a) 0)
X      (<= (+ (math-numdigs (nth 1 a)) (nth 2 a))
X	  (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec))))
)
X
(defun math-nearly-zerop (a b)   ; [P N R] [Public]
X  (setq a (math-float a))
X  (setq b (math-float b))
X  (if (eq (car a) 'cplx)
X      (and (math-nearly-zerop-float (nth 1 a) b)
X	   (math-nearly-zerop-float (nth 2 a) b))
X    (if (eq (car a) 'polar)
X	(math-nearly-zerop-float (nth 1 a) b)
X      (math-nearly-zerop-float a b)))
)
X
;;; This implementation could be improved, accuracy-wise.
(defun math-hypot (a b)
X  (cond ((Math-zerop a) (math-abs b))
X	((Math-zerop b) (math-abs a))
X	((not (Math-scalarp a))
X	 (if (math-infinitep a)
X	     (if (math-infinitep b)
X		 (if (equal a b)
X		     a
X		   '(var nan var-nan))
X	       a)
X	   (calc-record-why 'scalarp a)
X	   (list 'calcFunc-hypot a b)))
X	((not (Math-scalarp b))
X	 (if (math-infinitep b)
X	     b
X	   (calc-record-why 'scalarp b)
X	   (list 'calcFunc-hypot a b)))
X	((and (Math-numberp a) (Math-numberp b))
X	 (math-with-extra-prec 1
X	   (math-sqrt (math-add (calcFunc-abssqr a) (calcFunc-abssqr b)))))
X	((eq (car-safe a) 'hms)
X	 (if (eq (car-safe b) 'hms)   ; this helps sdev's of hms forms
X	     (math-to-hms (math-hypot (math-from-hms a 'deg)
X				      (math-from-hms b 'deg)))
X	   (math-to-hms (math-hypot (math-from-hms a 'deg) b))))
X	((eq (car-safe b) 'hms)
X	 (math-to-hms (math-hypot a (math-from-hms b 'deg))))
X	(t nil))
)
(fset 'calcFunc-hypot (symbol-function 'math-hypot))
X
(defun calcFunc-sqr (x)
X  (math-pow x 2)
)
X
X
X
(defun math-nth-root (a n)
X  (cond ((= n 2) (math-sqrt a))
X	((Math-zerop a) a)
X	((Math-negp a) nil)
X	((Math-integerp a)
X	 (let ((root (math-nth-root-integer a n)))
X	   (if (car root)
X	       (cdr root)
X	     (and (not calc-symbolic-mode)
X		  (math-nth-root-float (math-float a) n
X				       (math-float (cdr root)))))))
X	((eq (car-safe a) 'frac)
X	 (let* ((num-root (math-nth-root-integer (nth 1 a) n))
X		(den-root (math-nth-root-integer (nth 2 a) n)))
X	   (if (and (car num-root) (car den-root))
X	       (list 'frac (cdr num-root) (cdr den-root))
X	     (and (not calc-symbolic-mode)
X		  (math-nth-root-float
X		   (math-float a) n
X		   (math-div-float (math-float (cdr num-root))
X				   (math-float (cdr den-root))))))))
X	((eq (car-safe a) 'float)
X	 (and (not calc-symbolic-mode)
X	      (math-nth-root-float a n)))
X	((eq (car-safe a) 'polar)
X	 (let ((root (math-nth-root (nth 1 a) n)))
X	   (and root (list 'polar root (math-div (nth 2 a) n)))))
X	(t nil))
)
X
(defun math-nth-root-float (a n &optional guess)
X  (math-inexact-result)
X  (math-with-extra-prec 1
X    (let ((nf (math-float n))
X	  (nfm1 (math-float (1- n))))
X      (math-nth-root-float-iter a (or guess
X				      (math-make-float
X				       1 (/ (+ (math-numdigs (nth 1 a))
X					       (nth 2 a)
X					       (/ n 2))
X					    n))))))
)
X
(defun math-nth-root-float-iter (a guess)   ; uses "n", "nf", "nfm1"
X  (math-working "root" guess)
X  (let ((g2 (math-div-float (math-add-float (math-mul nfm1 guess)
X					    (math-div-float
X					     a (math-ipow guess (1- n))))
X			    nf)))
X    (if (math-nearly-equal-float g2 guess)
X	g2
X      (math-nth-root-float-iter a g2)))
)
X
(defun math-nth-root-integer (a n &optional guess)   ; [I I S]
X  (math-nth-root-int-iter a (or guess
X				(math-scale-int 1 (/ (+ (math-numdigs a)
X							(1- n))
X						     n))))
)
X
(defun math-nth-root-int-iter (a guess)   ; uses "n"
X  (math-working "root" guess)
X  (let* ((q (math-idivmod a (math-ipow guess (1- n))))
X	 (s (math-add (car q) (math-mul (1- n) guess)))
X	 (g2 (math-idivmod s n)))
X    (if (Math-natnum-lessp (car g2) guess)
X	(math-nth-root-int-iter a (car g2))
X      (cons (and (equal (car g2) guess)
X		 (eq (cdr q) 0)
X		 (eq (cdr g2) 0))
X	    guess)))
)
X
(defun calcFunc-nroot (x n)
X  (calcFunc-pow x (if (integerp n)
X		      (math-make-frac 1 n)
X		    (math-div 1 n)))
)
X
X
X
X
;;;; Transcendental functions.
X
;;; All of these functions are defined on the complex plane.
;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
X
;;; Most functions increase calc-internal-prec by 2 digits, then round
;;; down afterward.  "-raw" functions use the current precision, require
;;; their arguments to be in float (or complex float) format, and always
;;; work in radians (where applicable).
X
(defun math-to-radians (a)   ; [N N]
X  (cond ((eq (car-safe a) 'hms)
X	 (math-from-hms a 'rad))
X	((memq calc-angle-mode '(deg hms))
X	 (math-mul a (math-pi-over-180)))
X	(t a))
)
X
(defun math-from-radians (a)   ; [N N]
X  (cond ((eq calc-angle-mode 'deg)
X	 (if (math-constp a)
X	     (math-div a (math-pi-over-180))
X	   (list 'calcFunc-deg a)))
X	((eq calc-angle-mode 'hms)
X	 (math-to-hms a 'rad))
X	(t a))
)
X
(defun math-to-radians-2 (a)   ; [N N]
X  (cond ((eq (car-safe a) 'hms)
X	 (math-from-hms a 'rad))
X	((memq calc-angle-mode '(deg hms))
X	 (if calc-symbolic-mode
X	     (math-div (math-mul a '(var pi var-pi)) 180)
X	   (math-mul a (math-pi-over-180))))
X	(t a))
)
X
(defun math-from-radians-2 (a)   ; [N N]
X  (cond ((memq calc-angle-mode '(deg hms))
X	 (if calc-symbolic-mode
X	     (math-div (math-mul 180 a) '(var pi var-pi))
X	   (math-div a (math-pi-over-180))))
X	(t a))
)
X
X
X
;;; Sine, cosine, and tangent.
X
(defun calcFunc-sin (x)   ; [N N] [Public]
X  (cond ((and (integerp x)
X	      (if (eq calc-angle-mode 'deg)
X		  (= (% x 90) 0)
X		(= x 0)))
X	 (aref [0 1 0 -1] (math-mod (/ x 90) 4)))
X	((Math-scalarp x)
X	 (math-with-extra-prec 2
X	   (math-sin-raw (math-to-radians (math-float x)))))
X	((eq (car x) 'sdev)
X	 (if (math-constp x)
X	     (math-with-extra-prec 2
X	       (let* ((xx (math-to-radians (math-float (nth 1 x))))
X		      (xs (math-to-radians (math-float (nth 2 x))))
X		      (sc (math-sin-cos-raw xx)))
X		 (math-make-sdev (car sc) (math-mul xs (cdr sc)))))
X	   (math-make-sdev (calcFunc-sin (nth 1 x))
X			   (math-mul (nth 2 x) (calcFunc-cos (nth 1 x))))))
X	((and (eq (car x) 'intv) (math-intv-constp x))
X	 (calcFunc-cos (math-sub x (math-quarter-circle nil))))
X	((equal x '(var nan var-nan))
X	 x)
X	(t (calc-record-why 'scalarp x)
X	   (list 'calcFunc-sin x)))
)
X
(defun calcFunc-cos (x)   ; [N N] [Public]
X  (cond ((and (integerp x)
X	      (if (eq calc-angle-mode 'deg)
X		  (= (% x 90) 0)
X		(= x 0)))
X	 (aref [1 0 -1 0] (math-mod (/ x 90) 4)))
X	((Math-scalarp x)
X	 (math-with-extra-prec 2
X	   (math-cos-raw (math-to-radians (math-float x)))))
X	((eq (car x) 'sdev)
X	 (if (math-constp x)
X	     (math-with-extra-prec 2
X	       (let* ((xx (math-to-radians (math-float (nth 1 x))))
X		      (xs (math-to-radians (math-float (nth 2 x))))
X		      (sc (math-sin-cos-raw xx)))
X		 (math-make-sdev (cdr sc) (math-mul xs (car sc)))))
X	   (math-make-sdev (calcFunc-cos (nth 1 x))
X			   (math-mul (nth 2 x) (calcFunc-sin (nth 1 x))))))
X	((and (eq (car x) 'intv) (math-intv-constp x))
X	 (math-with-extra-prec 2
X	   (let* ((xx (math-to-radians (math-float x)))
X		  (na (math-floor (math-div (nth 2 xx) (math-pi))))
X		  (nb (math-floor (math-div (nth 3 xx) (math-pi))))
X		  (span (math-sub nb na)))
X	     (if (memq span '(0 1))
X		 (let ((int (math-sort-intv (nth 1 x)
X					    (math-cos-raw (nth 2 xx))
X					    (math-cos-raw (nth 3 xx)))))
X		   (if (eq span 1)
X		       (if (math-evenp na)
X			   (math-make-intv (logior (nth 1 x) 2)
X					   -1
X					   (nth 3 int))
X			 (math-make-intv (logior (nth 1 x) 1)
X					 (nth 2 int)
X					 1))
X		     int))
X	       (list 'intv 3 -1 1)))))
X	((equal x '(var nan var-nan))
X	 x)
X	(t (calc-record-why 'scalarp x)
X	   (list 'calcFunc-cos x)))
)
X
(defun calcFunc-sincos (x)   ; [V N] [Public]
X  (if (Math-scalarp x)
X      (math-with-extra-prec 2
X	(let ((sc (math-sin-cos-raw (math-to-radians (math-float x)))))
X	  (list 'vec (cdr sc) (car sc))))    ; the vector [cos, sin]
X    (list 'vec (calcFunc-sin x) (calcFunc-cos x)))
)
X
(defun calcFunc-tan (x)   ; [N N] [Public]
X  (cond ((and (integerp x)
X	      (if (eq calc-angle-mode 'deg)
X		  (= (% x 180) 0)
X		(= x 0)))
X	 0)
X	((Math-scalarp x)
X	 (math-with-extra-prec 2
X	   (math-tan-raw (math-to-radians (math-float x)))))
X	((eq (car x) 'sdev)
X	 (if (math-constp x)
X	     (math-with-extra-prec 2
X	       (let* ((xx (math-to-radians (math-float (nth 1 x))))
X		      (xs (math-to-radians (math-float (nth 2 x))))
X		      (sc (math-sin-cos-raw xx)))
X		 (if (and (math-zerop (cdr sc)) (not calc-infinite-mode))
X		     (progn
X		       (calc-record-why "*Division by zero")
X		       (list 'calcFunc-tan x))
X		   (math-make-sdev (math-div-float (car sc) (cdr sc))
X				   (math-div-float xs (math-sqr (cdr sc)))))))
X	   (math-make-sdev (calcFunc-tan (nth 1 x))
X			   (math-div (nth 2 x)
X				     (math-sqr (calcFunc-cos (nth 1 x)))))))
X	((and (eq (car x) 'intv) (math-intv-constp x))
X	 (or (math-with-extra-prec 2
X	       (let* ((xx (math-to-radians (math-float x)))
X		      (na (math-floor (math-div (math-sub (nth 2 xx)
X							  (math-pi-over-2))
X						(math-pi))))
X		      (nb (math-floor (math-div (math-sub (nth 3 xx)
X							  (math-pi-over-2))
X						(math-pi)))))
X		 (and (equal na nb)
X		      (math-sort-intv (nth 1 x)
X				      (math-tan-raw (nth 2 xx))
X				      (math-tan-raw (nth 3 xx))))))
X	     '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
X	((equal x '(var nan var-nan))
X	 x)
X	(t (calc-record-why 'scalarp x)
X	   (list 'calcFunc-tan x)))
)
X
(defun math-sin-raw (x)   ; [N N]
X  (cond ((eq (car x) 'cplx)
X	 (let* ((expx (math-exp-raw (nth 2 x)))
X		(expmx (math-div-float '(float 1 0) expx))
X		(sc (math-sin-cos-raw (nth 1 x))))
X	   (list 'cplx
X		 (math-mul-float (car sc)
X				 (math-mul-float (math-add-float expx expmx)
X						 '(float 5 -1)))
X		 (math-mul-float (cdr sc)
X				 (math-mul-float (math-sub-float expx expmx)
X						 '(float 5 -1))))))
X	((eq (car x) 'polar)
X	 (math-polar (math-sin-raw (math-complex x))))
X	((Math-integer-negp (nth 1 x))
X	 (math-neg-float (math-sin-raw (math-neg-float x))))
X	((math-lessp-float '(float 7 0) x)  ; avoid inf loops due to roundoff
X	 (math-sin-raw (math-mod x (math-two-pi))))
X	(t (math-sin-raw-2 x x)))
)
X
(defun math-cos-raw (x)   ; [N N]
X  (if (eq (car-safe x) 'polar)
X      (math-polar (math-cos-raw (math-complex x)))
X    (math-sin-raw (math-sub (math-pi-over-2) x)))
)
X
;;; This could use a smarter method:  Reduce x as in math-sin-raw, then
;;;   compute either sin(x) or cos(x), whichever is smaller, and compute
;;;   the other using the identity sin(x)^2 + cos(x)^2 = 1.
(defun math-sin-cos-raw (x)   ; [F.F F]  (result is (sin x . cos x))
X  (cons (math-sin-raw x) (math-cos-raw x))
)
X
(defun math-tan-raw (x)   ; [N N]
X  (cond ((eq (car x) 'cplx)
X	 (let* ((x (math-mul x '(float 2 0)))
X		(expx (math-exp-raw (nth 2 x)))
X		(expmx (math-div-float '(float 1 0) expx))
X		(sc (math-sin-cos-raw (nth 1 x)))
X		(d (math-add-float (cdr sc)
X				   (math-mul-float (math-add-float expx expmx)
X						   '(float 5 -1)))))
X	   (and (not (eq (nth 1 d) 0))
X		(list 'cplx
X		      (math-div-float (car sc) d)
X		      (math-div-float (math-mul-float (math-sub-float expx
X								      expmx)
X						      '(float 5 -1)) d)))))
X	((eq (car x) 'polar)
X	 (math-polar (math-tan-raw (math-complex x))))
X	(t
X	 (let ((sc (math-sin-cos-raw x)))
X	   (if (eq (nth 1 (cdr sc)) 0)
X	       (math-div (car sc) 0)
X	     (math-div-float (car sc) (cdr sc))))))
)
X
(defun math-sin-raw-2 (x orgx)   ; This avoids poss of inf recursion.  [F F]
X  (let ((xmpo2 (math-sub-float (math-pi-over-2) x)))
X    (cond ((Math-integer-negp (nth 1 xmpo2))
X	   (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi))
X					   orgx)))
X	  ((math-lessp-float (math-pi-over-4) x)
X	   (math-cos-raw-2 xmpo2 orgx))
X	  ((math-lessp-float x (math-neg (math-pi-over-4)))
X	   (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx)))
X	  ((math-nearly-zerop-float x orgx) '(float 0 0))
X	  (calc-symbolic-mode (signal 'inexact-result nil))
X	  (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x))))))
)
X
(defun math-cos-raw-2 (x orgx)   ; [F F]
X  (cond ((math-nearly-zerop-float x orgx) '(float 1 0))
X	(calc-symbolic-mode (signal 'inexact-result nil))
X	(t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
X	     (math-sin-series
X	      (math-add-float '(float 1 0)
X			      (math-mul-float xnegsqr '(float 5 -1)))
X	      24 5 xnegsqr xnegsqr))))
)
X
(defun math-sin-series (sum nfac n x xnegsqr)
X  (math-working "sin" sum)
X  (let* ((nextx (math-mul-float x xnegsqr))
X	 (nextsum (math-add-float sum (math-div-float nextx
X						      (math-float nfac)))))
X    (if (math-nearly-equal-float sum nextsum)
X	sum
X      (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
X		       (+ n 2) nextx xnegsqr)))
)
X
X
;;; Inverse sine, cosine, tangent.
X
(defun calcFunc-arcsin (x)   ; [N N] [Public]
X  (cond ((eq x 0) 0)
X	((and (eq x 1) (eq calc-angle-mode 'deg)) 90)
X	((and (eq x -1) (eq calc-angle-mode 'deg)) -90)
X	(calc-symbolic-mode (signal 'inexact-result nil))
X	((Math-numberp x)
X	 (math-with-extra-prec 2
X	   (math-from-radians (math-arcsin-raw (math-float x)))))
X	((eq (car x) 'sdev)
X	 (math-make-sdev (calcFunc-arcsin (nth 1 x))
X			 (math-from-radians
X			  (math-div (nth 2 x)
X				    (math-sqrt
X				     (math-sub 1 (math-sqr (nth 1 x))))))))
X	((eq (car x) 'intv)
X	 (math-sort-intv (nth 1 x)
X			 (calcFunc-arcsin (nth 2 x))
X			 (calcFunc-arcsin (nth 3 x))))
X	((equal x '(var nan var-nan))
X	 x)
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-arcsin x)))
)
X
(defun calcFunc-arccos (x)   ; [N N] [Public]
X  (cond ((eq x 1) 0)
X	((and (eq x 0) (eq calc-angle-mode 'deg)) 90)
X	((and (eq x -1) (eq calc-angle-mode 'deg)) 180)
X	(calc-symbolic-mode (signal 'inexact-result nil))
X	((Math-numberp x)
X	 (math-with-extra-prec 2
X	   (math-from-radians (math-arccos-raw (math-float x)))))
X	((eq (car x) 'sdev)
X	 (math-make-sdev (calcFunc-arccos (nth 1 x))
X			 (math-from-radians
X			  (math-div (nth 2 x)
X				    (math-sqrt
X				     (math-sub 1 (math-sqr (nth 1 x))))))))
X	((eq (car x) 'intv)
X	 (math-sort-intv (nth 1 x)
X			 (calcFunc-arccos (nth 2 x))
X			 (calcFunc-arccos (nth 3 x))))
X	((equal x '(var nan var-nan))
X	 x)
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-arccos x)))
)
X
(defun calcFunc-arctan (x)   ; [N N] [Public]
X  (cond ((eq x 0) 0)
X	((and (eq x 1) (eq calc-angle-mode 'deg)) 45)
X	((and (eq x -1) (eq calc-angle-mode 'deg)) -45)
X	((Math-numberp x)
X	 (math-with-extra-prec 2
X	   (math-from-radians (math-arctan-raw (math-float x)))))
X	((eq (car x) 'sdev)
X	 (math-make-sdev (calcFunc-arctan (nth 1 x))
X			 (math-from-radians
X			  (math-div (nth 2 x)
X				    (math-add 1 (math-sqr (nth 1 x)))))))
X	((eq (car x) 'intv)
X	 (math-sort-intv (nth 1 x)
X			 (calcFunc-arctan (nth 2 x))
X			 (calcFunc-arctan (nth 3 x))))
X	((equal x '(var inf var-inf))
X	 (math-quarter-circle t))
X	((equal x '(neg (var inf var-inf)))
X	 (math-neg (math-quarter-circle t)))
X	((equal x '(var nan var-nan))
X	 x)
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-arctan x)))
)
X
(defun math-arcsin-raw (x)   ; [N N]
X  (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x)))))
X    (if (or (memq (car x) '(cplx polar))
X	    (memq (car a) '(cplx polar)))
X	(math-with-extra-prec 2   ; use extra precision for difficult case
X	  (math-mul '(cplx 0 -1)
X		    (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a))))
X      (math-arctan2-raw x a)))
)
X
(defun math-arccos-raw (x)   ; [N N]
X  (math-sub (math-pi-over-2) (math-arcsin-raw x))
)
X
(defun math-arctan-raw (x)   ; [N N]
X  (cond ((memq (car x) '(cplx polar))
X	 (math-with-extra-prec 2   ; extra-extra
X	   (math-div (math-sub
X		      (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x)))
X		      (math-ln-raw (math-add 1 (math-mul '(cplx 0 -1) x))))
X		     '(cplx 0 2))))
X	((Math-integer-negp (nth 1 x))
X	 (math-neg-float (math-arctan-raw (math-neg-float x))))
X	((math-zerop x) x)
X	(calc-symbolic-mode (signal 'inexact-result nil))
X	((math-equal-int x 1) (math-pi-over-4))
X	((math-equal-int x -1) (math-neg (math-pi-over-4)))
X	((math-lessp-float '(float 414214 -6) x)  ; if x > sqrt(2) - 1, reduce
X	 (if (math-lessp-float '(float 1 0) x)
X	     (math-sub-float (math-mul-float (math-pi) '(float 5 -1))
X			     (math-arctan-raw (math-div-float '(float 1 0) x)))
X	   (math-sub-float (math-mul-float (math-pi) '(float 25 -2))
X			   (math-arctan-raw (math-div-float
X					     (math-sub-float '(float 1 0) x)
X					     (math-add-float '(float 1 0)
X							     x))))))
X	(t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x)))))
)
X
(defun math-arctan-series (sum n x xnegsqr)
X  (math-working "arctan" sum)
X  (let* ((nextx (math-mul-float x xnegsqr))
X	 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
X    (if (math-nearly-equal-float sum nextsum)
X	sum
X      (math-arctan-series nextsum (+ n 2) nextx xnegsqr)))
)
X
(defun calcFunc-arctan2 (y x)   ; [F R R] [Public]
X  (if (Math-anglep y)
X      (if (Math-anglep x)
X	  (math-with-extra-prec 2
X	    (math-from-radians (math-arctan2-raw (math-float y)
X						 (math-float x))))
X	(calc-record-why 'anglep x)
X	(list 'calcFunc-arctan2 y x))
X    (if (and (or (math-infinitep x) (math-anglep x))
X	     (or (math-infinitep y) (math-anglep y)))
X	(progn
X	  (if (math-posp x)
X	      (setq x 1)
X	    (if (math-negp x)
X		(setq x -1)
X	      (or (math-zerop x)
X		  (setq x nil))))
X	  (if (math-posp y)
X	      (setq y 1)
X	    (if (math-negp y)
X		(setq y -1)
X	      (or (math-zerop y)
X		  (setq y nil))))
X	  (if (and y x)
X	      (calcFunc-arctan2 y x)
X	    '(var nan var-nan)))
X      (calc-record-why 'anglep y)
X      (list 'calcFunc-arctan2 y x)))
)
X
(defun math-arctan2-raw (y x)   ; [F R R]
X  (cond ((math-zerop y)
X	 (if (math-negp x) (math-pi)
X	   (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0)))
X	((math-zerop x)
X	 (if (math-posp y)
X	     (math-pi-over-2)
X	   (math-neg (math-pi-over-2))))
X	((math-posp x)
X	 (math-arctan-raw (math-div-float y x)))
X	((math-posp y)
X	 (math-add-float (math-arctan-raw (math-div-float y x))
X			 (math-pi)))
X	(t
X	 (math-sub-float (math-arctan-raw (math-div-float y x))
X			 (math-pi))))
)
X
(defun calcFunc-arcsincos (x)   ; [V N] [Public]
X  (if (and (Math-vectorp x)
X	   (= (length x) 3))
X      (calcFunc-arctan2 (nth 2 x) (nth 1 x))
X    (math-reject-arg x "*Two-element vector expected"))
)
X
X
X
;;; Exponential function.
X
(defun calcFunc-exp (x)   ; [N N] [Public]
X  (cond ((eq x 0) 1)
X	((and (memq x '(1 -1)) calc-symbolic-mode)
X	 (if (eq x 1) '(var e var-e) (math-div 1 '(var e var-e))))
X	((Math-numberp x)
X	 (math-with-extra-prec 2 (math-exp-raw (math-float x))))
X	((eq (car-safe x) 'sdev)
X	 (let ((ex (calcFunc-exp (nth 1 x))))
X	   (math-make-sdev ex (math-mul (nth 2 x) ex))))
X	((eq (car-safe x) 'intv)
X	 (math-make-intv (nth 1 x) (calcFunc-exp (nth 2 x))
X			 (calcFunc-exp (nth 3 x))))
X	((equal x '(var inf var-inf))
X	 x)
X	((equal x '(neg (var inf var-inf)))
X	 0)
X	((equal x '(var nan var-nan))
X	 x)
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-exp x)))
)
X
(defun calcFunc-expm1 (x)   ; [N N] [Public]
X  (cond ((eq x 0) 0)
X	((math-zerop x) '(float 0 0))
X	(calc-symbolic-mode (signal 'inexact-result nil))
X	((Math-numberp x)
X	 (math-with-extra-prec 2
X	   (let ((x (math-float x)))
X	     (if (and (eq (car x) 'float)
X		      (math-lessp-float x '(float 1 0))
X		      (math-lessp-float '(float -1 0) x))
X		 (math-exp-minus-1-raw x)
X	       (math-add (math-exp-raw x) -1)))))
X	((eq (car-safe x) 'sdev)
X	 (if (math-constp x)
X	     (let ((ex (calcFunc-expm1 (nth 1 x))))
X	       (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1))))
X	   (math-make-sdev (calcFunc-expm1 (nth 1 x))
X			   (math-mul (nth 2 x) (calcFunc-exp (nth 1 x))))))
X	((eq (car-safe x) 'intv)
X	 (math-make-intv (nth 1 x)
X			 (calcFunc-expm1 (nth 2 x))
X			 (calcFunc-expm1 (nth 3 x))))
X	((equal x '(var inf var-inf))
X	 x)
X	((equal x '(neg (var inf var-inf)))
X	 -1)
X	((equal x '(var nan var-nan))
X	 x)
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-expm1 x)))
)
X
(defun calcFunc-exp10 (x)   ; [N N] [Public]
X  (if (eq x 0)
X      1
X    (math-pow '(float 1 1) x))
)
X
(defun math-exp-raw (x)   ; [N N]
X  (cond ((math-zerop x) '(float 1 0))
X	(calc-symbolic-mode (signal 'inexact-result nil))
X	((eq (car x) 'cplx)
X	 (let ((expx (math-exp-raw (nth 1 x)))
X	       (sc (math-sin-cos-raw (nth 2 x))))
X	   (list 'cplx
X		 (math-mul-float expx (cdr sc))
X		 (math-mul-float expx (car sc)))))
X	((eq (car x) 'polar)
X	 (let ((xc (math-complex x)))
X	   (list 'polar
X		 (math-exp-raw (nth 1 xc))
X		 (math-from-radians (nth 2 xc)))))
X	((or (math-lessp-float '(float 5 -1) x)
X	     (math-lessp-float x '(float -5 -1)))
X	 (if (math-lessp-float '(float 921035 1) x)
X	     (math-overflow)
X	   (if (math-lessp-float x '(float -921035 1))
X	       (math-underflow)))
X	 (let* ((two-x (math-mul-float x '(float 2 0)))
X		(hint (math-scale-int (nth 1 two-x) (nth 2 two-x)))
X		(hfrac (math-sub-float x (math-mul-float (math-float hint)
X							 '(float 5 -1)))))
X	   (math-mul-float (math-ipow (math-sqrt-e) hint)
X			   (math-add-float '(float 1 0)
X					   (math-exp-minus-1-raw hfrac)))))
X	(t (math-add-float '(float 1 0) (math-exp-minus-1-raw x))))
)
X
(defun math-exp-minus-1-raw (x)   ; [F F]
X  (math-exp-series x 2 3 x x)
)
X
(defun math-exp-series (sum nfac n xpow x)
X  (math-working "exp" sum)
X  (let* ((nextx (math-mul-float xpow x))
X	 (nextsum (math-add-float sum (math-div-float nextx
X						      (math-float nfac)))))
X    (if (math-nearly-equal-float sum nextsum)
X	sum
X      (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x)))
)
X
X
X
;;; Logarithms.
X
(defun calcFunc-ln (x)   ; [N N] [Public]
X  (cond ((math-zerop x)
X	 (if calc-infinite-mode
X	     '(neg (var inf var-inf))
X	   (math-reject-arg x "*Logarithm of zero")))
X	((eq x 1) 0)
X	((Math-numberp x)
X	 (math-with-extra-prec 2 (math-ln-raw (math-float x))))
X	((eq (car-safe x) 'sdev)
X	 (math-make-sdev (calcFunc-ln (nth 1 x))
X			 (math-div (nth 2 x) (nth 1 x))))
X	((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
X					  (Math-zerop (nth 2 x))
X					  (not (math-intv-constp x))))
X	 (let ((calc-infinite-mode t))
X	   (math-make-intv (nth 1 x) (calcFunc-ln (nth 2 x))
X			   (calcFunc-ln (nth 3 x)))))
X	((equal x '(var e var-e))
X	 1)
X	((and (eq (car-safe x) '^)
X	      (equal (nth 1 x) '(var e var-e))
X	      (math-known-realp (nth 2 x)))
X	 (nth 2 x))
X	((math-infinitep x)
X	 (if (equal x '(var nan var-nan))
X	     x
X	   '(var inf var-inf)))
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-ln x)))
)
X
(defun calcFunc-log10 (x)   ; [N N] [Public]
X  (cond ((math-equal-int x 1)
X	 (if (math-floatp x) '(float 0 0) 0))
X	((and (Math-integerp x)
X	      (math-posp x)
X	      (let ((res (math-integer-log x 10)))
X		(and (car res)
X		     (setq x (cdr res)))))
X	 x)
X	((and (eq (car-safe x) 'frac)
X	      (eq (nth 1 x) 1)
X	      (let ((res (math-integer-log (nth 2 x) 10)))
X		(and (car res)
X		     (setq x (- (cdr res))))))
X	 x)
X	((math-zerop x)
X	 (if calc-infinite-mode
X	     '(neg (var inf var-inf))
X	   (math-reject-arg x "*Logarithm of zero")))
X	(calc-symbolic-mode (signal 'inexact-result nil))
X	((Math-numberp x)
X	 (math-with-extra-prec 2
X	   (let ((xf (math-float x)))
X	     (if (eq (nth 1 xf) 0)
X		 (math-reject-arg x "*Logarithm of zero"))
X	     (if (Math-integer-posp (nth 1 xf))
X		 (if (eq (nth 1 xf) 1)    ; log10(1*10^n) = n
X		     (math-float (nth 2 xf))
X		   (let ((xdigs (1- (math-numdigs (nth 1 xf)))))
X		     (math-add-float
X		      (math-div-float (math-ln-raw-2
X				       (list 'float (nth 1 xf) (- xdigs)))
X				      (math-ln-10))
X		      (math-float (+ (nth 2 xf) xdigs)))))
X	       (math-div (calcFunc-ln xf) (math-ln-10))))))
X	((eq (car-safe x) 'sdev)
X	 (math-make-sdev (calcFunc-log10 (nth 1 x))
X			 (math-div (nth 2 x)
X				   (math-mul (nth 1 x) (math-ln-10)))))
X	((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
X					  (not (math-intv-constp x))))
X	 (math-make-intv (nth 1 x)
X			 (calcFunc-log10 (nth 2 x))
X			 (calcFunc-log10 (nth 3 x))))
X	((math-infinitep x)
X	 (if (equal x '(var nan var-nan))
X	     x
X	   '(var inf var-inf)))
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-log10 x)))
)
X
(defun calcFunc-log (x &optional b)   ; [N N N] [Public]
X  (cond ((or (null b) (equal b '(var e var-e)))
X	 (math-normalize (list 'calcFunc-ln x)))
X	((or (eq b 10) (equal b '(float 1 1)))
X	 (math-normalize (list 'calcFunc-log10 x)))
X	((math-zerop x)
X	 (if calc-infinite-mode
X	     (math-div (calcFunc-ln x) (calcFunc-ln b))
X	   (math-reject-arg x "*Logarithm of zero")))
X	((math-zerop b)
X	 (if calc-infinite-mode
X	     (math-div (calcFunc-ln x) (calcFunc-ln b))
X	   (math-reject-arg b "*Logarithm of zero")))
X	((math-equal-int b 1)
X	 (if calc-infinite-mode
X	     (math-div (calcFunc-ln x) 0)
X	   (math-reject-arg b "*Logarithm base one")))
X	((math-equal-int x 1)
X	 (if (or (math-floatp a) (math-floatp b)) '(float 0 0) 0))
X	((and (Math-ratp x) (Math-ratp b)
X	      (math-posp x) (math-posp b)
X	      (let* ((sign 1) (inv nil)
X		     (xx (if (Math-lessp 1 x)
X			     x
X			   (setq sign -1)
X			   (math-div 1 x)))
X		     (bb (if (Math-lessp 1 b)
X			     b
X			   (setq sign (- sign))
X			   (math-div 1 b)))
X		     (res (if (Math-lessp xx bb)
X			      (setq inv (math-integer-log bb xx))
X			    (math-integer-log xx bb))))
X		(and (car res)
X		     (setq x (if inv
X				 (math-div 1 (* sign (cdr res)))
X			       (* sign (cdr res)))))))
X	 x)
X	(calc-symbolic-mode (signal 'inexact-result nil))
X	((and (Math-numberp x) (Math-numberp b))
X	 (math-with-extra-prec 2
X	   (math-div (math-ln-raw (math-float x))
X		     (math-log-base-raw b))))
X	((and (eq (car-safe x) 'sdev)
X	      (Math-numberp b))
X	 (math-make-sdev (calcFunc-log (nth 1 x) b)
X			 (math-div (nth 2 x)
X				   (math-mul (nth 1 x)
X					     (math-log-base-raw b)))))
X	((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
X					  (not (math-intv-constp x)))
X	      (math-realp b))
X	 (math-make-intv (nth 1 x)
X			 (calcFunc-log (nth 2 x) b)
X			 (calcFunc-log (nth 3 x) b)))
X	((or (eq (car-safe x) 'intv) (eq (car-safe b) 'intv))
X	 (math-div (calcFunc-ln x) (calcFunc-ln b)))
X	((or (math-infinitep x)
X	     (math-infinitep b))
X	 (math-div (calcFunc-ln x) (calcFunc-ln b)))
X	(t (if (Math-numberp b)
X	       (calc-record-why 'numberp x)
X	     (calc-record-why 'numberp b))
X	   (list 'calcFunc-log x b)))
)
X
(defun calcFunc-alog (x &optional b)
X  (cond ((or (null b) (equal b '(var e var-e)))
X	 (math-normalize (list 'calcFunc-exp x)))
X	(t (math-pow b x)))
)
X
(defun calcFunc-ilog (x b)
X  (if (and (math-natnump x) (not (eq x 0))
X	   (math-natnump b) (not (eq b 0)))
X      (if (eq b 1)
X	  (math-reject-arg x "*Logarithm base one")
X	(if (Math-natnum-lessp x b)
X	    0
X	  (cdr (math-integer-log x b))))
X    (math-floor (calcFunc-log x b)))
)
X
(defun math-integer-log (x b)
X  (let ((pows (list b))
X	(pow (math-sqr b))
X	next
X	sum n)
X    (while (not (Math-lessp x pow))
X      (setq pows (cons pow pows)
X	    pow (math-sqr pow)))
X    (setq n (lsh 1 (1- (length pows)))
X	  sum n
X	  pow (car pows))
X    (while (and (setq pows (cdr pows))
X		(Math-lessp pow x))
X      (setq n (/ n 2)
X	    next (math-mul pow (car pows)))
X      (or (Math-lessp x next)
X	  (setq pow next
X		sum (+ sum n))))
X    (cons (equal pow x) sum))
)
X
X
(defun math-log-base-raw (b)   ; [N N]
X  (if (not (and (equal (car math-log-base-cache) b)
X		(eq (nth 1 math-log-base-cache) calc-internal-prec)))
X      (setq math-log-base-cache (list b calc-internal-prec
X				      (math-ln-raw (math-float b)))))
X  (nth 2 math-log-base-cache)
)
(setq math-log-base-cache nil)
X
(defun calcFunc-lnp1 (x)   ; [N N] [Public]
X  (cond ((Math-equal-int x -1)
X	 (if calc-infinite-mode
X	     '(neg (var inf var-inf))
X	   (math-reject-arg x "*Logarithm of zero")))
X	((eq x 0) 0)
X	((math-zerop x) '(float 0 0))
X	(calc-symbolic-mode (signal 'inexact-result nil))
X	((Math-numberp x)
X	 (math-with-extra-prec 2
X	   (let ((x (math-float x)))
X	     (if (and (eq (car x) 'float)
X		      (math-lessp-float x '(float 5 -1))
X		      (math-lessp-float '(float -5 -1) x))
X		 (math-ln-plus-1-raw x)
X	       (math-ln-raw (math-add-float x '(float 1 0)))))))
X	((eq (car-safe x) 'sdev)
X	 (math-make-sdev (calcFunc-lnp1 (nth 1 x))
X			 (math-div (nth 2 x) (math-add (nth 1 x) 1))))
X	((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
X					  (not (math-intv-constp x))))
X	 (math-make-intv (nth 1 x)
X			 (calcFunc-lnp1 (nth 2 x))
X			 (calcFunc-lnp1 (nth 3 x))))
X	((math-infinitep x)
X	 (if (equal x '(var nan var-nan))
X	     x
X	   '(var inf var-inf)))
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-lnp1 x)))
)
X
(defun math-ln-raw (x)    ; [N N] --- must be float format!
X  (cond ((eq (car-safe x) 'cplx)
X	 (list 'cplx
X	       (math-mul-float (math-ln-raw
X				(math-add-float (math-sqr-float (nth 1 x))
X						(math-sqr-float (nth 2 x))))
X			       '(float 5 -1))
X	       (math-arctan2-raw (nth 2 x) (nth 1 x))))
X	((eq (car x) 'polar)
X	 (math-polar (list 'cplx
X			   (math-ln-raw (nth 1 x))
X			   (math-to-radians (nth 2 x)))))
X	((Math-equal-int x 1)
X	 '(float 0 0))
X	(calc-symbolic-mode (signal 'inexact-result nil))
X	((math-posp (nth 1 x))    ; positive and real
X	 (let ((xdigs (1- (math-numdigs (nth 1 x)))))
X	   (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs)))
X			   (math-mul-float (math-float (+ (nth 2 x) xdigs))
X					   (math-ln-10)))))
X	((math-zerop x)
X	 (math-reject-arg x "*Logarithm of zero"))
X	((eq calc-complex-mode 'polar)    ; negative and real
X	 (math-polar
X	  (list 'cplx   ; negative and real
X		(math-ln-raw (math-neg-float x))
X		(math-pi))))
X	(t (list 'cplx   ; negative and real
X		 (math-ln-raw (math-neg-float x))
X		 (math-pi))))
)
X
(defun math-ln-raw-2 (x)    ; [F F]
X  (cond ((math-lessp-float '(float 14 -1) x)
X	 (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1)))
X			 (math-ln-2)))
X	(t    ; now .7 < x <= 1.4
X	 (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
X					(math-add-float x '(float 1 0))))))
)
X
(defun math-ln-raw-3 (x)   ; [F F]
X  (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
X		  '(float 2 0))
)
X
;;; Compute ln((1+x)/(1-x))
(defun math-ln-raw-series (sum n x xsqr)
X  (math-working "log" sum)
X  (let* ((nextx (math-mul-float x xsqr))
X	 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
X    (if (math-nearly-equal-float sum nextsum)
X	sum
X      (math-ln-raw-series nextsum (+ n 2) nextx xsqr)))
)
X
(defun math-ln-plus-1-raw (x)
X  (math-lnp1-series x 2 x (math-neg x))
)
X
(defun math-lnp1-series (sum n xpow x)
X  (math-working "lnp1" sum)
X  (let* ((nextx (math-mul-float xpow x))
X	 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
X    (if (math-nearly-equal-float sum nextsum)
X	sum
X      (math-lnp1-series nextsum (1+ n) nextx x)))
)
X
(math-defcache math-ln-10 (float (bigpos 018 684 045 994 092 585 302 2) -21)
X  (math-ln-raw-2 '(float 1 1)))
X
(math-defcache math-ln-2 (float (bigpos 417 309 945 559 180 147 693) -21)
X  (math-ln-raw-3 (math-float '(frac 1 3))))
X
X
X
;;; Hyperbolic functions.
X
(defun calcFunc-sinh (x)   ; [N N] [Public]
X  (cond ((eq x 0) 0)
X	(math-expand-formulas
X	 (math-normalize
X	  (list '/ (list '- (list 'calcFunc-exp x)
X			 (list 'calcFunc-exp (list 'neg x))) 2)))
X	((Math-numberp x)
X	 (if calc-symbolic-mode (signal 'inexact-result nil))
X	 (math-with-extra-prec 2
X	   (let ((expx (math-exp-raw (math-float x))))
X	     (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1)))))
X	((eq (car-safe x) 'sdev)
X	 (math-make-sdev (calcFunc-sinh (nth 1 x))
X			 (math-mul (nth 2 x) (calcFunc-cosh (nth 1 x)))))
X	((eq (car x) 'intv)
X	 (math-sort-intv (nth 1 x)
X			 (calcFunc-sinh (nth 2 x))
X			 (calcFunc-sinh (nth 3 x))))
X	((or (equal x '(var inf var-inf))
X	     (equal x '(neg (var inf var-inf)))
X	     (equal x '(var nan var-nan)))
X	 x)
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-sinh x)))
)
(put 'calcFunc-sinh 'math-expandable t)
X
(defun calcFunc-cosh (x)   ; [N N] [Public]
X  (cond ((eq x 0) 1)
X	(math-expand-formulas
X	 (math-normalize
X	  (list '/ (list '+ (list 'calcFunc-exp x)
X			 (list 'calcFunc-exp (list 'neg x))) 2)))
X	((Math-numberp x)
X	 (if calc-symbolic-mode (signal 'inexact-result nil))
X	 (math-with-extra-prec 2
X	   (let ((expx (math-exp-raw (math-float x))))
X	     (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1)))))
X	((eq (car-safe x) 'sdev)
X	 (math-make-sdev (calcFunc-cosh (nth 1 x))
X			 (math-mul (nth 2 x)
X				   (calcFunc-sinh (nth 1 x)))))
X	((and (eq (car x) 'intv) (math-intv-constp x))
X	 (setq x (math-abs x))
X	 (math-sort-intv (nth 1 x)
X			 (calcFunc-cosh (nth 2 x))
X			 (calcFunc-cosh (nth 3 x))))
X	((or (equal x '(var inf var-inf))
X	     (equal x '(neg (var inf var-inf)))
X	     (equal x '(var nan var-nan)))
X	 (math-abs x))
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-cosh x)))
)
(put 'calcFunc-cosh 'math-expandable t)
X
(defun calcFunc-tanh (x)   ; [N N] [Public]
X  (cond ((eq x 0) 0)
X	(math-expand-formulas
X	 (math-normalize
X	  (let ((expx (list 'calcFunc-exp x))
X		(expmx (list 'calcFunc-exp (list 'neg x))))
X	    (math-normalize
X	     (list '/ (list '- expx expmx) (list '+ expx expmx))))))
X	((Math-numberp x)
X	 (if calc-symbolic-mode (signal 'inexact-result nil))
X	 (math-with-extra-prec 2
X	   (let* ((expx (calcFunc-exp (math-float x)))
X		  (expmx (math-div 1 expx)))
X	     (math-div (math-sub expx expmx)
X		       (math-add expx expmx)))))
X	((eq (car-safe x) 'sdev)
X	 (math-make-sdev (calcFunc-tanh (nth 1 x))
X			 (math-div (nth 2 x)
X				   (math-sqr (calcFunc-cosh (nth 1 x))))))
X	((eq (car x) 'intv)
X	 (math-sort-intv (nth 1 x)
X			 (calcFunc-tanh (nth 2 x))
X			 (calcFunc-tanh (nth 3 x))))
X	((equal x '(var inf var-inf))
X	 1)
X	((equal x '(neg (var inf var-inf)))
X	 -1)
X	((equal x '(var nan var-nan))
X	 x)
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-tanh x)))
)
(put 'calcFunc-tanh 'math-expandable t)
X
(defun calcFunc-arcsinh (x)   ; [N N] [Public]
X  (cond ((eq x 0) 0)
X	(math-expand-formulas
X	 (math-normalize
X	  (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
X					      (list '+ (list '^ x 2) 1))))))
X	((Math-numberp x)
X	 (if calc-symbolic-mode (signal 'inexact-result nil))
X	 (math-with-extra-prec 2
X	   (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x)
X							     '(float 1 0)))))))
X	((eq (car-safe x) 'sdev)
X	 (math-make-sdev (calcFunc-arcsinh (nth 1 x))
X			 (math-div (nth 2 x)
X				   (math-sqrt
X				    (math-add (math-sqr (nth 1 x)) 1)))))
X	((eq (car x) 'intv)
X	 (math-sort-intv (nth 1 x)
X			 (calcFunc-arcsinh (nth 2 x))
X			 (calcFunc-arcsinh (nth 3 x))))
X	((or (equal x '(var inf var-inf))
X	     (equal x '(neg (var inf var-inf)))
X	     (equal x '(var nan var-nan)))
X	 x)
X	(t (calc-record-why 'numberp x)
X	   (list 'calcFunc-arcsinh x)))
)
(put 'calcFunc-arcsinh 'math-expandable t)
X
(defun calcFunc-arccosh (x)   ; [N N] [Public]
X  (cond ((eq x 1) 0)
X	((and (eq x -1) calc-symbolic-mode)
X	 '(var pi var-pi))
X	((and (eq x 0) calc-symbolic-mode)
X	 (math-div (math-mul '(var pi var-pi) '(var i var-i)) 2))
X	(math-expand-formulas
X	 (math-normalize
X	  (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
X					      (list '- (list '^ x 2) 1))))))
X	((Math-numberp x)
X	 (if calc-symbolic-mode (signal 'inexact-result nil))
X	 (if (Math-equal-int x -1)
X	     (math-imaginary (math-pi))
X	   (math-with-extra-prec 2
X	     (if (or t    ; need to do this even in the real case!
X		     (memq (car-safe x) '(cplx polar)))
X		 (let ((xp1 (math-add 1 x)))  ; this gets the branch cuts right
X		   (math-ln-raw
X		    (math-add x (math-mul xp1
X					  (math-sqrt-raw
X					   (math-div (math-sub
X						      x
X						      '(float 1 0))
X						     xp1))))))
X	       (math-ln-raw
X		(math-add x (math-sqrt-raw (math-add (math-sqr x)
X						     '(float -1 0)))))))))
X	((eq (car-safe x) 'sdev)
X	 (math-make-sdev (calcFunc-arccosh (nth 1 x))
X			 (math-div (nth 2 x)
X				   (math-sqrt
X				    (math-add (math-sqr (nth 1 x)) -1)))))
X	((eq (car x) 'intv)
X	 (math-sort-intv (nth 1 x)
X			 (calcFunc-arccosh (nth 2 x))
X			 (calcFunc-arccosh (nth 3 x))))
X	((or (equal x '(var inf var-inf))
X	     (equal x '(neg (var inf var-inf)))
X	     (equal x '(var nan var-nan)))
SHAR_EOF
true || echo 'restore of calc-math.el failed'
fi
echo 'End of  part 21'
echo 'File calc-math.el is continued in part 22'
echo 22 > _shar_seq_.tmp
exit 0
exit 0 # Just in case...
-- 
Kent Landfield                   INTERNET: kent@sparky.IMD.Sterling.COM
Sterling Software, IMD           UUCP:     uunet!sparky!kent
Phone:    (402) 291-8300         FAX:      (402) 291-4362
Please send comp.sources.misc-related mail to kent@uunet.uu.net.