view lisp/calc/calc-math.el @ 110410:f2e111723c3a

Merge changes made in Gnus trunk. Reimplement nnimap, and do tweaks to the rest of the code to support that. * gnus-int.el (gnus-finish-retrieve-group-infos) (gnus-retrieve-group-data-early): New functions. * gnus-range.el (gnus-range-nconcat): New function. * gnus-start.el (gnus-get-unread-articles): Support early retrieval of data. (gnus-read-active-for-groups): Support finishing the early retrieval of data. * gnus-sum.el (gnus-summary-move-article): Pass the move-to group name if the move is internal, so that nnimap can do fast internal moves. * gnus.el (gnus-article-special-mark-lists): Add uid/active tuples, for nnimap usage. * nnimap.el: Rewritten. * nnmail.el (nnmail-inhibit-default-split-group): New internal variable to allow the mail splitting to not return a default group. This is useful for nnimap, which will leave unmatched mail in the inbox. * utf7.el (utf7-encode): Autoload. Implement shell connection. * nnimap.el (nnimap-open-shell-stream): New function. (nnimap-open-connection): Use it. Get the number of lines by using BODYSTRUCTURE. (nnimap-transform-headers): Get the number of lines in each message. (nnimap-retrieve-headers): Query for BODYSTRUCTURE so that we get the number of lines. Not all servers return UIDNEXT. Work past this problem. Remove junk from end of file. Fix typo in "bogus" section. Make capabilties be case-insensitive. Require cl when compiling. Don't bug out if the LIST command doesn't have any parameters. 2010-09-17 Knut Anders Hatlen <[e-mail]> (tiny change) * nnimap.el (nnimap-get-groups): Don't bug out if the LIST command doesn't have any parameters. (mm-text-html-renderer): Document gnus-article-html. 2010-09-17 Julien Danjou <[e-mail]> (tiny fix) * mm-decode.el (mm-text-html-renderer): Document gnus-article-html. * dgnushack.el: Define netrc-credentials. If the user doesn't have a /etc/services, supply some sensible port defaults. Have `unseen-or-unread' select an unread unseen article first. (nntp-open-server): Return whether the open was successful or not. Throughout all files, replace (save-excursion (set-buffer ...)) with (with-current-buffer ... ). Save result so that it doesn't say "failed" all the time. Add ~/.authinfo to the default, since that's probably most useful for users. Don't use the "finish" method when we're reading from the agent. Add some more nnimap-relevant agent stuff to nnagent.el. * nnimap.el (nnimap-with-process-buffer): Removed. Revert one line that was changed by mistake in the last checkin. (nnimap-open-connection): Don't error out when we can't make a connection nnimap-related changes to avoid bugging out if we can't contact a server. * gnus-start.el (gnus-get-unread-articles): Don't try to scan groups from methods that are denied. * nnimap.el (nnimap-possibly-change-group): Return nil if we can't log in. (nnimap-finish-retrieve-group-infos): Make sure we're not waiting for nothing. * gnus-sum.el (gnus-select-newsgroup): Indent.

;;; calc-math.el --- mathematical functions for Calc

;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
;;   2005, 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc.

;; Author: David Gillespie <[e-mail]>
;; Maintainer: Jay Belanger <[e-mail]>

;; This file is part of GNU Emacs.

;; GNU Emacs is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.

;; GNU Emacs is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.

;; You should have received a copy of the GNU General Public License
;; along with GNU Emacs.  If not, see <http://www.gnu.org/licenses/>.

;;; Commentary:

;;; Code:

;; This file is autoloaded from calc-ext.el.

(require 'calc-ext)
(require 'calc-macs)


;;; Find out how many 9s in 9.9999... will give distinct Emacs floats,
;;; then back off by one.

(defvar math-emacs-precision
  (let* ((n 1)
         (x 9)
         (xx (+ x (* 9 (expt 10 (- n))))))
    (while (/= x xx)
      (progn
        (setq n (1+ n))
        (setq x xx)
        (setq xx (+ x (* 9 (expt 10 (- n)))))))
    (1- n))
  "The number of digits in an Emacs float.")

;;; Find the largest power of 10 which is an Emacs float, 
;;; then back off by one so that any float d.dddd...eN 
;;; is an Emacs float, for acceptable d.dddd....

(defvar math-largest-emacs-expt
  (let ((x 1)
        (pow 1e2))
    ;; The following loop is for efficiency; it should stop when 
    ;; 10^(2x) is too large.  This could be indicated by a range 
    ;; error when computing 10^(2x) or an infinite value for 10^(2x).
    (while (and
            pow
            (< pow 1.0e+INF))
      (setq x (* 2 x))
      (setq pow (condition-case nil
                    (expt 10.0 (* 2 x))
                  (error nil))))
    ;; The following loop should stop when 10^(x+1) is too large.
    (setq pow (condition-case nil
                    (expt 10.0 (1+ x))
                  (error nil)))
    (while (and
            pow
            (< pow 1.0e+INF))
      (setq x (1+ x))
      (setq pow (condition-case nil
                    (expt 10.0 (1+ x))
                  (error nil))))
    (1- x))
  "The largest exponent which Calc will convert to an Emacs float.")

(defvar math-smallest-emacs-expt
  (let ((x -1))
    (while (condition-case nil
               (> (expt 10.0 x) 0.0)
             (error nil))
      (setq x (* 2 x)))
    (setq x (/ x 2))
    (while (condition-case nil
               (> (expt 10.0 x) 0.0)
             (error nil))
      (setq x (1- x)))
    (+ x 2))
    "The smallest exponent which Calc will convert to an Emacs float.")

(defun math-use-emacs-fn (fn x)
  "Use the native Emacs function FN to evaluate the Calc number X.
If this can't be done, return NIL."
  (and
   (<= calc-internal-prec math-emacs-precision)
   (math-realp x)
   (let* ((fx (math-float x))
          (xpon (+ (nth 2 x) (1- (math-numdigs (nth 1 x))))))
     (and (<= math-smallest-emacs-expt xpon)
          (<= xpon math-largest-emacs-expt)
          (condition-case nil
              (math-read-number
               (number-to-string
                (funcall fn 
			 (string-to-number 
			  (let 
                              ((calc-number-radix 10)
                               (calc-twos-complement-mode nil)
                               (calc-float-format (list 'float calc-internal-prec))
                               (calc-group-digits nil)
                               (calc-point-char "."))
			    (math-format-number (math-float x)))))))
            (error nil))))))

(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-is-hyperbolic)
	   (if calc-symbolic-mode
	       (calc-pop-push-record 0 "phi" '(var phi var-phi))
	     (calc-pop-push-record 0 "phi" (math-phi)))
	 (if calc-symbolic-mode
	     (calc-pop-push-record 0 "gmma" '(var gamma var-gamma))
	   (calc-pop-push-record 0 "gmma" (math-gamma-const))))
     (if (calc-is-hyperbolic)
	 (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-sec (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (calc-unary-op "sech" 'calcFunc-sech arg)
     (calc-unary-op "sec" 'calcFunc-sec arg))))

(defun calc-sech (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-sec 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-csc (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (calc-unary-op "csch" 'calcFunc-csch arg)
     (calc-unary-op "csc" 'calcFunc-csc arg))))

(defun calc-csch (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-csc 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-cot (arg)
  (interactive "P")
  (calc-slow-wrapper
   (if (calc-is-hyperbolic)
       (calc-unary-op "coth" 'calcFunc-coth arg)
     (calc-unary-op "cot" 'calcFunc-cot arg))))

(defun calc-coth (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-cot 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!
(defun math-isqrt (a)
  (cond ((Math-zerop a) a)
	((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)))
    (if (= (% len 2) 0)
	(let* ((top (nthcdr (- len 2) a)))
	  (math-isqrt-bignum-iter
	   a
	   (math-scale-bignum-digit-size
	    (math-bignum-big
	     (1+ (math-isqrt-small
		  (+ (* (nth 1 top) math-bignum-digit-size) (car top)))))
	    (1- (/ len 2)))))
      (let* ((top (nth (1- len) a)))
	(math-isqrt-bignum-iter
	 a
	 (math-scale-bignum-digit-size
	  (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))
	 (g2 (math-div2-bignum s))
	 (comp (math-compare-bignum g2 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-digit-size (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 1000000) 10000)
                 ((>= 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.
;;; [T N] if possible, else [F N] if possible, else [C N].  [Public]
(defun math-sqrt (a)
  (or
   (and (Math-zerop a) a)
   (and (math-known-nonposp a)
	(math-imaginary (math-sqrt (math-neg a))))
   (and (integerp a)
	(let ((sqrt (math-isqrt-small a)))
	  (if (= (* sqrt sqrt) a)
	      sqrt
	    (if calc-symbolic-mode
		(list 'calcFunc-sqrt a)
	      (math-sqrt-float (math-float a) (math-float sqrt))))))
   (and (eq (car-safe a) 'bigpos)
	(let* ((res (math-isqrt-bignum (cdr a)))
	       (sqrt (math-normalize (cons 'bigpos (cdr res)))))
	  (if (car res)
	      sqrt
	    (if calc-symbolic-mode
		(list 'calcFunc-sqrt a)
	      (math-sqrt-float (math-float a) (math-float sqrt))))))
   (and (eq (car-safe a) 'frac)
	(let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a)))))
	       (num-sqrt (math-normalize (cons 'bigpos (cdr num-res))))
	       (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a)))))
	       (den-sqrt (math-normalize (cons 'bigpos (cdr den-res)))))
	  (if (and (car num-res) (car den-res))
	      (list 'frac num-sqrt den-sqrt)
	    (if calc-symbolic-mode
		(if (or (car num-res) (car den-res))
		    (math-div (if (car num-res)
				  num-sqrt (list 'calcFunc-sqrt (nth 1 a)))
			      (if (car den-res)
				  den-sqrt (list 'calcFunc-sqrt (nth 2 a))))
		  (list 'calcFunc-sqrt a))
	      (math-sqrt-float (math-float a)
			       (math-div (math-float num-sqrt) den-sqrt))))))
   (and (eq (car-safe a) 'float)
	(if calc-symbolic-mode
	    (if (= (% (nth 2 a) 2) 0)
		(let ((res (math-isqrt-bignum
			    (cdr (Math-bignum-test (nth 1 a))))))
		  (if (car res)
		      (math-make-float (math-normalize
					(cons 'bigpos (cdr res)))
				       (/ (nth 2 a) 2))
		    (signal 'inexact-result nil)))
	      (signal 'inexact-result nil))
	  (math-sqrt-float a)))
   (and (eq (car-safe a) 'cplx)
	(math-with-extra-prec 2
	  (let* ((d (math-abs a))
		 (imag (math-sqrt (math-mul (math-sub d (nth 1 a))
					    '(float 5 -1)))))
	    (list 'cplx
		  (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1)))
		  (if (math-negp (nth 2 a)) (math-neg imag) imag)))))
   (and (eq (car-safe a) 'polar)
	(list 'polar
	      (math-sqrt (nth 1 a))
	      (math-mul (nth 2 a) '(float 5 -1))))
   (and (eq (car-safe a) 'sdev)
	(let ((sqrt (math-sqrt (nth 1 a))))
	  (math-make-sdev sqrt
			  (math-div (nth 2 a) (math-mul sqrt 2)))))
   (and (eq (car-safe a) 'intv)
	(not (math-negp (nth 2 a)))
	(math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a))))
   (and (eq (car-safe a) '*)
	(or (math-known-nonnegp (nth 1 a))
	    (math-known-nonnegp (nth 2 a)))
	(math-mul (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
   (and (eq (car-safe a) '/)
	(or (and (math-known-nonnegp (nth 2 a))
		 (math-div (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
	    (and (math-known-nonnegp (nth 1 a))
		 (not (math-equal-int (nth 1 a) 1))
		 (math-mul (math-sqrt (nth 1 a))
			   (math-sqrt (math-div 1 (nth 2 a)))))))
   (and (eq (car-safe a) '^)
	(math-known-evenp (nth 2 a))
	(math-known-realp (nth 1 a))
	(math-abs (math-pow (nth 1 a) (math-div (nth 2 a) 2))))
   (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)
    (cond
     ((math-use-emacs-fn 'sqrt a))
     (t
      (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)) (nth 1 b)))
	   (or (eq ediff 0)
	       (and (not (consp ediff))
		    (< ediff 10)
		    (> ediff -10)
		    (= (math-numdigs (nth 1 a)) calc-internal-prec))))
	  ((= ediff 1)
	   (setq ediff (math-add (Math-integer-neg (nth 1 b))
				 (math-scale-int (nth 1 a) 1)))
	   (and (not (consp ediff))
		(< ediff 10)
		(> ediff -10)
		(= (math-numdigs (nth 1 b)) calc-internal-prec)))
	  ((= ediff -1)
	   (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)))
  (if (eq (car b) 'polar) (setq b (math-complex b)))
  (if (eq (car a) 'cplx)
      (if (eq (car b) 'cplx)
	  (and (or (math-nearly-equal-float (nth 1 a) (nth 1 b))
		   (and (math-nearly-zerop-float (nth 1 a) (nth 2 a))
			(math-nearly-zerop-float (nth 1 b) (nth 2 b))))
	       (or (math-nearly-equal-float (nth 2 a) (nth 2 b))
		   (and (math-nearly-zerop-float (nth 2 a) (nth 1 a))
			(math-nearly-zerop-float (nth 2 b) (nth 1 b)))))
	(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))
	((not (Math-scalarp a))
	 (if (math-infinitep a)
	     (if (math-infinitep b)
		 (if (equal a b)
		     a
		   '(var nan var-nan))
	       a)
	   (calc-record-why 'scalarp a)
	   (list 'calcFunc-hypot a b)))
	((not (Math-scalarp b))
	 (if (math-infinitep b)
	     b
	   (calc-record-why 'scalarp b)
	   (list 'calcFunc-hypot a b)))
	((and (Math-numberp a) (Math-numberp b))
	 (math-with-extra-prec 1
	   (math-sqrt (math-add (calcFunc-abssqr a) (calcFunc-abssqr b)))))
	((eq (car-safe a) 'hms)
	 (if (eq (car-safe b) 'hms)   ; this helps sdev's of hms forms
	     (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))
	((Math-zerop a) a)
	((Math-negp a) nil)
	((Math-integerp a)
	 (let ((root (math-nth-root-integer a n)))
	   (if (car root)
	       (cdr root)
	     (and (not calc-symbolic-mode)
		  (math-nth-root-float (math-float a) n
				       (math-float (cdr root)))))))
	((eq (car-safe a) 'frac)
	 (let* ((num-root (math-nth-root-integer (nth 1 a) n))
		(den-root (math-nth-root-integer (nth 2 a) n)))
	   (if (and (car num-root) (car den-root))
	       (list 'frac (cdr num-root) (cdr den-root))
	     (and (not calc-symbolic-mode)
		  (math-nth-root-float
		   (math-float a) n
		   (math-div-float (math-float (cdr num-root))
				   (math-float (cdr den-root))))))))
	((eq (car-safe a) 'float)
	 (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)))

;; The variables math-nrf-n, math-nrf-nf and math-nrf-nfm1 are local
;; to math-nth-root-float, but are used by math-nth-root-float-iter,
;; which is called by math-nth-root-float.
(defvar math-nrf-n)
(defvar math-nrf-nf)
(defvar math-nrf-nfm1)

(defun math-nth-root-float (a math-nrf-n &optional guess)
  (math-inexact-result)
  (math-with-extra-prec 1
    (let ((math-nrf-nf (math-float math-nrf-n))
	  (math-nrf-nfm1 (math-float (1- math-nrf-n))))
      (math-nth-root-float-iter a (or guess
				      (math-make-float
				       1 (/ (+ (math-numdigs (nth 1 a))
					       (nth 2 a)
					       (/ math-nrf-n 2))
					    math-nrf-n)))))))

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

;; The variable math-nri-n is local to math-nth-root-integer, but
;; is used by math-nth-root-int-iter, which is called by
;; math-nth-root-int.
(defvar math-nri-n)

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

(defun math-nth-root-int-iter (a guess)
  (math-working "root" guess)
  (let* ((q (math-idivmod a (math-ipow guess (1- math-nri-n))))
	 (s (math-add (car q) (math-mul (1- math-nri-n) guess)))
	 (g2 (math-idivmod s math-nri-n)))
    (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.

;;; All of these functions are defined on the complex plane.
;;; (Branch cuts, etc. follow Steele's Common Lisp book.)

;;; 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).

(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.

(defun calcFunc-sin (x)   ; [N N] [Public]
  (cond ((and (integerp x)
	      (if (eq calc-angle-mode 'deg)
		  (= (% x 90) 0)
		(= x 0)))
	 (aref [0 1 0 -1] (math-mod (/ x 90) 4)))
	((Math-scalarp x)
	 (math-with-extra-prec 2
	   (math-sin-raw (math-to-radians (math-float x)))))
	((eq (car x) 'sdev)
	 (if (math-constp x)
	     (math-with-extra-prec 2
	       (let* ((xx (math-to-radians (math-float (nth 1 x))))
		      (xs (math-to-radians (math-float (nth 2 x))))
		      (sc (math-sin-cos-raw xx)))
		 (math-make-sdev (car sc) (math-mul xs (cdr sc)))))
	   (math-make-sdev (calcFunc-sin (nth 1 x))
			   (math-mul (nth 2 x) (calcFunc-cos (nth 1 x))))))
	((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)
		(= x 0)))
	 (aref [1 0 -1 0] (math-mod (/ x 90) 4)))
	((Math-scalarp x)
	 (math-with-extra-prec 2
	   (math-cos-raw (math-to-radians (math-float x)))))
	((eq (car x) 'sdev)
	 (if (math-constp x)
	     (math-with-extra-prec 2
	       (let* ((xx (math-to-radians (math-float (nth 1 x))))
		      (xs (math-to-radians (math-float (nth 2 x))))
		      (sc (math-sin-cos-raw xx)))
		 (math-make-sdev (cdr sc) (math-mul xs (car sc)))))
	   (math-make-sdev (calcFunc-cos (nth 1 x))
			   (math-mul (nth 2 x) (calcFunc-sin (nth 1 x))))))
	((and (eq (car x) 'intv) (math-intv-constp x))
	 (math-with-extra-prec 2
	   (let* ((xx (math-to-radians (math-float x)))
		  (na (math-floor (math-div (nth 2 xx) (math-pi))))
		  (nb (math-floor (math-div (nth 3 xx) (math-pi))))
		  (span (math-sub nb na)))
	     (if (memq span '(0 1))
		 (let ((int (math-sort-intv (nth 1 x)
					    (math-cos-raw (nth 2 xx))
					    (math-cos-raw (nth 3 xx)))))
		   (if (eq span 1)
		       (if (math-evenp na)
			   (math-make-intv (logior (nth 1 x) 2)
					   -1
					   (nth 3 int))
			 (math-make-intv (logior (nth 1 x) 1)
					 (nth 2 int)
					 1))
		     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)
		(= x 0)))
	 0)
	((Math-scalarp x)
	 (math-with-extra-prec 2
	   (math-tan-raw (math-to-radians (math-float x)))))
	((eq (car x) 'sdev)
	 (if (math-constp x)
	     (math-with-extra-prec 2
	       (let* ((xx (math-to-radians (math-float (nth 1 x))))
		      (xs (math-to-radians (math-float (nth 2 x))))
		      (sc (math-sin-cos-raw xx)))
		 (if (and (math-zerop (cdr sc)) (not calc-infinite-mode))
		     (progn
		       (calc-record-why "*Division by zero")
		       (list 'calcFunc-tan x))
		   (math-make-sdev (math-div-float (car sc) (cdr sc))
				   (math-div-float xs (math-sqr (cdr sc)))))))
	   (math-make-sdev (calcFunc-tan (nth 1 x))
			   (math-div (nth 2 x)
				     (math-sqr (calcFunc-cos (nth 1 x)))))))
	((and (eq (car x) 'intv) (math-intv-constp x))
	 (or (math-with-extra-prec 2
	       (let* ((xx (math-to-radians (math-float x)))
		      (na (math-floor (math-div (math-sub (nth 2 xx)
							  (math-pi-over-2))
						(math-pi))))
		      (nb (math-floor (math-div (math-sub (nth 3 xx)
							  (math-pi-over-2))
						(math-pi)))))
		 (and (equal na nb)
		      (math-sort-intv (nth 1 x)
				      (math-tan-raw (nth 2 xx))
				      (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 calcFunc-sec (x)
  (cond ((and (integerp x)
              (eq calc-angle-mode 'deg)
              (= (% x 180) 0))
         (if (= (% x 360) 0)
             1
           -1))
        ((and (integerp x)
              (eq calc-angle-mode 'rad)
              (= x 0))
         1)
        ((Math-scalarp x)
         (math-with-extra-prec 2
           (math-sec-raw (math-to-radians (math-float x)))))
        ((eq (car x) 'sdev)
         (if (math-constp x)
             (math-with-extra-prec 2
               (let* ((xx (math-to-radians (math-float (nth 1 x))))
                      (xs (math-to-radians (math-float (nth 2 x))))
                      (sc (math-sin-cos-raw xx)))
                 (if (and (math-zerop (cdr sc))
                          (not calc-infinite-mode))
                     (progn
                       (calc-record-why "*Division by zero")
                       (list 'calcFunc-sec x))
                   (math-make-sdev (math-div-float '(float 1 0) (cdr sc))
                                   (math-div-float
                                    (math-mul xs (car sc))
                                    (math-sqr (cdr sc)))))))
           (math-make-sdev (calcFunc-sec (nth 1 x))
                           (math-div 
                            (math-mul (nth 2 x)
                                      (calcFunc-sin (nth 1 x)))
                            (math-sqr (calcFunc-cos (nth 1 x)))))))
        ((and (eq (car x) 'intv)
              (math-intv-constp x))
         (math-with-extra-prec 2
           (let* ((xx (math-to-radians (math-float x)))
                  (na (math-floor (math-div (math-sub (nth 2 xx)
                                                      (math-pi-over-2))
                                            (math-pi))))
                  (nb (math-floor (math-div (math-sub (nth 3 xx)
                                                      (math-pi-over-2))
                                            (math-pi))))
                  (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2))))
                  (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2))))
                  (span (math-sub nbb naa)))
             (if (not (equal na nb))
                 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
               (let ((int (math-sort-intv (nth 1 x)
                                          (math-sec-raw (nth 2 xx))
                                          (math-sec-raw (nth 3 xx)))))
                 (if (eq span 1)
                     (if (math-evenp (math-div (math-add naa 1) 2))
                         (math-make-intv (logior (nth 1 int) 2)
                                         1
                                         (nth 3 int))
                       (math-make-intv (logior (nth 1 int) 1)
                                       (nth 2 int)
                                       -1))
                   int))))))
        ((equal x '(var nan var-nan))
         x)
        (t (calc-record-why 'scalarp x)
           (list 'calcFunc-sec x))))

(defun calcFunc-csc (x)
  (cond ((and (integerp x)
              (eq calc-angle-mode 'deg)
              (= (% (- x 90) 180) 0))
         (if (= (% (- x 90) 360) 0)
             1
           -1))
        ((Math-scalarp x)
         (math-with-extra-prec 2
           (math-csc-raw (math-to-radians (math-float x)))))
        ((eq (car x) 'sdev)
         (if (math-constp x)
             (math-with-extra-prec 2
               (let* ((xx (math-to-radians (math-float (nth 1 x))))
                      (xs (math-to-radians (math-float (nth 2 x))))
                      (sc (math-sin-cos-raw xx)))
                 (if (and (math-zerop (car sc))
                          (not calc-infinite-mode))
                     (progn
                       (calc-record-why "*Division by zero")
                       (list 'calcFunc-csc x))
                   (math-make-sdev (math-div-float '(float 1 0) (car sc))
                                   (math-div-float
                                    (math-mul xs (cdr sc))
                                    (math-sqr (car sc)))))))
           (math-make-sdev (calcFunc-csc (nth 1 x))
                           (math-div 
                            (math-mul (nth 2 x)
                                      (calcFunc-cos (nth 1 x)))
                            (math-sqr (calcFunc-sin (nth 1 x)))))))
        ((and (eq (car x) 'intv)
              (math-intv-constp x))
         (math-with-extra-prec 2
           (let* ((xx (math-to-radians (math-float x)))
                  (na (math-floor (math-div (nth 2 xx) (math-pi))))
                  (nb (math-floor (math-div (nth 3 xx) (math-pi))))
                  (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2))))
                  (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2))))
                  (span (math-sub nbb naa)))
             (if (not (equal na nb))
                 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
               (let ((int (math-sort-intv (nth 1 x)
                                          (math-csc-raw (nth 2 xx))
                                          (math-csc-raw (nth 3 xx)))))
                 (if (eq span 1)
                     (if (math-evenp (math-div naa 2))
                         (math-make-intv (logior (nth 1 int) 2)
                                         1
                                         (nth 3 int))
                       (math-make-intv (logior (nth 1 int) 1)
                                       (nth 2 int)
                                       -1))
                   int))))))
        ((equal x '(var nan var-nan))
         x)
        (t (calc-record-why 'scalarp x)
           (list 'calcFunc-csc x))))

(defun calcFunc-cot (x)   ; [N N] [Public]
  (cond ((and (integerp x)
	      (if (eq calc-angle-mode 'deg)
		  (= (% (- x 90) 180) 0)
		(= x 0)))
	 0)
	((Math-scalarp x)
	 (math-with-extra-prec 2
	   (math-cot-raw (math-to-radians (math-float x)))))
	((eq (car x) 'sdev)
	 (if (math-constp x)
	     (math-with-extra-prec 2
	       (let* ((xx (math-to-radians (math-float (nth 1 x))))
		      (xs (math-to-radians (math-float (nth 2 x))))
		      (sc (math-sin-cos-raw xx)))
		 (if (and (math-zerop (car sc)) (not calc-infinite-mode))
		     (progn
		       (calc-record-why "*Division by zero")
		       (list 'calcFunc-cot x))
		   (math-make-sdev (math-div-float (cdr sc) (car sc))
				   (math-div-float xs (math-sqr (car sc)))))))
	   (math-make-sdev (calcFunc-cot (nth 1 x))
			   (math-div (nth 2 x)
				     (math-sqr (calcFunc-sin (nth 1 x)))))))
	((and (eq (car x) 'intv) (math-intv-constp x))
	 (or (math-with-extra-prec 2
	       (let* ((xx (math-to-radians (math-float x)))
		      (na (math-floor (math-div (nth 2 xx) (math-pi))))
		      (nb (math-floor (math-div (nth 3 xx) (math-pi)))))
		 (and (equal na nb)
		      (math-sort-intv (nth 1 x)
				      (math-cot-raw (nth 2 xx))
				      (math-cot-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-cot x))))

(defun math-sin-raw (x &optional orgx)   ; [N N]
  (cond ((eq (car x) 'cplx)
	 (let* ((expx (math-exp-raw (nth 2 x)))
		(expmx (math-div-float '(float 1 0) expx))
		(sc (math-sin-cos-raw (nth 1 x))))
	   (list 'cplx
		 (math-mul-float (car sc)
				 (math-mul-float (math-add-float expx expmx)
						 '(float 5 -1)))
		 (math-mul-float (cdr sc)
				 (math-mul-float (math-sub-float expx expmx)
						 '(float 5 -1))))))
	((eq (car x) 'polar)
	 (math-polar (math-sin-raw (math-complex x))))
	((Math-integer-negp (nth 1 x))
	 (math-neg-float (math-sin-raw (math-neg-float x) (if orgx orgx x))))
	((math-lessp-float '(float 7 0) x)  ; avoid inf loops due to roundoff
	 (math-sin-raw (math-mod x (math-two-pi)) (if orgx orgx x)))
	(t (math-sin-raw-2 x (if orgx orgx 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) x)))

(defun math-sec-raw (x)   ; [N N]
  (cond ((eq (car x) 'cplx)
	 (let* ((x (math-mul x '(float 1 0)))
                (expx (math-exp-raw (nth 2 x)))
		(expmx (math-div-float '(float 1 0) expx))
                (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
                (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
		(sc (math-sin-cos-raw (nth 1 x)))
		(d (math-add-float 
                    (math-mul-float (math-sqr (car sc))
                                    (math-sqr sh))
                    (math-mul-float (math-sqr (cdr sc))
                                    (math-sqr ch)))))
	   (and (not (eq (nth 1 d) 0))
		(list 'cplx
		      (math-div-float (math-mul-float (cdr sc) ch) d)
		      (math-div-float (math-mul-float (car sc) sh) d)))))
	((eq (car x) 'polar)
	 (math-polar (math-sec-raw (math-complex x))))
	(t
	 (let ((cs (math-cos-raw x)))
           (if (eq cs 0)
               (math-div 1 0)
	     (math-div-float '(float 1 0) cs))))))

(defun math-csc-raw (x)   ; [N N]
  (cond ((eq (car x) 'cplx)
	 (let* ((x (math-mul x '(float 1 0)))
                (expx (math-exp-raw (nth 2 x)))
		(expmx (math-div-float '(float 1 0) expx))
                (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
                (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
		(sc (math-sin-cos-raw (nth 1 x)))
		(d (math-add-float 
                    (math-mul-float (math-sqr (car sc))
                                    (math-sqr ch))
                    (math-mul-float (math-sqr (cdr sc))
                                    (math-sqr sh)))))
	   (and (not (eq (nth 1 d) 0))
		(list 'cplx
		      (math-div-float (math-mul-float (car sc) ch) d)
		      (math-div-float (math-mul-float (cdr sc) sh) d)))))
	((eq (car x) 'polar)
	 (math-polar (math-csc-raw (math-complex x))))
	(t
	 (let ((sn (math-sin-raw x)))
           (if (eq sn 0)
               (math-div 1 0)
	     (math-div-float '(float 1 0) sn))))))

(defun math-cot-raw (x)   ; [N N]
  (cond ((eq (car x) 'cplx)
	 (let* ((x (math-mul x '(float 1 0)))
                (expx (math-exp-raw (nth 2 x)))
		(expmx (math-div-float '(float 1 0) expx))
                (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
                (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
		(sc (math-sin-cos-raw (nth 1 x)))
		(d (math-add-float 
                    (math-sqr (car sc))
                    (math-sqr sh))))
	   (and (not (eq (nth 1 d) 0))
		(list 'cplx
		      (math-div-float 
                       (math-mul-float (car sc) (cdr sc))
                       d)
                      (math-neg
                       (math-div-float 
                        (math-mul-float sh ch) 
                        d))))))
	((eq (car x) 'polar)
	 (math-polar (math-cot-raw (math-complex x))))
	(t
	 (let ((sc (math-sin-cos-raw x)))
	   (if (eq (nth 1 (car sc)) 0)
	       (math-div (cdr sc) 0)
	     (math-div-float (cdr sc) (car sc)))))))


;;; 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)))
		(expmx (math-div-float '(float 1 0) expx))
		(sc (math-sin-cos-raw (nth 1 x)))
		(d (math-add-float (cdr sc)
				   (math-mul-float (math-add-float expx expmx)
						   '(float 5 -1)))))
	   (and (not (eq (nth 1 d) 0))
		(list 'cplx
		      (math-div-float (car sc) d)
		      (math-div-float (math-mul-float (math-sub-float expx
								      expmx)
						      '(float 5 -1)) d)))))
	((eq (car x) 'polar)
	 (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))
					   orgx)))
	  ((math-lessp-float (math-pi-over-4) x)
	   (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-with-extra-prec -1 (math-nearly-zerop-float x orgx)) 
           '(float 0 0))
          ((math-use-emacs-fn 'sin x))
	  (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-with-extra-prec -1 (math-nearly-zerop-float x orgx))
         '(float 1 0))
        ((math-use-emacs-fn 'cos x))
	(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]
  (cond ((eq x 0) 0)
	((and (eq x 1) (eq calc-angle-mode 'deg)) 90)
	((and (eq x -1) (eq calc-angle-mode 'deg)) -90)
	(calc-symbolic-mode (signal 'inexact-result nil))
	((Math-numberp x)
	 (math-with-extra-prec 2
	   (math-from-radians (math-arcsin-raw (math-float x)))))
	((eq (car x) 'sdev)
	 (math-make-sdev (calcFunc-arcsin (nth 1 x))
			 (math-from-radians
			  (math-div (nth 2 x)
				    (math-sqrt
				     (math-sub 1 (math-sqr (nth 1 x))))))))
	((eq (car x) 'intv)
	 (math-sort-intv (nth 1 x)
			 (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)
	(calc-symbolic-mode (signal 'inexact-result nil))
	((Math-numberp x)
	 (math-with-extra-prec 2
	   (math-from-radians (math-arccos-raw (math-float x)))))
	((eq (car x) 'sdev)
	 (math-make-sdev (calcFunc-arccos (nth 1 x))
			 (math-from-radians
			  (math-div (nth 2 x)
				    (math-sqrt
				     (math-sub 1 (math-sqr (nth 1 x))))))))
	((eq (car x) 'intv)
	 (math-sort-intv (nth 1 x)
			 (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)
	((Math-numberp x)
	 (math-with-extra-prec 2
	   (math-from-radians (math-arctan-raw (math-float x)))))
	((eq (car x) 'sdev)
	 (math-make-sdev (calcFunc-arctan (nth 1 x))
			 (math-from-radians
			  (math-div (nth 2 x)
				    (math-add 1 (math-sqr (nth 1 x)))))))
	((eq (car x) 'intv)
	 (math-sort-intv (nth 1 x)
			 (calcFunc-arctan (nth 2 x))
			 (calcFunc-arctan (nth 3 x))))
	((equal x '(var inf var-inf))
	 (math-quarter-circle t))
	((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-ln-raw (math-add 1 (math-mul '(cplx 0 1) x)))
		      (math-ln-raw (math-add 1 (math-mul '(cplx 0 -1) x))))
		     '(cplx 0 2))))
	((Math-integer-negp (nth 1 x))
	 (math-neg-float (math-arctan-raw (math-neg-float x))))
	((math-zerop x) x)
        ((math-use-emacs-fn 'atan x))
	(calc-symbolic-mode (signal 'inexact-result nil))
	((math-equal-int x 1) (math-pi-over-4))
	((math-equal-int x -1) (math-neg (math-pi-over-4)))
	((math-lessp-float '(float 414214 -6) x)  ; if x > sqrt(2) - 1, reduce
	 (if (math-lessp-float '(float 1 0) x)
	     (math-sub-float (math-mul-float (math-pi) '(float 5 -1))
			     (math-arctan-raw (math-div-float '(float 1 0) x)))
	   (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
	    (math-from-radians (math-arctan2-raw (math-float y)
						 (math-float x))))
	(calc-record-why 'anglep x)
	(list 'calcFunc-arctan2 y x))
    (if (and (or (math-infinitep x) (math-anglep x))
	     (or (math-infinitep y) (math-anglep y)))
	(progn
	  (if (math-posp x)
	      (setq x 1)
	    (if (math-negp x)
		(setq x -1)
	      (or (math-zerop x)
		  (setq x nil))))
	  (if (math-posp y)
	      (setq y 1)
	    (if (math-negp y)
		(setq y -1)
	      (or (math-zerop y)
		  (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-zerop x)
	 (if (math-posp y)
	     (math-pi-over-2)
	   (math-neg (math-pi-over-2))))
	((math-posp x)
	 (math-arctan-raw (math-div-float y x)))
	((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.

(defun calcFunc-exp (x)   ; [N N] [Public]
  (cond ((eq x 0) 1)
	((and (memq x '(1 -1)) calc-symbolic-mode)
	 (if (eq x 1) '(var e var-e) (math-div 1 '(var e var-e))))
	((Math-numberp x)
	 (math-with-extra-prec 2 (math-exp-raw (math-float x))))
	((eq (car-safe x) 'sdev)
	 (let ((ex (calcFunc-exp (nth 1 x))))
	   (math-make-sdev ex (math-mul (nth 2 x) ex))))
	((eq (car-safe x) 'intv)
	 (math-make-intv (nth 1 x) (calcFunc-exp (nth 2 x))
			 (calcFunc-exp (nth 3 x))))
	((equal x '(var inf var-inf))
	 x)
	((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))
	((Math-numberp x)
	 (math-with-extra-prec 2
	   (let ((x (math-float x)))
	     (if (and (eq (car x) 'float)
		      (math-lessp-float x '(float 1 0))
		      (math-lessp-float '(float -1 0) x))
		 (math-exp-minus-1-raw x)
	       (math-add (math-exp-raw x) -1)))))
	((eq (car-safe x) 'sdev)
	 (if (math-constp x)
	     (let ((ex (calcFunc-expm1 (nth 1 x))))
	       (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1))))
	   (math-make-sdev (calcFunc-expm1 (nth 1 x))
			   (math-mul (nth 2 x) (calcFunc-exp (nth 1 x))))))
	((eq (car-safe x) 'intv)
	 (math-make-intv (nth 1 x)
			 (calcFunc-expm1 (nth 2 x))
			 (calcFunc-expm1 (nth 3 x))))
	((equal x '(var inf var-inf))
	 x)
	((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)
	 (let ((expx (math-exp-raw (nth 1 x)))
	       (sc (math-sin-cos-raw (nth 2 x))))
	   (list 'cplx
		 (math-mul-float expx (cdr sc))
		 (math-mul-float expx (car sc)))))
	((eq (car x) 'polar)
	 (let ((xc (math-complex x)))
	   (list 'polar
		 (math-exp-raw (nth 1 xc))
		 (math-from-radians (nth 2 xc)))))
        ((math-use-emacs-fn 'exp x))
	((or (math-lessp-float '(float 5 -1) x)
	     (math-lessp-float x '(float -5 -1)))
	 (if (math-lessp-float '(float 921035 1) x)
	     (math-overflow)
	   (if (math-lessp-float x '(float -921035 1))
	       (math-underflow)))
	 (let* ((two-x (math-mul-float x '(float 2 0)))
		(hint (math-scale-int (nth 1 two-x) (nth 2 two-x)))
		(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.

(defun calcFunc-ln (x)   ; [N N] [Public]
  (cond ((math-zerop x)
	 (if calc-infinite-mode
	     '(neg (var inf var-inf))
	   (math-reject-arg x "*Logarithm of zero")))
	((eq x 1) 0)
	((Math-numberp x)
	 (math-with-extra-prec 2 (math-ln-raw (math-float x))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-ln (nth 1 x))
			 (math-div (nth 2 x) (nth 1 x))))
	((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
					  (Math-zerop (nth 2 x))
					  (not (math-intv-constp x))))
	 (let ((calc-infinite-mode t))
	   (math-make-intv (nth 1 x) (calcFunc-ln (nth 2 x))
			   (calcFunc-ln (nth 3 x)))))
	((equal x '(var e var-e))
	 1)
	((and (eq (car-safe x) '^)
	      (equal (nth 1 x) '(var e var-e))
	      (math-known-realp (nth 2 x)))
	 (nth 2 x))
	((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-posp x)
	      (let ((res (math-integer-log x 10)))
		(and (car res)
		     (setq x (cdr res)))))
	 x)
	((and (eq (car-safe x) 'frac)
	      (eq (nth 1 x) 1)
	      (let ((res (math-integer-log (nth 2 x) 10)))
		(and (car res)
		     (setq x (- (cdr res))))))
	 x)
	((math-zerop x)
	 (if calc-infinite-mode
	     '(neg (var inf var-inf))
	   (math-reject-arg x "*Logarithm of zero")))
	(calc-symbolic-mode (signal 'inexact-result nil))
	((Math-numberp x)
	 (math-with-extra-prec 2
	   (let ((xf (math-float x)))
	     (if (eq (nth 1 xf) 0)
		 (math-reject-arg x "*Logarithm of zero"))
	     (if (Math-integer-posp (nth 1 xf))
		 (if (eq (nth 1 xf) 1)    ; log10(1*10^n) = n
		     (math-float (nth 2 xf))
		   (let ((xdigs (1- (math-numdigs (nth 1 xf)))))
		     (math-add-float
		      (math-div-float (math-ln-raw-2
				       (list 'float (nth 1 xf) (- xdigs)))
				      (math-ln-10))
		      (math-float (+ (nth 2 xf) xdigs)))))
	       (math-div (calcFunc-ln xf) (math-ln-10))))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-log10 (nth 1 x))
			 (math-div (nth 2 x)
				   (math-mul (nth 1 x) (math-ln-10)))))
	((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
					  (not (math-intv-constp x))))
	 (math-make-intv (nth 1 x)
			 (calcFunc-log10 (nth 2 x))
			 (calcFunc-log10 (nth 3 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-normalize (list 'calcFunc-log10 x)))
	((math-zerop x)
	 (if calc-infinite-mode
	     (math-div (calcFunc-ln x) (calcFunc-ln b))
	   (math-reject-arg x "*Logarithm of zero")))
	((math-zerop b)
	 (if calc-infinite-mode
	     (math-div (calcFunc-ln x) (calcFunc-ln b))
	   (math-reject-arg b "*Logarithm of zero")))
	((math-equal-int b 1)
	 (if calc-infinite-mode
	     (math-div (calcFunc-ln x) 0)
	   (math-reject-arg b "*Logarithm base one")))
	((math-equal-int x 1)
	 (if (math-floatp b) '(float 0 0) 0))
	((and (Math-ratp x) (Math-ratp b)
	      (math-posp x) (math-posp b)
	      (let* ((sign 1) (inv nil)
		     (xx (if (Math-lessp 1 x)
			     x
			   (setq sign -1)
			   (math-div 1 x)))
		     (bb (if (Math-lessp 1 b)
			     b
			   (setq sign (- sign))
			   (math-div 1 b)))
		     (res (if (Math-lessp xx bb)
			      (setq inv (math-integer-log bb xx))
			    (math-integer-log xx bb))))
		(and (car res)
		     (setq x (if inv
				 (math-div 1 (* sign (cdr res)))
			       (* sign (cdr res)))))))
	 x)
	(calc-symbolic-mode (signal 'inexact-result nil))
	((and (Math-numberp x) (Math-numberp b))
	 (math-with-extra-prec 2
	   (math-div (math-ln-raw (math-float x))
		     (math-log-base-raw b))))
	((and (eq (car-safe x) 'sdev)
	      (Math-numberp b))
	 (math-make-sdev (calcFunc-log (nth 1 x) b)
			 (math-div (nth 2 x)
				   (math-mul (nth 1 x)
					     (math-log-base-raw b)))))
	((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
					  (not (math-intv-constp x)))
	      (math-realp b))
	 (math-make-intv (nth 1 x)
			 (calcFunc-log (nth 2 x) b)
			 (calcFunc-log (nth 3 x) b)))
	((or (eq (car-safe x) 'intv) (eq (car-safe b) 'intv))
	 (math-div (calcFunc-ln x) (calcFunc-ln b)))
	((or (math-infinitep x)
	     (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
	sum n)
    (while (not (Math-lessp x pow))
      (setq pows (cons pow pows)
	    pow (math-sqr pow)))
    (setq n (lsh 1 (1- (length pows)))
	  sum n
	  pow (car pows))
    (while (and (setq pows (cdr pows))
		(Math-lessp pow x))
      (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)))


(defvar math-log-base-cache nil)
(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))

(defun calcFunc-lnp1 (x)   ; [N N] [Public]
  (cond ((Math-equal-int x -1)
	 (if calc-infinite-mode
	     '(neg (var inf var-inf))
	   (math-reject-arg x "*Logarithm of zero")))
	((eq x 0) 0)
	((math-zerop x) '(float 0 0))
	(calc-symbolic-mode (signal 'inexact-result nil))
	((Math-numberp x)
	 (math-with-extra-prec 2
	   (let ((x (math-float x)))
	     (if (and (eq (car x) 'float)
		      (math-lessp-float x '(float 5 -1))
		      (math-lessp-float '(float -5 -1) x))
		 (math-ln-plus-1-raw x)
	       (math-ln-raw (math-add-float x '(float 1 0)))))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-lnp1 (nth 1 x))
			 (math-div (nth 2 x) (math-add (nth 1 x) 1))))
	((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
					  (not (math-intv-constp x))))
	 (math-make-intv (nth 1 x)
			 (calcFunc-lnp1 (nth 2 x))
			 (calcFunc-lnp1 (nth 3 x))))
	((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
				(math-add-float (math-sqr-float (nth 1 x))
						(math-sqr-float (nth 2 x))))
			       '(float 5 -1))
	       (math-arctan2-raw (nth 2 x) (nth 1 x))))
	((eq (car x) 'polar)
	 (math-polar (list 'cplx
			   (math-ln-raw (nth 1 x))
			   (math-to-radians (nth 2 x)))))
	((Math-equal-int x 1)
	 '(float 0 0))
	(calc-symbolic-mode (signal 'inexact-result nil))
	((math-posp (nth 1 x))    ; positive and real
         (cond 
          ((math-use-emacs-fn 'log x))
          (t
           (let ((xdigs (1- (math-numdigs (nth 1 x)))))
             (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs)))
                             (math-mul-float (math-float (+ (nth 2 x) xdigs))
                                             (math-ln-10)))))))
	((math-zerop x)
	 (math-reject-arg x "*Logarithm of zero"))
	((eq calc-complex-mode 'polar)    ; negative and real
	 (math-polar
	  (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))))

(defconst math-approx-ln-10
  (math-read-number-simple "2.302585092994045684018")
  "An approximation for ln(10).")
     
(math-defcache math-ln-10 math-approx-ln-10
  (math-ln-raw-2 '(float 1 1)))

(defconst math-approx-ln-2
  (math-read-number-simple "0.693147180559945309417")
  "An approximation for ln(2).")

(math-defcache math-ln-2 math-approx-ln-2
  (math-ln-raw-3 (math-float '(frac 1 3))))



;;; Hyperbolic functions.

(defun calcFunc-sinh (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	(math-expand-formulas
	 (math-normalize
	  (list '/ (list '- (list 'calcFunc-exp x)
			 (list 'calcFunc-exp (list 'neg x))) 2)))
	((Math-numberp x)
	 (if calc-symbolic-mode (signal 'inexact-result nil))
	 (math-with-extra-prec 2
	   (let ((expx (math-exp-raw (math-float x))))
	     (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1)))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-sinh (nth 1 x))
			 (math-mul (nth 2 x) (calcFunc-cosh (nth 1 x)))))
	((eq (car x) 'intv)
	 (math-sort-intv (nth 1 x)
			 (calcFunc-sinh (nth 2 x))
			 (calcFunc-sinh (nth 3 x))))
	((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
	 (math-normalize
	  (list '/ (list '+ (list 'calcFunc-exp x)
			 (list 'calcFunc-exp (list 'neg x))) 2)))
	((Math-numberp x)
	 (if calc-symbolic-mode (signal 'inexact-result nil))
	 (math-with-extra-prec 2
	   (let ((expx (math-exp-raw (math-float x))))
	     (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1)))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-cosh (nth 1 x))
			 (math-mul (nth 2 x)
				   (calcFunc-sinh (nth 1 x)))))
	((and (eq (car x) 'intv) (math-intv-constp x))
	 (setq x (math-abs x))
	 (math-sort-intv (nth 1 x)
			 (calcFunc-cosh (nth 2 x))
			 (calcFunc-cosh (nth 3 x))))
	((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
	 (math-normalize
	  (let ((expx (list 'calcFunc-exp x))
		(expmx (list 'calcFunc-exp (list 'neg x))))
	    (math-normalize
	     (list '/ (list '- expx expmx) (list '+ expx expmx))))))
	((Math-numberp x)
	 (if calc-symbolic-mode (signal 'inexact-result nil))
	 (math-with-extra-prec 2
	   (let* ((expx (calcFunc-exp (math-float x)))
		  (expmx (math-div 1 expx)))
	     (math-div (math-sub expx expmx)
		       (math-add expx expmx)))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-tanh (nth 1 x))
			 (math-div (nth 2 x)
				   (math-sqr (calcFunc-cosh (nth 1 x))))))
	((eq (car x) 'intv)
	 (math-sort-intv (nth 1 x)
			 (calcFunc-tanh (nth 2 x))
			 (calcFunc-tanh (nth 3 x))))
	((equal x '(var inf var-inf))
	 1)
	((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-sech (x)   ; [N N] [Public]
  (cond ((eq x 0) 1)
	(math-expand-formulas
	 (math-normalize
	  (list '/ 2 (list '+ (list 'calcFunc-exp x)
                           (list 'calcFunc-exp (list 'neg x))))))
	((Math-numberp x)
	 (if calc-symbolic-mode (signal 'inexact-result nil))
	 (math-with-extra-prec 2
	   (let ((expx (math-exp-raw (math-float x))))
	     (math-div '(float 2 0) (math-add expx (math-div 1 expx))))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-sech (nth 1 x))
			 (math-mul (nth 2 x)
                                   (math-mul (calcFunc-sech (nth 1 x))
                                             (calcFunc-tanh (nth 1 x))))))
	((and (eq (car x) 'intv) (math-intv-constp x))
	 (setq x (math-abs x))
	 (math-sort-intv (nth 1 x)
			 (calcFunc-sech (nth 2 x))
			 (calcFunc-sech (nth 3 x))))
	((or (equal x '(var inf var-inf))
	     (equal x '(neg (var inf var-inf))))
         0)
        ((equal x '(var nan var-nan))
         x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-sech x))))
(put 'calcFunc-sech 'math-expandable t)

(defun calcFunc-csch (x)   ; [N N] [Public]
  (cond ((eq x 0) (math-div 1 0))
	(math-expand-formulas
	 (math-normalize
	  (list '/ 2 (list '- (list 'calcFunc-exp x)
                           (list 'calcFunc-exp (list 'neg x))))))
	((Math-numberp x)
	 (if calc-symbolic-mode (signal 'inexact-result nil))
	 (math-with-extra-prec 2
	   (let ((expx (math-exp-raw (math-float x))))
	     (math-div '(float 2 0) (math-add expx (math-div -1 expx))))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-csch (nth 1 x))
			 (math-mul (nth 2 x) 
                                   (math-mul (calcFunc-csch (nth 1 x))
                                             (calcFunc-coth (nth 1 x))))))
	((eq (car x) 'intv)
         (if (and (Math-negp (nth 2 x))
                  (Math-posp (nth 3 x)))
             '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
           (math-sort-intv (nth 1 x)
                           (calcFunc-csch (nth 2 x))
                           (calcFunc-csch (nth 3 x)))))
	((or (equal x '(var inf var-inf))
	     (equal x '(neg (var inf var-inf))))
         0)
        ((equal x '(var nan var-nan))
         x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-csch x))))
(put 'calcFunc-csch 'math-expandable t)

(defun calcFunc-coth (x)   ; [N N] [Public]
  (cond ((eq x 0) (math-div 1 0))
	(math-expand-formulas
	 (math-normalize
	  (let ((expx (list 'calcFunc-exp x))
		(expmx (list 'calcFunc-exp (list 'neg x))))
	    (math-normalize
	     (list '/ (list '+ expx expmx) (list '- expx expmx))))))
	((Math-numberp x)
	 (if calc-symbolic-mode (signal 'inexact-result nil))
	 (math-with-extra-prec 2
	   (let* ((expx (calcFunc-exp (math-float x)))
		  (expmx (math-div 1 expx)))
	     (math-div (math-add expx expmx)
		       (math-sub expx expmx)))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-coth (nth 1 x))
			 (math-div (nth 2 x)
				   (math-sqr (calcFunc-sinh (nth 1 x))))))
	((eq (car x) 'intv)
         (if (and (Math-negp (nth 2 x))
                  (Math-posp (nth 3 x)))
             '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
           (math-sort-intv (nth 1 x)
                           (calcFunc-coth (nth 2 x))
                           (calcFunc-coth (nth 3 x)))))
	((equal x '(var inf var-inf))
	 1)
	((equal x '(neg (var inf var-inf)))
	 -1)
	((equal x '(var nan var-nan))
	 x)
	(t (calc-record-why 'numberp x)
	   (list 'calcFunc-coth x))))
(put 'calcFunc-coth 'math-expandable t)

(defun calcFunc-arcsinh (x)   ; [N N] [Public]
  (cond ((eq x 0) 0)
	(math-expand-formulas
	 (math-normalize
	  (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
					      (list '+ (list '^ x 2) 1))))))
	((Math-numberp x)
	 (if calc-symbolic-mode (signal 'inexact-result nil))
	 (math-with-extra-prec 2
	   (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x)
							     '(float 1 0)))))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-arcsinh (nth 1 x))
			 (math-div (nth 2 x)
				   (math-sqrt
				    (math-add (math-sqr (nth 1 x)) 1)))))
	((eq (car x) 'intv)
	 (math-sort-intv (nth 1 x)
			 (calcFunc-arcsinh (nth 2 x))
			 (calcFunc-arcsinh (nth 3 x))))
	((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)
	 '(var pi var-pi))
	((and (eq x 0) calc-symbolic-mode)
	 (math-div (math-mul '(var pi var-pi) '(var i var-i)) 2))
	(math-expand-formulas
	 (math-normalize
	  (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
					      (list '- (list '^ x 2) 1))))))
	((Math-numberp x)
	 (if calc-symbolic-mode (signal 'inexact-result nil))
	 (if (Math-equal-int x -1)
	     (math-imaginary (math-pi))
	   (math-with-extra-prec 2
	     (if (or t    ; need to do this even in the real case!
		     (memq (car-safe x) '(cplx polar)))
		 (let ((xp1 (math-add 1 x)))  ; this gets the branch cuts right
		   (math-ln-raw
		    (math-add x (math-mul xp1
					  (math-sqrt-raw
					   (math-div (math-sub
						      x
						      '(float 1 0))
						     xp1))))))
	       (math-ln-raw
		(math-add x (math-sqrt-raw (math-add (math-sqr x)
						     '(float -1 0)))))))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-arccosh (nth 1 x))
			 (math-div (nth 2 x)
				   (math-sqrt
				    (math-add (math-sqr (nth 1 x)) -1)))))
	((eq (car x) 'intv)
	 (math-sort-intv (nth 1 x)
			 (calcFunc-arccosh (nth 2 x))
			 (calcFunc-arccosh (nth 3 x))))
	((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)
	 '(var inf var-inf))
	((and (Math-equal-int x -1) calc-infinite-mode)
	 '(neg (var inf var-inf)))
	(math-expand-formulas
	 (list '/ (list '-
			(list 'calcFunc-ln (list '+ 1 x))
			(list 'calcFunc-ln (list '- 1 x))) 2))
	((Math-numberp x)
	 (if calc-symbolic-mode (signal 'inexact-result nil))
	 (math-with-extra-prec 2
	   (if (or (memq (car-safe x) '(cplx polar))
		   (Math-lessp 1 x))
	       (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x))
				   (math-ln-raw (math-sub '(float 1 0) x)))
			 '(float 5 -1))
	     (if (and (math-equal-int x 1) calc-infinite-mode)
		 '(var inf var-inf)
	       (if (and (math-equal-int x -1) calc-infinite-mode)
		   '(neg (var inf var-inf))
		 (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x)
						  (math-sub 1 x)))
			   '(float 5 -1)))))))
	((eq (car-safe x) 'sdev)
	 (math-make-sdev (calcFunc-arctanh (nth 1 x))
			 (math-div (nth 2 x)
				   (math-sub 1 (math-sqr (nth 1 x))))))
	((eq (car x) 'intv)
	 (math-sort-intv (nth 1 x)
			 (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]
  (cond ((or (Math-numberp a)
	     (eq (car a) 'intv))
	 (math-with-extra-prec 2
	   (math-mul a (math-pi-over-180))))
	((eq (car a) 'hms)
	 (math-from-hms a 'rad))
	((eq (car a) 'sdev)
	 (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)
	     (eq (car a) 'intv))
	 (math-with-extra-prec 2
	   (math-div a (math-pi-over-180))))
	((eq (car a) 'hms)
	 (math-from-hms a 'deg))
	((eq (car a) 'sdev)
	 (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)

(provide 'calc-math)

;; arch-tag: c7367e8e-d0b8-4f70-8577-2fb3f31dbb4c
;;; calc-math.el ends here