comparison lisp/calc/calc-math.el @ 41044:4549dec29728

(calcFunc-sqrt, calcFunc-hypot): Use `defalias' instead of `fset' and `symbol-function'. Style cleanup; don't put closing parens on their own line, add "foo.el ends here" to each file, and update copyright date.

41043:21a6b9fea031

;; Calculator for GNU Emacs, part II [calc-math.el]
;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc.
;; Written by Dave Gillespie, [e-mail].

;; This file is part of GNU Emacs.

;; GNU Emacs is distributed in the hope that it will be useful,

(defun calc-sqrt (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-unary-op "^2" 'calcFunc-sqr arg)
     (calc-unary-op "sqrt" 'calcFunc-sqrt arg)))
)

(defun calc-isqrt (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-unary-op "^2" 'calcFunc-sqr arg)
     (calc-unary-op "isqt" 'calcFunc-isqrt arg)))
)


(defun calc-hypot (arg)
  (interactive "P")
  (calc-slow-wrapper
   (calc-binary-op "hypt" 'calcFunc-hypot arg))
)

(defun calc-ln (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-exp arg)
)

(defun calc-log10 (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-ln arg)
)

(defun calc-log (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-binary-op "alog" 'calcFunc-alog arg)
     (calc-binary-op "log" 'calcFunc-log arg)))
)

(defun calc-ilog (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-binary-op "alog" 'calcFunc-alog arg)
     (calc-binary-op "ilog" 'calcFunc-ilog arg)))
)

(defun calc-lnp1 (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-expm1 arg)
)

(defun calc-exp (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (if (calc-is-inverse)
	   (calc-unary-op "lg10" 'calcFunc-log10 arg)
	 (calc-unary-op "10^" 'calcFunc-exp10 arg))
     (if (calc-is-inverse)
	 (calc-unary-op "ln" 'calcFunc-ln arg)
       (calc-unary-op "exp" 'calcFunc-exp arg))))
)

(defun calc-expm1 (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-unary-op "ln+1" 'calcFunc-lnp1 arg)
     (calc-unary-op "ex-1" 'calcFunc-expm1 arg)))
)

(defun calc-pi ()
  (interactive)
  (calc-slow-wrapper
   (if (calc-is-inverse)

	 (if calc-symbolic-mode
	     (calc-pop-push-record 0 "e" '(var e var-e))
	   (calc-pop-push-record 0 "e" (math-e)))
       (if calc-symbolic-mode
	   (calc-pop-push-record 0 "pi" '(var pi var-pi))
	 (calc-pop-push-record 0 "pi" (math-pi))))))
)

(defun calc-sin (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (if (calc-is-inverse)
	   (calc-unary-op "asnh" 'calcFunc-arcsinh arg)
	 (calc-unary-op "sinh" 'calcFunc-sinh arg))
     (if (calc-is-inverse)
	 (calc-unary-op "asin" 'calcFunc-arcsin arg)
       (calc-unary-op "sin" 'calcFunc-sin arg))))
)

(defun calc-arcsin (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-sin arg)
)

(defun calc-sinh (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-sin arg)
)

(defun calc-arcsinh (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-hyperbolic-func)
  (calc-sin arg)
)

(defun calc-cos (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (if (calc-is-inverse)
	   (calc-unary-op "acsh" 'calcFunc-arccosh arg)
	 (calc-unary-op "cosh" 'calcFunc-cosh arg))
     (if (calc-is-inverse)
	 (calc-unary-op "acos" 'calcFunc-arccos arg)
       (calc-unary-op "cos" 'calcFunc-cos arg))))
)

(defun calc-arccos (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-cos arg)
)

(defun calc-cosh (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-cos arg)
)

(defun calc-arccosh (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-hyperbolic-func)
  (calc-cos arg)
)

(defun calc-sincos ()
  (interactive)
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1)))
     (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1)))))
)

(defun calc-tan (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (if (calc-is-inverse)
	   (calc-unary-op "atnh" 'calcFunc-arctanh arg)
	 (calc-unary-op "tanh" 'calcFunc-tanh arg))
     (if (calc-is-inverse)
	 (calc-unary-op "atan" 'calcFunc-arctan arg)
       (calc-unary-op "tan" 'calcFunc-tan arg))))
)

(defun calc-arctan (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-tan arg)
)

(defun calc-tanh (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-tan arg)
)

(defun calc-arctanh (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-hyperbolic-func)
  (calc-tan arg)
)

(defun calc-arctan2 ()
  (interactive)
  (calc-slow-wrapper
   (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2))))
)

(defun calc-conj (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "conj" 'calcFunc-conj arg))
)

(defun calc-imaginary ()
  (interactive)
  (calc-slow-wrapper
   (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1))))
)



(defun calc-to-degrees (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op ">deg" 'calcFunc-deg arg))
)

(defun calc-to-radians (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op ">rad" 'calcFunc-rad arg))
)


(defun calc-degrees-mode (arg)
  (interactive "p")
  (cond ((= arg 1)
	 (calc-wrapper
	  (calc-change-mode 'calc-angle-mode 'deg)
	  (message "Angles measured in degrees.")))
	((= arg 2) (calc-radians-mode))
	((= arg 3) (calc-hms-mode))
	(t (error "Prefix argument out of range")))
)

(defun calc-radians-mode ()
  (interactive)
  (calc-wrapper
   (calc-change-mode 'calc-angle-mode 'rad)
   (message "Angles measured in radians."))
)


;;; 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!

	((not (math-natnump a))
	 (math-reject-arg a 'natnump))
	((integerp a)
	 (math-isqrt-small a))
	(t
	 (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a)))))))
)

(defun calcFunc-isqrt (a)
  (if (math-realp a)
      (math-isqrt (math-floor a))
    (math-floor (math-sqrt a)))
)


;;; This returns (flag . result) where the flag is T if A is a perfect square.
(defun math-isqrt-bignum (a)   ; [P.l L]
  (let ((len (length a)))

      (let* ((top (nth (1- len) a)))
	(math-isqrt-bignum-iter
	 a
	 (math-scale-bignum-3
	  (list (1+ (math-isqrt-small top)))
	  (/ len 2))))))
)

(defun math-isqrt-bignum-iter (a guess)   ; [l L l]
  (math-working "isqrt" (cons 'bigpos guess))
  (let* ((q (math-div-bignum a guess))
	 (s (math-add-bignum (car q) guess))

    (if (< comp 0)
	(math-isqrt-bignum-iter a g2)
      (cons (and (= comp 0)
		 (math-zerop-bignum (cdr q))
		 (= (% (car s) 2) 0))
	    guess)))
)

(defun math-zerop-bignum (a)
  (and (eq (car a) 0)
       (progn
	 (while (eq (car (setq a (cdr a))) 0))
	 (null a)))
)

(defun math-scale-bignum-3 (a n)   ; [L L S]
  (while (> n 0)
    (setq a (cons 0 a)
	  n (1- n)))
  a
)

(defun math-isqrt-small (a)   ; A > 0.  [S S]
  (let ((g (cond ((>= a 10000) 1000)
		 ((>= a 100) 100)
		 (t 10)))
	g2)
    (while (< (setq g2 (/ (+ g (/ a g)) 2)) g)
      (setq g g2))
    g)
)




;;; Compute the square root of a number.

   (let ((inf (math-infinitep a)))
     (and inf
	  (math-mul (math-sqrt (math-infinite-dir a inf)) inf)))
   (progn
     (calc-record-why 'numberp a)
     (list 'calcFunc-sqrt a)))
)
(fset 'calcFunc-sqrt (symbol-function 'math-sqrt))

(defun math-infinite-dir (a &optional inf)
  (or inf (setq inf (math-infinitep a)))
  (math-normalize (math-expr-subst a inf 1))
)

(defun math-sqrt-float (a &optional guess)   ; [F F F]
  (if calc-symbolic-mode
      (signal 'inexact-result nil)
    (math-with-extra-prec 1 (math-sqrt-raw a guess)))
)

(defun math-sqrt-raw (a &optional guess)   ; [F F F]
  (if (not (Math-posp a))
      (math-sqrt a)
    (if (null guess)
	(let ((ldiff (- (math-numdigs (nth 1 a)) 6)))
	  (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff)))
	  (setq guess (math-make-float (math-isqrt-small
					(math-scale-int (nth 1 a) (- ldiff)))
				       (/ (+ (nth 2 a) ldiff) 2)))))
    (math-sqrt-float-iter a guess))
)

(defun math-sqrt-float-iter (a guess)   ; [F F F]
  (math-working "sqrt" guess)
  (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess))
			    '(float 5 -1))))
     (if (math-nearly-equal-float g2 guess)
	 g2
       (math-sqrt-float-iter a g2)))
)

;;; True if A and B differ only in the last digit of precision.  [P F F]
(defun math-nearly-equal-float (a b)
  (let ((ediff (- (nth 2 a) (nth 2 b))))
    (cond ((= ediff 0)   ;; Expanded out for speed

	   (setq ediff (math-add (Math-integer-neg (nth 1 a))
				 (math-scale-int (nth 1 b) 1)))
	   (and (not (consp ediff))
		(< ediff 10)
		(> ediff -10)
		(= (math-numdigs (nth 1 a)) calc-internal-prec)))))
)

(defun math-nearly-equal (a b)   ;  [P N N] [Public]
  (setq a (math-float a))
  (setq b (math-float b))
  (if (eq (car a) 'polar) (setq a (math-complex a)))

	(and (math-nearly-equal-float (nth 1 a) b)
	     (math-nearly-zerop-float (nth 2 a) b)))
      (if (eq (car b) 'cplx)
	  (and (math-nearly-equal-float a (nth 1 b))
	       (math-nearly-zerop-float a (nth 2 b)))
	(math-nearly-equal-float a b)))
)

;;; True if A is nearly zero compared to B.  [P F F]
(defun math-nearly-zerop-float (a b)
  (or (eq (nth 1 a) 0)
      (<= (+ (math-numdigs (nth 1 a)) (nth 2 a))
	  (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec))))
)

(defun math-nearly-zerop (a b)   ; [P N R] [Public]
  (setq a (math-float a))
  (setq b (math-float b))
  (if (eq (car a) 'cplx)
      (and (math-nearly-zerop-float (nth 1 a) b)
	   (math-nearly-zerop-float (nth 2 a) b))
    (if (eq (car a) 'polar)
	(math-nearly-zerop-float (nth 1 a) b)
      (math-nearly-zerop-float a b)))
)

;;; This implementation could be improved, accuracy-wise.
(defun math-hypot (a b)
  (cond ((Math-zerop a) (math-abs b))
	((Math-zerop b) (math-abs a))

	     (math-to-hms (math-hypot (math-from-hms a 'deg)
				      (math-from-hms b 'deg)))
	   (math-to-hms (math-hypot (math-from-hms a 'deg) b))))
	((eq (car-safe b) 'hms)
	 (math-to-hms (math-hypot a (math-from-hms b 'deg))))
	(t nil))
)
(fset 'calcFunc-hypot (symbol-function 'math-hypot))

(defun calcFunc-sqr (x)
  (math-pow x 2)
)



(defun math-nth-root (a n)
  (cond ((= n 2) (math-sqrt a))

	 (and (not calc-symbolic-mode)
	      (math-nth-root-float a n)))
	((eq (car-safe a) 'polar)
	 (let ((root (math-nth-root (nth 1 a) n)))
	   (and root (list 'polar root (math-div (nth 2 a) n)))))
	(t nil))
)

(defun math-nth-root-float (a n &optional guess)
  (math-inexact-result)
  (math-with-extra-prec 1
    (let ((nf (math-float n))

      (math-nth-root-float-iter a (or guess
				      (math-make-float
				       1 (/ (+ (math-numdigs (nth 1 a))
					       (nth 2 a)
					       (/ n 2))
					    n))))))
)

(defun math-nth-root-float-iter (a guess)   ; uses "n", "nf", "nfm1"
  (math-working "root" guess)
  (let ((g2 (math-div-float (math-add-float (math-mul nfm1 guess)
					    (math-div-float
					     a (math-ipow guess (1- n))))
			    nf)))
    (if (math-nearly-equal-float g2 guess)
	g2
      (math-nth-root-float-iter a g2)))
)

(defun math-nth-root-integer (a n &optional guess)   ; [I I S]
  (math-nth-root-int-iter a (or guess
				(math-scale-int 1 (/ (+ (math-numdigs a)
							(1- n))
						     n))))
)

(defun math-nth-root-int-iter (a guess)   ; uses "n"
  (math-working "root" guess)
  (let* ((q (math-idivmod a (math-ipow guess (1- n))))
	 (s (math-add (car q) (math-mul (1- n) guess)))

    (if (Math-natnum-lessp (car g2) guess)
	(math-nth-root-int-iter a (car g2))
      (cons (and (equal (car g2) guess)
		 (eq (cdr q) 0)
		 (eq (cdr g2) 0))
	    guess)))
)

(defun calcFunc-nroot (x n)
  (calcFunc-pow x (if (integerp n)
		      (math-make-frac 1 n)
		    (math-div 1 n)))
)




;;;; Transcendental functions.

(defun math-to-radians (a)   ; [N N]
  (cond ((eq (car-safe a) 'hms)
	 (math-from-hms a 'rad))
	((memq calc-angle-mode '(deg hms))
	 (math-mul a (math-pi-over-180)))
	(t a))
)

(defun math-from-radians (a)   ; [N N]
  (cond ((eq calc-angle-mode 'deg)
	 (if (math-constp a)
	     (math-div a (math-pi-over-180))
	   (list 'calcFunc-deg a)))
	((eq calc-angle-mode 'hms)
	 (math-to-hms a 'rad))
	(t a))
)

(defun math-to-radians-2 (a)   ; [N N]
  (cond ((eq (car-safe a) 'hms)
	 (math-from-hms a 'rad))
	((memq calc-angle-mode '(deg hms))
	 (if calc-symbolic-mode
	     (math-div (math-mul a '(var pi var-pi)) 180)
	   (math-mul a (math-pi-over-180))))
	(t a))
)

(defun math-from-radians-2 (a)   ; [N N]
  (cond ((memq calc-angle-mode '(deg hms))
	 (if calc-symbolic-mode
	     (math-div (math-mul 180 a) '(var pi var-pi))
	   (math-div a (math-pi-over-180))))
	(t a))
)



;;; Sine, cosine, and tangent.

	((and (eq (car x) 'intv) (math-intv-constp x))
	 (calcFunc-cos (math-sub x (math-quarter-circle nil))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'scalarp x)
	   (list 'calcFunc-sin x)))
)

(defun calcFunc-cos (x)   ; [N N] [Public]
  (cond ((and (integerp x)
	      (if (eq calc-angle-mode 'deg)
		  (= (% x 90) 0)

		     int))
	       (list 'intv 3 -1 1)))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'scalarp x)
	   (list 'calcFunc-cos x)))
)

(defun calcFunc-sincos (x)   ; [V N] [Public]
  (if (Math-scalarp x)
      (math-with-extra-prec 2
	(let ((sc (math-sin-cos-raw (math-to-radians (math-float x)))))
	  (list 'vec (cdr sc) (car sc))))    ; the vector [cos, sin]
    (list 'vec (calcFunc-sin x) (calcFunc-cos x)))
)

(defun calcFunc-tan (x)   ; [N N] [Public]
  (cond ((and (integerp x)
	      (if (eq calc-angle-mode 'deg)
		  (= (% x 180) 0)

				      (math-tan-raw (nth 3 xx))))))
	     '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'scalarp x)
	   (list 'calcFunc-tan x)))
)

(defun math-sin-raw (x)   ; [N N]
  (cond ((eq (car x) 'cplx)
	 (let* ((expx (math-exp-raw (nth 2 x)))
		(expmx (math-div-float '(float 1 0) expx))

	 (math-polar (math-sin-raw (math-complex x))))
	((Math-integer-negp (nth 1 x))
	 (math-neg-float (math-sin-raw (math-neg-float x))))
	((math-lessp-float '(float 7 0) x)  ; avoid inf loops due to roundoff
	 (math-sin-raw (math-mod x (math-two-pi))))
	(t (math-sin-raw-2 x x)))
)

(defun math-cos-raw (x)   ; [N N]
  (if (eq (car-safe x) 'polar)
      (math-polar (math-cos-raw (math-complex x)))
    (math-sin-raw (math-sub (math-pi-over-2) 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))
  (cons (math-sin-raw x) (math-cos-raw x))
)

(defun math-tan-raw (x)   ; [N N]
  (cond ((eq (car x) 'cplx)
	 (let* ((x (math-mul x '(float 2 0)))
		(expx (math-exp-raw (nth 2 x)))

	 (math-polar (math-tan-raw (math-complex x))))
	(t
	 (let ((sc (math-sin-cos-raw x)))
	   (if (eq (nth 1 (cdr sc)) 0)
	       (math-div (car sc) 0)
	     (math-div-float (car sc) (cdr sc))))))
)

(defun math-sin-raw-2 (x orgx)   ; This avoids poss of inf recursion.  [F F]
  (let ((xmpo2 (math-sub-float (math-pi-over-2) x)))
    (cond ((Math-integer-negp (nth 1 xmpo2))
	   (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi))

	   (math-cos-raw-2 xmpo2 orgx))
	  ((math-lessp-float x (math-neg (math-pi-over-4)))
	   (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx)))
	  ((math-nearly-zerop-float x orgx) '(float 0 0))
	  (calc-symbolic-mode (signal 'inexact-result nil))
	  (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x))))))
)

(defun math-cos-raw-2 (x orgx)   ; [F F]
  (cond ((math-nearly-zerop-float x orgx) '(float 1 0))
	(calc-symbolic-mode (signal 'inexact-result nil))
	(t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
	     (math-sin-series
	      (math-add-float '(float 1 0)
			      (math-mul-float xnegsqr '(float 5 -1)))
	      24 5 xnegsqr xnegsqr))))
)

(defun math-sin-series (sum nfac n x xnegsqr)
  (math-working "sin" sum)
  (let* ((nextx (math-mul-float x xnegsqr))
	 (nextsum (math-add-float sum (math-div-float nextx
						      (math-float nfac)))))
    (if (math-nearly-equal-float sum nextsum)
	sum
      (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
		       (+ n 2) nextx xnegsqr)))
)


;;; Inverse sine, cosine, tangent.

(defun calcFunc-arcsin (x)   ; [N N] [Public]

			 (calcFunc-arcsin (nth 2 x))
			 (calcFunc-arcsin (nth 3 x))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arcsin x)))
)

(defun calcFunc-arccos (x)   ; [N N] [Public]
  (cond ((eq x 1) 0)
	((and (eq x 0) (eq calc-angle-mode 'deg)) 90)
	((and (eq x -1) (eq calc-angle-mode 'deg)) 180)

			 (calcFunc-arccos (nth 2 x))
			 (calcFunc-arccos (nth 3 x))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arccos x)))
)

(defun calcFunc-arctan (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	((and (eq x 1) (eq calc-angle-mode 'deg)) 45)
	((and (eq x -1) (eq calc-angle-mode 'deg)) -45)

	((equal x '(neg (var inf var-inf)))
	 (math-neg (math-quarter-circle t)))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arctan x)))
)

(defun math-arcsin-raw (x)   ; [N N]
  (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x)))))
    (if (or (memq (car x) '(cplx polar))
	    (memq (car a) '(cplx polar)))
	(math-with-extra-prec 2   ; use extra precision for difficult case
	  (math-mul '(cplx 0 -1)
		    (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a))))
      (math-arctan2-raw x a)))
)

(defun math-arccos-raw (x)   ; [N N]
  (math-sub (math-pi-over-2) (math-arcsin-raw x))
)

(defun math-arctan-raw (x)   ; [N N]
  (cond ((memq (car x) '(cplx polar))
	 (math-with-extra-prec 2   ; extra-extra
	   (math-div (math-sub

	   (math-sub-float (math-mul-float (math-pi) '(float 25 -2))
			   (math-arctan-raw (math-div-float
					     (math-sub-float '(float 1 0) x)
					     (math-add-float '(float 1 0)
							     x))))))
	(t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x)))))
)

(defun math-arctan-series (sum n x xnegsqr)
  (math-working "arctan" sum)
  (let* ((nextx (math-mul-float x xnegsqr))
	 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
    (if (math-nearly-equal-float sum nextsum)
	sum
      (math-arctan-series nextsum (+ n 2) nextx xnegsqr)))
)

(defun calcFunc-arctan2 (y x)   ; [F R R] [Public]
  (if (Math-anglep y)
      (if (Math-anglep x)
	  (math-with-extra-prec 2

		  (setq y nil))))
	  (if (and y x)
	      (calcFunc-arctan2 y x)
	    '(var nan var-nan)))
      (calc-record-why 'anglep y)
      (list 'calcFunc-arctan2 y x)))
)

(defun math-arctan2-raw (y x)   ; [F R R]
  (cond ((math-zerop y)
	 (if (math-negp x) (math-pi)
	   (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0)))

	((math-posp y)
	 (math-add-float (math-arctan-raw (math-div-float y x))
			 (math-pi)))
	(t
	 (math-sub-float (math-arctan-raw (math-div-float y x))
			 (math-pi))))
)

(defun calcFunc-arcsincos (x)   ; [V N] [Public]
  (if (and (Math-vectorp x)
	   (= (length x) 3))
      (calcFunc-arctan2 (nth 2 x) (nth 1 x))
    (math-reject-arg x "*Two-element vector expected"))
)



;;; Exponential function.

	((equal x '(neg (var inf var-inf)))
	 0)
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-exp x)))
)

(defun calcFunc-expm1 (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	((math-zerop x) '(float 0 0))
	(calc-symbolic-mode (signal 'inexact-result nil))

	((equal x '(neg (var inf var-inf)))
	 -1)
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-expm1 x)))
)

(defun calcFunc-exp10 (x)   ; [N N] [Public]
  (if (eq x 0)
      1
    (math-pow '(float 1 1) x))
)

(defun math-exp-raw (x)   ; [N N]
  (cond ((math-zerop x) '(float 1 0))
	(calc-symbolic-mode (signal 'inexact-result nil))
	((eq (car x) 'cplx)

		(hfrac (math-sub-float x (math-mul-float (math-float hint)
							 '(float 5 -1)))))
	   (math-mul-float (math-ipow (math-sqrt-e) hint)
			   (math-add-float '(float 1 0)
					   (math-exp-minus-1-raw hfrac)))))
	(t (math-add-float '(float 1 0) (math-exp-minus-1-raw x))))
)

(defun math-exp-minus-1-raw (x)   ; [F F]
  (math-exp-series x 2 3 x x)
)

(defun math-exp-series (sum nfac n xpow x)
  (math-working "exp" sum)
  (let* ((nextx (math-mul-float xpow x))
	 (nextsum (math-add-float sum (math-div-float nextx
						      (math-float nfac)))))
    (if (math-nearly-equal-float sum nextsum)
	sum
      (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x)))
)



;;; Logarithms.

	((math-infinitep x)
	 (if (equal x '(var nan var-nan))
	     x
	   '(var inf var-inf)))
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-ln x)))
)

(defun calcFunc-log10 (x)   ; [N N] [Public]
  (cond ((math-equal-int x 1)
	 (if (math-floatp x) '(float 0 0) 0))
	((and (Math-integerp x)

	((math-infinitep x)
	 (if (equal x '(var nan var-nan))
	     x
	   '(var inf var-inf)))
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-log10 x)))
)

(defun calcFunc-log (x &optional b)   ; [N N N] [Public]
  (cond ((or (null b) (equal b '(var e var-e)))
	 (math-normalize (list 'calcFunc-ln x)))
	((or (eq b 10) (equal b '(float 1 1)))

	     (math-infinitep b))
	 (math-div (calcFunc-ln x) (calcFunc-ln b)))
	(t (if (Math-numberp b)
	       (calc-record-why 'numberp x)
	     (calc-record-why 'numberp b))
	   (list 'calcFunc-log x b)))
)

(defun calcFunc-alog (x &optional b)
  (cond ((or (null b) (equal b '(var e var-e)))
	 (math-normalize (list 'calcFunc-exp x)))
	(t (math-pow b x)))
)

(defun calcFunc-ilog (x b)
  (if (and (math-natnump x) (not (eq x 0))
	   (math-natnump b) (not (eq b 0)))
      (if (eq b 1)
	  (math-reject-arg x "*Logarithm base one")
	(if (Math-natnum-lessp x b)
	    0
	  (cdr (math-integer-log x b))))
    (math-floor (calcFunc-log x b)))
)

(defun math-integer-log (x b)
  (let ((pows (list b))
	(pow (math-sqr b))
	next

      (setq n (/ n 2)
	    next (math-mul pow (car pows)))
      (or (Math-lessp x next)
	  (setq pow next
		sum (+ sum n))))
    (cons (equal pow x) sum))
)


(defun math-log-base-raw (b)   ; [N N]
  (if (not (and (equal (car math-log-base-cache) b)
		(eq (nth 1 math-log-base-cache) calc-internal-prec)))
      (setq math-log-base-cache (list b calc-internal-prec
				      (math-ln-raw (math-float b)))))
  (nth 2 math-log-base-cache)
)
(setq math-log-base-cache nil)

(defun calcFunc-lnp1 (x)   ; [N N] [Public]
  (cond ((Math-equal-int x -1)
	 (if calc-infinite-mode

	((math-infinitep x)
	 (if (equal x '(var nan var-nan))
	     x
	   '(var inf var-inf)))
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-lnp1 x)))
)

(defun math-ln-raw (x)    ; [N N] --- must be float format!
  (cond ((eq (car-safe x) 'cplx)
	 (list 'cplx
	       (math-mul-float (math-ln-raw

	  (list 'cplx   ; negative and real
		(math-ln-raw (math-neg-float x))
		(math-pi))))
	(t (list 'cplx   ; negative and real
		 (math-ln-raw (math-neg-float x))
		 (math-pi))))
)

(defun math-ln-raw-2 (x)    ; [F F]
  (cond ((math-lessp-float '(float 14 -1) x)
	 (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1)))
			 (math-ln-2)))
	(t    ; now .7 < x <= 1.4
	 (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
					(math-add-float x '(float 1 0))))))
)

(defun math-ln-raw-3 (x)   ; [F F]
  (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
		  '(float 2 0))
)

;;; Compute ln((1+x)/(1-x))
(defun math-ln-raw-series (sum n x xsqr)
  (math-working "log" sum)
  (let* ((nextx (math-mul-float x xsqr))
	 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
    (if (math-nearly-equal-float sum nextsum)
	sum
      (math-ln-raw-series nextsum (+ n 2) nextx xsqr)))
)

(defun math-ln-plus-1-raw (x)
  (math-lnp1-series x 2 x (math-neg x))
)

(defun math-lnp1-series (sum n xpow x)
  (math-working "lnp1" sum)
  (let* ((nextx (math-mul-float xpow x))
	 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
    (if (math-nearly-equal-float sum nextsum)
	sum
      (math-lnp1-series nextsum (1+ n) nextx x)))
)

(math-defcache math-ln-10 (float (bigpos 018 684 045 994 092 585 302 2) -21)
  (math-ln-raw-2 '(float 1 1)))

(math-defcache math-ln-2 (float (bigpos 417 309 945 559 180 147 693) -21)

	((or (equal x '(var inf var-inf))
	     (equal x '(neg (var inf var-inf)))
	     (equal x '(var nan var-nan)))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-sinh x)))
)
(put 'calcFunc-sinh 'math-expandable t)

(defun calcFunc-cosh (x)   ; [N N] [Public]
  (cond ((eq x 0) 1)
	(math-expand-formulas

	((or (equal x '(var inf var-inf))
	     (equal x '(neg (var inf var-inf)))
	     (equal x '(var nan var-nan)))
	 (math-abs x))
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-cosh x)))
)
(put 'calcFunc-cosh 'math-expandable t)

(defun calcFunc-tanh (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	(math-expand-formulas

	((equal x '(neg (var inf var-inf)))
	 -1)
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-tanh x)))
)
(put 'calcFunc-tanh 'math-expandable t)

(defun calcFunc-arcsinh (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	(math-expand-formulas

	((or (equal x '(var inf var-inf))
	     (equal x '(neg (var inf var-inf)))
	     (equal x '(var nan var-nan)))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arcsinh x)))
)
(put 'calcFunc-arcsinh 'math-expandable t)

(defun calcFunc-arccosh (x)   ; [N N] [Public]
  (cond ((eq x 1) 0)
	((and (eq x -1) calc-symbolic-mode)

	((or (equal x '(var inf var-inf))
	     (equal x '(neg (var inf var-inf)))
	     (equal x '(var nan var-nan)))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arccosh x)))
)
(put 'calcFunc-arccosh 'math-expandable t)

(defun calcFunc-arctanh (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	((and (Math-equal-int x 1) calc-infinite-mode)

			 (calcFunc-arctanh (nth 2 x))
			 (calcFunc-arctanh (nth 3 x))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arctanh x)))
)
(put 'calcFunc-arctanh 'math-expandable t)


;;; Convert A from HMS or degrees to radians.
(defun calcFunc-rad (a)   ; [R R] [Public]

	 (math-make-sdev (calcFunc-rad (nth 1 a))
			 (calcFunc-rad (nth 2 a))))
	(math-expand-formulas
	 (math-div (math-mul a '(var pi var-pi)) 180))
	((math-infinitep a) a)
	(t (list 'calcFunc-rad a)))
)
(put 'calcFunc-rad 'math-expandable t)

;;; Convert A from HMS or radians to degrees.
(defun calcFunc-deg (a)   ; [R R] [Public]
  (cond ((or (Math-numberp a)

	 (math-make-sdev (calcFunc-deg (nth 1 a))
			 (calcFunc-deg (nth 2 a))))
	(math-expand-formulas
	 (math-div (math-mul 180 a) '(var pi var-pi)))
	((math-infinitep a) a)
	(t (list 'calcFunc-deg a)))
)
(put 'calcFunc-deg 'math-expandable t)

41044:4549dec29728

;; Calculator for GNU Emacs, part II [calc-math.el]
;; Copyright (C) 1990, 1991, 1992, 1993, 2001 Free Software Foundation, Inc.
;; Written by Dave Gillespie, [e-mail].

;; This file is part of GNU Emacs.

;; GNU Emacs is distributed in the hope that it will be useful,

(defun calc-sqrt (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-unary-op "^2" 'calcFunc-sqr arg)
     (calc-unary-op "sqrt" 'calcFunc-sqrt arg))))

(defun calc-isqrt (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-unary-op "^2" 'calcFunc-sqr arg)
     (calc-unary-op "isqt" 'calcFunc-isqrt arg))))


(defun calc-hypot (arg)
  (interactive "P")
  (calc-slow-wrapper
   (calc-binary-op "hypt" 'calcFunc-hypot arg)))

(defun calc-ln (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-exp arg))

(defun calc-log10 (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-ln arg))

(defun calc-log (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-binary-op "alog" 'calcFunc-alog arg)
     (calc-binary-op "log" 'calcFunc-log arg))))

(defun calc-ilog (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-binary-op "alog" 'calcFunc-alog arg)
     (calc-binary-op "ilog" 'calcFunc-ilog arg))))

(defun calc-lnp1 (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-expm1 arg))

(defun calc-exp (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (if (calc-is-inverse)
	   (calc-unary-op "lg10" 'calcFunc-log10 arg)
	 (calc-unary-op "10^" 'calcFunc-exp10 arg))
     (if (calc-is-inverse)
	 (calc-unary-op "ln" 'calcFunc-ln arg)
       (calc-unary-op "exp" 'calcFunc-exp arg)))))

(defun calc-expm1 (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-unary-op "ln+1" 'calcFunc-lnp1 arg)
     (calc-unary-op "ex-1" 'calcFunc-expm1 arg))))

(defun calc-pi ()
  (interactive)
  (calc-slow-wrapper
   (if (calc-is-inverse)

	 (if calc-symbolic-mode
	     (calc-pop-push-record 0 "e" '(var e var-e))
	   (calc-pop-push-record 0 "e" (math-e)))
       (if calc-symbolic-mode
	   (calc-pop-push-record 0 "pi" '(var pi var-pi))
	 (calc-pop-push-record 0 "pi" (math-pi)))))))

(defun calc-sin (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (if (calc-is-inverse)
	   (calc-unary-op "asnh" 'calcFunc-arcsinh arg)
	 (calc-unary-op "sinh" 'calcFunc-sinh arg))
     (if (calc-is-inverse)
	 (calc-unary-op "asin" 'calcFunc-arcsin arg)
       (calc-unary-op "sin" 'calcFunc-sin arg)))))

(defun calc-arcsin (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-sin arg))

(defun calc-sinh (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-sin arg))

(defun calc-arcsinh (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-hyperbolic-func)
  (calc-sin arg))

(defun calc-cos (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (if (calc-is-inverse)
	   (calc-unary-op "acsh" 'calcFunc-arccosh arg)
	 (calc-unary-op "cosh" 'calcFunc-cosh arg))
     (if (calc-is-inverse)
	 (calc-unary-op "acos" 'calcFunc-arccos arg)
       (calc-unary-op "cos" 'calcFunc-cos arg)))))

(defun calc-arccos (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-cos arg))

(defun calc-cosh (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-cos arg))

(defun calc-arccosh (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-hyperbolic-func)
  (calc-cos arg))

(defun calc-sincos ()
  (interactive)
  (calc-slow-wrapper
   (if (calc-is-inverse)
       (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1)))
     (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1))))))

(defun calc-tan (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (if (calc-is-inverse)
	   (calc-unary-op "atnh" 'calcFunc-arctanh arg)
	 (calc-unary-op "tanh" 'calcFunc-tanh arg))
     (if (calc-is-inverse)
	 (calc-unary-op "atan" 'calcFunc-arctan arg)
       (calc-unary-op "tan" 'calcFunc-tan arg)))))

(defun calc-arctan (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-tan arg))

(defun calc-tanh (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-tan arg))

(defun calc-arctanh (arg)
  (interactive "P")
  (calc-invert-func)
  (calc-hyperbolic-func)
  (calc-tan arg))

(defun calc-arctan2 ()
  (interactive)
  (calc-slow-wrapper
   (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2)))))

(defun calc-conj (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op "conj" 'calcFunc-conj arg)))

(defun calc-imaginary ()
  (interactive)
  (calc-slow-wrapper
   (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))))



(defun calc-to-degrees (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op ">deg" 'calcFunc-deg arg)))

(defun calc-to-radians (arg)
  (interactive "P")
  (calc-wrapper
   (calc-unary-op ">rad" 'calcFunc-rad arg)))


(defun calc-degrees-mode (arg)
  (interactive "p")
  (cond ((= arg 1)
	 (calc-wrapper
	  (calc-change-mode 'calc-angle-mode 'deg)
	  (message "Angles measured in degrees.")))
	((= arg 2) (calc-radians-mode))
	((= arg 3) (calc-hms-mode))
	(t (error "Prefix argument out of range"))))

(defun calc-radians-mode ()
  (interactive)
  (calc-wrapper
   (calc-change-mode 'calc-angle-mode 'rad)
   (message "Angles measured in radians.")))


;;; 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!

	((not (math-natnump a))
	 (math-reject-arg a 'natnump))
	((integerp a)
	 (math-isqrt-small a))
	(t
	 (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a))))))))

(defun calcFunc-isqrt (a)
  (if (math-realp a)
      (math-isqrt (math-floor a))
    (math-floor (math-sqrt a))))


;;; This returns (flag . result) where the flag is T if A is a perfect square.
(defun math-isqrt-bignum (a)   ; [P.l L]
  (let ((len (length a)))

      (let* ((top (nth (1- len) a)))
	(math-isqrt-bignum-iter
	 a
	 (math-scale-bignum-3
	  (list (1+ (math-isqrt-small top)))
	  (/ len 2)))))))

(defun math-isqrt-bignum-iter (a guess)   ; [l L l]
  (math-working "isqrt" (cons 'bigpos guess))
  (let* ((q (math-div-bignum a guess))
	 (s (math-add-bignum (car q) guess))

    (if (< comp 0)
	(math-isqrt-bignum-iter a g2)
      (cons (and (= comp 0)
		 (math-zerop-bignum (cdr q))
		 (= (% (car s) 2) 0))
	    guess))))

(defun math-zerop-bignum (a)
  (and (eq (car a) 0)
       (progn
	 (while (eq (car (setq a (cdr a))) 0))
	 (null a))))

(defun math-scale-bignum-3 (a n)   ; [L L S]
  (while (> n 0)
    (setq a (cons 0 a)
	  n (1- n)))
  a)

(defun math-isqrt-small (a)   ; A > 0.  [S S]
  (let ((g (cond ((>= a 10000) 1000)
		 ((>= a 100) 100)
		 (t 10)))
	g2)
    (while (< (setq g2 (/ (+ g (/ a g)) 2)) g)
      (setq g g2))
    g))




;;; Compute the square root of a number.

   (let ((inf (math-infinitep a)))
     (and inf
	  (math-mul (math-sqrt (math-infinite-dir a inf)) inf)))
   (progn
     (calc-record-why 'numberp a)
     (list 'calcFunc-sqrt a))))
(defalias calcFunc-sqrt 'math-sqrt)

(defun math-infinite-dir (a &optional inf)
  (or inf (setq inf (math-infinitep a)))
  (math-normalize (math-expr-subst a inf 1)))

(defun math-sqrt-float (a &optional guess)   ; [F F F]
  (if calc-symbolic-mode
      (signal 'inexact-result nil)
    (math-with-extra-prec 1 (math-sqrt-raw a guess))))

(defun math-sqrt-raw (a &optional guess)   ; [F F F]
  (if (not (Math-posp a))
      (math-sqrt a)
    (if (null guess)
	(let ((ldiff (- (math-numdigs (nth 1 a)) 6)))
	  (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff)))
	  (setq guess (math-make-float (math-isqrt-small
					(math-scale-int (nth 1 a) (- ldiff)))
				       (/ (+ (nth 2 a) ldiff) 2)))))
    (math-sqrt-float-iter a guess)))

(defun math-sqrt-float-iter (a guess)   ; [F F F]
  (math-working "sqrt" guess)
  (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess))
			    '(float 5 -1))))
     (if (math-nearly-equal-float g2 guess)
	 g2
       (math-sqrt-float-iter a g2))))

;;; True if A and B differ only in the last digit of precision.  [P F F]
(defun math-nearly-equal-float (a b)
  (let ((ediff (- (nth 2 a) (nth 2 b))))
    (cond ((= ediff 0)   ;; Expanded out for speed

	   (setq ediff (math-add (Math-integer-neg (nth 1 a))
				 (math-scale-int (nth 1 b) 1)))
	   (and (not (consp ediff))
		(< ediff 10)
		(> ediff -10)
		(= (math-numdigs (nth 1 a)) calc-internal-prec))))))

(defun math-nearly-equal (a b)   ;  [P N N] [Public]
  (setq a (math-float a))
  (setq b (math-float b))
  (if (eq (car a) 'polar) (setq a (math-complex a)))

	(and (math-nearly-equal-float (nth 1 a) b)
	     (math-nearly-zerop-float (nth 2 a) b)))
      (if (eq (car b) 'cplx)
	  (and (math-nearly-equal-float a (nth 1 b))
	       (math-nearly-zerop-float a (nth 2 b)))
	(math-nearly-equal-float a b))))

;;; True if A is nearly zero compared to B.  [P F F]
(defun math-nearly-zerop-float (a b)
  (or (eq (nth 1 a) 0)
      (<= (+ (math-numdigs (nth 1 a)) (nth 2 a))
	  (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec)))))

(defun math-nearly-zerop (a b)   ; [P N R] [Public]
  (setq a (math-float a))
  (setq b (math-float b))
  (if (eq (car a) 'cplx)
      (and (math-nearly-zerop-float (nth 1 a) b)
	   (math-nearly-zerop-float (nth 2 a) b))
    (if (eq (car a) 'polar)
	(math-nearly-zerop-float (nth 1 a) b)
      (math-nearly-zerop-float a b))))

;;; This implementation could be improved, accuracy-wise.
(defun math-hypot (a b)
  (cond ((Math-zerop a) (math-abs b))
	((Math-zerop b) (math-abs a))

	     (math-to-hms (math-hypot (math-from-hms a 'deg)
				      (math-from-hms b 'deg)))
	   (math-to-hms (math-hypot (math-from-hms a 'deg) b))))
	((eq (car-safe b) 'hms)
	 (math-to-hms (math-hypot a (math-from-hms b 'deg))))
	(t nil)))
(defalias calcFunc-hypot 'math-hypot)

(defun calcFunc-sqr (x)
  (math-pow x 2))



(defun math-nth-root (a n)
  (cond ((= n 2) (math-sqrt a))

	 (and (not calc-symbolic-mode)
	      (math-nth-root-float a n)))
	((eq (car-safe a) 'polar)
	 (let ((root (math-nth-root (nth 1 a) n)))
	   (and root (list 'polar root (math-div (nth 2 a) n)))))
	(t nil)))

(defun math-nth-root-float (a n &optional guess)
  (math-inexact-result)
  (math-with-extra-prec 1
    (let ((nf (math-float n))

      (math-nth-root-float-iter a (or guess
				      (math-make-float
				       1 (/ (+ (math-numdigs (nth 1 a))
					       (nth 2 a)
					       (/ n 2))
					    n)))))))

(defun math-nth-root-float-iter (a guess)   ; uses "n", "nf", "nfm1"
  (math-working "root" guess)
  (let ((g2 (math-div-float (math-add-float (math-mul nfm1 guess)
					    (math-div-float
					     a (math-ipow guess (1- n))))
			    nf)))
    (if (math-nearly-equal-float g2 guess)
	g2
      (math-nth-root-float-iter a g2))))

(defun math-nth-root-integer (a n &optional guess)   ; [I I S]
  (math-nth-root-int-iter a (or guess
				(math-scale-int 1 (/ (+ (math-numdigs a)
							(1- n))
						     n)))))

(defun math-nth-root-int-iter (a guess)   ; uses "n"
  (math-working "root" guess)
  (let* ((q (math-idivmod a (math-ipow guess (1- n))))
	 (s (math-add (car q) (math-mul (1- n) guess)))

    (if (Math-natnum-lessp (car g2) guess)
	(math-nth-root-int-iter a (car g2))
      (cons (and (equal (car g2) guess)
		 (eq (cdr q) 0)
		 (eq (cdr g2) 0))
	    guess))))

(defun calcFunc-nroot (x n)
  (calcFunc-pow x (if (integerp n)
		      (math-make-frac 1 n)
		    (math-div 1 n))))




;;;; Transcendental functions.

(defun math-to-radians (a)   ; [N N]
  (cond ((eq (car-safe a) 'hms)
	 (math-from-hms a 'rad))
	((memq calc-angle-mode '(deg hms))
	 (math-mul a (math-pi-over-180)))
	(t a)))

(defun math-from-radians (a)   ; [N N]
  (cond ((eq calc-angle-mode 'deg)
	 (if (math-constp a)
	     (math-div a (math-pi-over-180))
	   (list 'calcFunc-deg a)))
	((eq calc-angle-mode 'hms)
	 (math-to-hms a 'rad))
	(t a)))

(defun math-to-radians-2 (a)   ; [N N]
  (cond ((eq (car-safe a) 'hms)
	 (math-from-hms a 'rad))
	((memq calc-angle-mode '(deg hms))
	 (if calc-symbolic-mode
	     (math-div (math-mul a '(var pi var-pi)) 180)
	   (math-mul a (math-pi-over-180))))
	(t a)))

(defun math-from-radians-2 (a)   ; [N N]
  (cond ((memq calc-angle-mode '(deg hms))
	 (if calc-symbolic-mode
	     (math-div (math-mul 180 a) '(var pi var-pi))
	   (math-div a (math-pi-over-180))))
	(t a)))



;;; Sine, cosine, and tangent.

	((and (eq (car x) 'intv) (math-intv-constp x))
	 (calcFunc-cos (math-sub x (math-quarter-circle nil))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'scalarp x)
	   (list 'calcFunc-sin x))))

(defun calcFunc-cos (x)   ; [N N] [Public]
  (cond ((and (integerp x)
	      (if (eq calc-angle-mode 'deg)
		  (= (% x 90) 0)

		     int))
	       (list 'intv 3 -1 1)))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'scalarp x)
	   (list 'calcFunc-cos x))))

(defun calcFunc-sincos (x)   ; [V N] [Public]
  (if (Math-scalarp x)
      (math-with-extra-prec 2
	(let ((sc (math-sin-cos-raw (math-to-radians (math-float x)))))
	  (list 'vec (cdr sc) (car sc))))    ; the vector [cos, sin]
    (list 'vec (calcFunc-sin x) (calcFunc-cos x))))

(defun calcFunc-tan (x)   ; [N N] [Public]
  (cond ((and (integerp x)
	      (if (eq calc-angle-mode 'deg)
		  (= (% x 180) 0)

				      (math-tan-raw (nth 3 xx))))))
	     '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'scalarp x)
	   (list 'calcFunc-tan x))))

(defun math-sin-raw (x)   ; [N N]
  (cond ((eq (car x) 'cplx)
	 (let* ((expx (math-exp-raw (nth 2 x)))
		(expmx (math-div-float '(float 1 0) expx))

	 (math-polar (math-sin-raw (math-complex x))))
	((Math-integer-negp (nth 1 x))
	 (math-neg-float (math-sin-raw (math-neg-float x))))
	((math-lessp-float '(float 7 0) x)  ; avoid inf loops due to roundoff
	 (math-sin-raw (math-mod x (math-two-pi))))
	(t (math-sin-raw-2 x x))))

(defun math-cos-raw (x)   ; [N N]
  (if (eq (car-safe x) 'polar)
      (math-polar (math-cos-raw (math-complex x)))
    (math-sin-raw (math-sub (math-pi-over-2) 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))
  (cons (math-sin-raw x) (math-cos-raw x)))

(defun math-tan-raw (x)   ; [N N]
  (cond ((eq (car x) 'cplx)
	 (let* ((x (math-mul x '(float 2 0)))
		(expx (math-exp-raw (nth 2 x)))

	 (math-polar (math-tan-raw (math-complex x))))
	(t
	 (let ((sc (math-sin-cos-raw x)))
	   (if (eq (nth 1 (cdr sc)) 0)
	       (math-div (car sc) 0)
	     (math-div-float (car sc) (cdr sc)))))))

(defun math-sin-raw-2 (x orgx)   ; This avoids poss of inf recursion.  [F F]
  (let ((xmpo2 (math-sub-float (math-pi-over-2) x)))
    (cond ((Math-integer-negp (nth 1 xmpo2))
	   (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi))

	   (math-cos-raw-2 xmpo2 orgx))
	  ((math-lessp-float x (math-neg (math-pi-over-4)))
	   (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx)))
	  ((math-nearly-zerop-float x orgx) '(float 0 0))
	  (calc-symbolic-mode (signal 'inexact-result nil))
	  (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x)))))))

(defun math-cos-raw-2 (x orgx)   ; [F F]
  (cond ((math-nearly-zerop-float x orgx) '(float 1 0))
	(calc-symbolic-mode (signal 'inexact-result nil))
	(t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
	     (math-sin-series
	      (math-add-float '(float 1 0)
			      (math-mul-float xnegsqr '(float 5 -1)))
	      24 5 xnegsqr xnegsqr)))))

(defun math-sin-series (sum nfac n x xnegsqr)
  (math-working "sin" sum)
  (let* ((nextx (math-mul-float x xnegsqr))
	 (nextsum (math-add-float sum (math-div-float nextx
						      (math-float nfac)))))
    (if (math-nearly-equal-float sum nextsum)
	sum
      (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
		       (+ n 2) nextx xnegsqr))))


;;; Inverse sine, cosine, tangent.

(defun calcFunc-arcsin (x)   ; [N N] [Public]

			 (calcFunc-arcsin (nth 2 x))
			 (calcFunc-arcsin (nth 3 x))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arcsin x))))

(defun calcFunc-arccos (x)   ; [N N] [Public]
  (cond ((eq x 1) 0)
	((and (eq x 0) (eq calc-angle-mode 'deg)) 90)
	((and (eq x -1) (eq calc-angle-mode 'deg)) 180)

			 (calcFunc-arccos (nth 2 x))
			 (calcFunc-arccos (nth 3 x))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arccos x))))

(defun calcFunc-arctan (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	((and (eq x 1) (eq calc-angle-mode 'deg)) 45)
	((and (eq x -1) (eq calc-angle-mode 'deg)) -45)

	((equal x '(neg (var inf var-inf)))
	 (math-neg (math-quarter-circle t)))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arctan x))))

(defun math-arcsin-raw (x)   ; [N N]
  (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x)))))
    (if (or (memq (car x) '(cplx polar))
	    (memq (car a) '(cplx polar)))
	(math-with-extra-prec 2   ; use extra precision for difficult case
	  (math-mul '(cplx 0 -1)
		    (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a))))
      (math-arctan2-raw x a))))

(defun math-arccos-raw (x)   ; [N N]
  (math-sub (math-pi-over-2) (math-arcsin-raw x)))

(defun math-arctan-raw (x)   ; [N N]
  (cond ((memq (car x) '(cplx polar))
	 (math-with-extra-prec 2   ; extra-extra
	   (math-div (math-sub

	   (math-sub-float (math-mul-float (math-pi) '(float 25 -2))
			   (math-arctan-raw (math-div-float
					     (math-sub-float '(float 1 0) x)
					     (math-add-float '(float 1 0)
							     x))))))
	(t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x))))))

(defun math-arctan-series (sum n x xnegsqr)
  (math-working "arctan" sum)
  (let* ((nextx (math-mul-float x xnegsqr))
	 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
    (if (math-nearly-equal-float sum nextsum)
	sum
      (math-arctan-series nextsum (+ n 2) nextx xnegsqr))))

(defun calcFunc-arctan2 (y x)   ; [F R R] [Public]
  (if (Math-anglep y)
      (if (Math-anglep x)
	  (math-with-extra-prec 2

		  (setq y nil))))
	  (if (and y x)
	      (calcFunc-arctan2 y x)
	    '(var nan var-nan)))
      (calc-record-why 'anglep y)
      (list 'calcFunc-arctan2 y x))))

(defun math-arctan2-raw (y x)   ; [F R R]
  (cond ((math-zerop y)
	 (if (math-negp x) (math-pi)
	   (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0)))

	((math-posp y)
	 (math-add-float (math-arctan-raw (math-div-float y x))
			 (math-pi)))
	(t
	 (math-sub-float (math-arctan-raw (math-div-float y x))
			 (math-pi)))))

(defun calcFunc-arcsincos (x)   ; [V N] [Public]
  (if (and (Math-vectorp x)
	   (= (length x) 3))
      (calcFunc-arctan2 (nth 2 x) (nth 1 x))
    (math-reject-arg x "*Two-element vector expected")))



;;; Exponential function.

	((equal x '(neg (var inf var-inf)))
	 0)
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-exp x))))

(defun calcFunc-expm1 (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	((math-zerop x) '(float 0 0))
	(calc-symbolic-mode (signal 'inexact-result nil))

	((equal x '(neg (var inf var-inf)))
	 -1)
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-expm1 x))))

(defun calcFunc-exp10 (x)   ; [N N] [Public]
  (if (eq x 0)
      1
    (math-pow '(float 1 1) x)))

(defun math-exp-raw (x)   ; [N N]
  (cond ((math-zerop x) '(float 1 0))
	(calc-symbolic-mode (signal 'inexact-result nil))
	((eq (car x) 'cplx)

		(hfrac (math-sub-float x (math-mul-float (math-float hint)
							 '(float 5 -1)))))
	   (math-mul-float (math-ipow (math-sqrt-e) hint)
			   (math-add-float '(float 1 0)
					   (math-exp-minus-1-raw hfrac)))))
	(t (math-add-float '(float 1 0) (math-exp-minus-1-raw x)))))

(defun math-exp-minus-1-raw (x)   ; [F F]
  (math-exp-series x 2 3 x x))

(defun math-exp-series (sum nfac n xpow x)
  (math-working "exp" sum)
  (let* ((nextx (math-mul-float xpow x))
	 (nextsum (math-add-float sum (math-div-float nextx
						      (math-float nfac)))))
    (if (math-nearly-equal-float sum nextsum)
	sum
      (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x))))



;;; Logarithms.

	((math-infinitep x)
	 (if (equal x '(var nan var-nan))
	     x
	   '(var inf var-inf)))
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-ln x))))

(defun calcFunc-log10 (x)   ; [N N] [Public]
  (cond ((math-equal-int x 1)
	 (if (math-floatp x) '(float 0 0) 0))
	((and (Math-integerp x)

	((math-infinitep x)
	 (if (equal x '(var nan var-nan))
	     x
	   '(var inf var-inf)))
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-log10 x))))

(defun calcFunc-log (x &optional b)   ; [N N N] [Public]
  (cond ((or (null b) (equal b '(var e var-e)))
	 (math-normalize (list 'calcFunc-ln x)))
	((or (eq b 10) (equal b '(float 1 1)))

	     (math-infinitep b))
	 (math-div (calcFunc-ln x) (calcFunc-ln b)))
	(t (if (Math-numberp b)
	       (calc-record-why 'numberp x)
	     (calc-record-why 'numberp b))
	   (list 'calcFunc-log x b))))

(defun calcFunc-alog (x &optional b)
  (cond ((or (null b) (equal b '(var e var-e)))
	 (math-normalize (list 'calcFunc-exp x)))
	(t (math-pow b x))))

(defun calcFunc-ilog (x b)
  (if (and (math-natnump x) (not (eq x 0))
	   (math-natnump b) (not (eq b 0)))
      (if (eq b 1)
	  (math-reject-arg x "*Logarithm base one")
	(if (Math-natnum-lessp x b)
	    0
	  (cdr (math-integer-log x b))))
    (math-floor (calcFunc-log x b))))

(defun math-integer-log (x b)
  (let ((pows (list b))
	(pow (math-sqr b))
	next

      (setq n (/ n 2)
	    next (math-mul pow (car pows)))
      (or (Math-lessp x next)
	  (setq pow next
		sum (+ sum n))))
    (cons (equal pow x) sum)))


(defun math-log-base-raw (b)   ; [N N]
  (if (not (and (equal (car math-log-base-cache) b)
		(eq (nth 1 math-log-base-cache) calc-internal-prec)))
      (setq math-log-base-cache (list b calc-internal-prec
				      (math-ln-raw (math-float b)))))
  (nth 2 math-log-base-cache))
(setq math-log-base-cache nil)

(defun calcFunc-lnp1 (x)   ; [N N] [Public]
  (cond ((Math-equal-int x -1)
	 (if calc-infinite-mode

	((math-infinitep x)
	 (if (equal x '(var nan var-nan))
	     x
	   '(var inf var-inf)))
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-lnp1 x))))

(defun math-ln-raw (x)    ; [N N] --- must be float format!
  (cond ((eq (car-safe x) 'cplx)
	 (list 'cplx
	       (math-mul-float (math-ln-raw

	  (list 'cplx   ; negative and real
		(math-ln-raw (math-neg-float x))
		(math-pi))))
	(t (list 'cplx   ; negative and real
		 (math-ln-raw (math-neg-float x))
		 (math-pi)))))

(defun math-ln-raw-2 (x)    ; [F F]
  (cond ((math-lessp-float '(float 14 -1) x)
	 (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1)))
			 (math-ln-2)))
	(t    ; now .7 < x <= 1.4
	 (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
					(math-add-float x '(float 1 0)))))))

(defun math-ln-raw-3 (x)   ; [F F]
  (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
		  '(float 2 0)))

;;; Compute ln((1+x)/(1-x))
(defun math-ln-raw-series (sum n x xsqr)
  (math-working "log" sum)
  (let* ((nextx (math-mul-float x xsqr))
	 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
    (if (math-nearly-equal-float sum nextsum)
	sum
      (math-ln-raw-series nextsum (+ n 2) nextx xsqr))))

(defun math-ln-plus-1-raw (x)
  (math-lnp1-series x 2 x (math-neg x)))

(defun math-lnp1-series (sum n xpow x)
  (math-working "lnp1" sum)
  (let* ((nextx (math-mul-float xpow x))
	 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
    (if (math-nearly-equal-float sum nextsum)
	sum
      (math-lnp1-series nextsum (1+ n) nextx x))))

(math-defcache math-ln-10 (float (bigpos 018 684 045 994 092 585 302 2) -21)
  (math-ln-raw-2 '(float 1 1)))

(math-defcache math-ln-2 (float (bigpos 417 309 945 559 180 147 693) -21)

	((or (equal x '(var inf var-inf))
	     (equal x '(neg (var inf var-inf)))
	     (equal x '(var nan var-nan)))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-sinh x))))
(put 'calcFunc-sinh 'math-expandable t)

(defun calcFunc-cosh (x)   ; [N N] [Public]
  (cond ((eq x 0) 1)
	(math-expand-formulas

	((or (equal x '(var inf var-inf))
	     (equal x '(neg (var inf var-inf)))
	     (equal x '(var nan var-nan)))
	 (math-abs x))
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-cosh x))))
(put 'calcFunc-cosh 'math-expandable t)

(defun calcFunc-tanh (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	(math-expand-formulas

	((equal x '(neg (var inf var-inf)))
	 -1)
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-tanh x))))
(put 'calcFunc-tanh 'math-expandable t)

(defun calcFunc-arcsinh (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	(math-expand-formulas

	((or (equal x '(var inf var-inf))
	     (equal x '(neg (var inf var-inf)))
	     (equal x '(var nan var-nan)))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arcsinh x))))
(put 'calcFunc-arcsinh 'math-expandable t)

(defun calcFunc-arccosh (x)   ; [N N] [Public]
  (cond ((eq x 1) 0)
	((and (eq x -1) calc-symbolic-mode)

	((or (equal x '(var inf var-inf))
	     (equal x '(neg (var inf var-inf)))
	     (equal x '(var nan var-nan)))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arccosh x))))
(put 'calcFunc-arccosh 'math-expandable t)

(defun calcFunc-arctanh (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	((and (Math-equal-int x 1) calc-infinite-mode)

			 (calcFunc-arctanh (nth 2 x))
			 (calcFunc-arctanh (nth 3 x))))
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-arctanh x))))
(put 'calcFunc-arctanh 'math-expandable t)


;;; Convert A from HMS or degrees to radians.
(defun calcFunc-rad (a)   ; [R R] [Public]

	 (math-make-sdev (calcFunc-rad (nth 1 a))
			 (calcFunc-rad (nth 2 a))))
	(math-expand-formulas
	 (math-div (math-mul a '(var pi var-pi)) 180))
	((math-infinitep a) a)
	(t (list 'calcFunc-rad a))))
(put 'calcFunc-rad 'math-expandable t)

;;; Convert A from HMS or radians to degrees.
(defun calcFunc-deg (a)   ; [R R] [Public]
  (cond ((or (Math-numberp a)

	 (math-make-sdev (calcFunc-deg (nth 1 a))
			 (calcFunc-deg (nth 2 a))))
	(math-expand-formulas
	 (math-div (math-mul 180 a) '(var pi var-pi)))
	((math-infinitep a) a)
	(t (list 'calcFunc-deg a))))
(put 'calcFunc-deg 'math-expandable t)

;;; calc-math.el ends here