lisp/calc/calc-math.el

   1 ;;; calc-math.el --- mathematical functions for Calc
   2
   3 ;; Copyright (C) 1990-1993, 2001-2011 Free Software Foundation, Inc.
   4
   5 ;; Author: David Gillespie <daveg@synaptics.com>
   6 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
   7
   8 ;; This file is part of GNU Emacs.
   9
  10 ;; GNU Emacs is free software: you can redistribute it and/or modify
  11 ;; it under the terms of the GNU General Public License as published by
  12 ;; the Free Software Foundation, either version 3 of the License, or
  13 ;; (at your option) any later version.
  14
  15 ;; GNU Emacs is distributed in the hope that it will be useful,
  16 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18 ;; GNU General Public License for more details.
  19
  20 ;; You should have received a copy of the GNU General Public License
  21 ;; along with GNU Emacs.  If not, see <http://www.gnu.org/licenses/>.
  22
  23 ;;; Commentary:
  24
  25 ;;; Code:
  26
  27 ;; This file is autoloaded from calc-ext.el.
  28
  29 (require 'calc-ext)
  30 (require 'calc-macs)
  31
  32
  33 ;;; Find out how many 9s in 9.9999... will give distinct Emacs floats,
  34 ;;; then back off by one.
  35
  36 (defvar math-emacs-precision
  37   (let* ((n 1)
  38          (x 9)
  39          (xx (+ x (* 9 (expt 10 (- n))))))
  40     (while (/= x xx)
  41       (progn
  42         (setq n (1+ n))
  43         (setq x xx)
  44         (setq xx (+ x (* 9 (expt 10 (- n)))))))
  45     (1- n))
  46   "The number of digits in an Emacs float.")
  47
  48 ;;; Find the largest power of 10 which is an Emacs float,
  49 ;;; then back off by one so that any float d.dddd...eN
  50 ;;; is an Emacs float, for acceptable d.dddd....
  51
  52 (defvar math-largest-emacs-expt
  53   (let ((x 1)
  54         (pow 1e2))
  55     ;; The following loop is for efficiency; it should stop when
  56     ;; 10^(2x) is too large.  This could be indicated by a range
  57     ;; error when computing 10^(2x) or an infinite value for 10^(2x).
  58     (while (and
  59             pow
  60             (< pow 1.0e+INF))
  61       (setq x (* 2 x))
  62       (setq pow (condition-case nil
  63                     (expt 10.0 (* 2 x))
  64                   (error nil))))
  65     ;; The following loop should stop when 10^(x+1) is too large.
  66     (setq pow (condition-case nil
  67                     (expt 10.0 (1+ x))
  68                   (error nil)))
  69     (while (and
  70             pow
  71             (< pow 1.0e+INF))
  72       (setq x (1+ x))
  73       (setq pow (condition-case nil
  74                     (expt 10.0 (1+ x))
  75                   (error nil))))
  76     (1- x))
  77   "The largest exponent which Calc will convert to an Emacs float.")
  78
  79 (defvar math-smallest-emacs-expt
  80   (let ((x -1))
  81     (while (condition-case nil
  82                (> (expt 10.0 x) 0.0)
  83              (error nil))
  84       (setq x (* 2 x)))
  85     (setq x (/ x 2))
  86     (while (condition-case nil
  87                (> (expt 10.0 x) 0.0)
  88              (error nil))
  89       (setq x (1- x)))
  90     (+ x 2))
  91     "The smallest exponent which Calc will convert to an Emacs float.")
  92
  93 (defun math-use-emacs-fn (fn x)
  94   "Use the native Emacs function FN to evaluate the Calc number X.
  95 If this can't be done, return NIL."
  96   (and
  97    (<= calc-internal-prec math-emacs-precision)
  98    (math-realp x)
  99    (let* ((fx (math-float x))
 100           (xpon (+ (nth 2 x) (1- (math-numdigs (nth 1 x))))))
 101      (and (<= math-smallest-emacs-expt xpon)
 102           (<= xpon math-largest-emacs-expt)
 103           (condition-case nil
 104               (math-read-number
 105                (number-to-string
 106                 (funcall fn
 107                          (string-to-number
 108                           (let
 109                               ((calc-number-radix 10)
 110                                (calc-twos-complement-mode nil)
 111                                (calc-float-format (list 'float calc-internal-prec))
 112                                (calc-group-digits nil)
 113                                (calc-point-char "."))
 114                             (math-format-number (math-float x)))))))
 115             (error nil))))))
 116
 117 (defun calc-sqrt (arg)
 118   (interactive "P")
 119   (calc-slow-wrapper
 120    (if (calc-is-inverse)
 121        (calc-unary-op "^2" 'calcFunc-sqr arg)
 122      (calc-unary-op "sqrt" 'calcFunc-sqrt arg))))
 123
 124 (defun calc-isqrt (arg)
 125   (interactive "P")
 126   (calc-slow-wrapper
 127    (if (calc-is-inverse)
 128        (calc-unary-op "^2" 'calcFunc-sqr arg)
 129      (calc-unary-op "isqt" 'calcFunc-isqrt arg))))
 130
 131
 132 (defun calc-hypot (arg)
 133   (interactive "P")
 134   (calc-slow-wrapper
 135    (calc-binary-op "hypt" 'calcFunc-hypot arg)))
 136
 137 (defun calc-ln (arg)
 138   (interactive "P")
 139   (calc-invert-func)
 140   (calc-exp arg))
 141
 142 (defun calc-log10 (arg)
 143   (interactive "P")
 144   (calc-hyperbolic-func)
 145   (calc-ln arg))
 146
 147 (defun calc-log (arg)
 148   (interactive "P")
 149   (calc-slow-wrapper
 150    (if (calc-is-inverse)
 151        (calc-binary-op "alog" 'calcFunc-alog arg)
 152      (calc-binary-op "log" 'calcFunc-log arg))))
 153
 154 (defun calc-ilog (arg)
 155   (interactive "P")
 156   (calc-slow-wrapper
 157    (if (calc-is-inverse)
 158        (calc-binary-op "alog" 'calcFunc-alog arg)
 159      (calc-binary-op "ilog" 'calcFunc-ilog arg))))
 160
 161 (defun calc-lnp1 (arg)
 162   (interactive "P")
 163   (calc-invert-func)
 164   (calc-expm1 arg))
 165
 166 (defun calc-exp (arg)
 167   (interactive "P")
 168   (calc-slow-wrapper
 169    (if (calc-is-hyperbolic)
 170        (if (calc-is-inverse)
 171            (calc-unary-op "lg10" 'calcFunc-log10 arg)
 172          (calc-unary-op "10^" 'calcFunc-exp10 arg))
 173      (if (calc-is-inverse)
 174          (calc-unary-op "ln" 'calcFunc-ln arg)
 175        (calc-unary-op "exp" 'calcFunc-exp arg)))))
 176
 177 (defun calc-expm1 (arg)
 178   (interactive "P")
 179   (calc-slow-wrapper
 180    (if (calc-is-inverse)
 181        (calc-unary-op "ln+1" 'calcFunc-lnp1 arg)
 182      (calc-unary-op "ex-1" 'calcFunc-expm1 arg))))
 183
 184 (defun calc-pi ()
 185   (interactive)
 186   (calc-slow-wrapper
 187    (if (calc-is-inverse)
 188        (if (calc-is-hyperbolic)
 189            (if calc-symbolic-mode
 190                (calc-pop-push-record 0 "phi" '(var phi var-phi))
 191              (calc-pop-push-record 0 "phi" (math-phi)))
 192          (if calc-symbolic-mode
 193              (calc-pop-push-record 0 "gmma" '(var gamma var-gamma))
 194            (calc-pop-push-record 0 "gmma" (math-gamma-const))))
 195      (if (calc-is-hyperbolic)
 196          (if calc-symbolic-mode
 197              (calc-pop-push-record 0 "e" '(var e var-e))
 198            (calc-pop-push-record 0 "e" (math-e)))
 199        (if calc-symbolic-mode
 200            (calc-pop-push-record 0 "pi" '(var pi var-pi))
 201          (calc-pop-push-record 0 "pi" (math-pi)))))))
 202
 203 (defun calc-sin (arg)
 204   (interactive "P")
 205   (calc-slow-wrapper
 206    (if (calc-is-hyperbolic)
 207        (if (calc-is-inverse)
 208            (calc-unary-op "asnh" 'calcFunc-arcsinh arg)
 209          (calc-unary-op "sinh" 'calcFunc-sinh arg))
 210      (if (calc-is-inverse)
 211          (calc-unary-op "asin" 'calcFunc-arcsin arg)
 212        (calc-unary-op "sin" 'calcFunc-sin arg)))))
 213
 214 (defun calc-arcsin (arg)
 215   (interactive "P")
 216   (calc-invert-func)
 217   (calc-sin arg))
 218
 219 (defun calc-sinh (arg)
 220   (interactive "P")
 221   (calc-hyperbolic-func)
 222   (calc-sin arg))
 223
 224 (defun calc-arcsinh (arg)
 225   (interactive "P")
 226   (calc-invert-func)
 227   (calc-hyperbolic-func)
 228   (calc-sin arg))
 229
 230 (defun calc-sec (arg)
 231   (interactive "P")
 232   (calc-slow-wrapper
 233    (if (calc-is-hyperbolic)
 234        (calc-unary-op "sech" 'calcFunc-sech arg)
 235      (calc-unary-op "sec" 'calcFunc-sec arg))))
 236
 237 (defun calc-sech (arg)
 238   (interactive "P")
 239   (calc-hyperbolic-func)
 240   (calc-sec arg))
 241
 242 (defun calc-cos (arg)
 243   (interactive "P")
 244   (calc-slow-wrapper
 245    (if (calc-is-hyperbolic)
 246        (if (calc-is-inverse)
 247            (calc-unary-op "acsh" 'calcFunc-arccosh arg)
 248          (calc-unary-op "cosh" 'calcFunc-cosh arg))
 249      (if (calc-is-inverse)
 250          (calc-unary-op "acos" 'calcFunc-arccos arg)
 251        (calc-unary-op "cos" 'calcFunc-cos arg)))))
 252
 253 (defun calc-arccos (arg)
 254   (interactive "P")
 255   (calc-invert-func)
 256   (calc-cos arg))
 257
 258 (defun calc-cosh (arg)
 259   (interactive "P")
 260   (calc-hyperbolic-func)
 261   (calc-cos arg))
 262
 263 (defun calc-arccosh (arg)
 264   (interactive "P")
 265   (calc-invert-func)
 266   (calc-hyperbolic-func)
 267   (calc-cos arg))
 268
 269 (defun calc-csc (arg)
 270   (interactive "P")
 271   (calc-slow-wrapper
 272    (if (calc-is-hyperbolic)
 273        (calc-unary-op "csch" 'calcFunc-csch arg)
 274      (calc-unary-op "csc" 'calcFunc-csc arg))))
 275
 276 (defun calc-csch (arg)
 277   (interactive "P")
 278   (calc-hyperbolic-func)
 279   (calc-csc arg))
 280
 281 (defun calc-sincos ()
 282   (interactive)
 283   (calc-slow-wrapper
 284    (if (calc-is-inverse)
 285        (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1)))
 286      (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1))))))
 287
 288 (defun calc-tan (arg)
 289   (interactive "P")
 290   (calc-slow-wrapper
 291    (if (calc-is-hyperbolic)
 292        (if (calc-is-inverse)
 293            (calc-unary-op "atnh" 'calcFunc-arctanh arg)
 294          (calc-unary-op "tanh" 'calcFunc-tanh arg))
 295      (if (calc-is-inverse)
 296          (calc-unary-op "atan" 'calcFunc-arctan arg)
 297        (calc-unary-op "tan" 'calcFunc-tan arg)))))
 298
 299 (defun calc-arctan (arg)
 300   (interactive "P")
 301   (calc-invert-func)
 302   (calc-tan arg))
 303
 304 (defun calc-tanh (arg)
 305   (interactive "P")
 306   (calc-hyperbolic-func)
 307   (calc-tan arg))
 308
 309 (defun calc-arctanh (arg)
 310   (interactive "P")
 311   (calc-invert-func)
 312   (calc-hyperbolic-func)
 313   (calc-tan arg))
 314
 315 (defun calc-cot (arg)
 316   (interactive "P")
 317   (calc-slow-wrapper
 318    (if (calc-is-hyperbolic)
 319        (calc-unary-op "coth" 'calcFunc-coth arg)
 320      (calc-unary-op "cot" 'calcFunc-cot arg))))
 321
 322 (defun calc-coth (arg)
 323   (interactive "P")
 324   (calc-hyperbolic-func)
 325   (calc-cot arg))
 326
 327 (defun calc-arctan2 ()
 328   (interactive)
 329   (calc-slow-wrapper
 330    (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2)))))
 331
 332 (defun calc-conj (arg)
 333   (interactive "P")
 334   (calc-wrapper
 335    (calc-unary-op "conj" 'calcFunc-conj arg)))
 336
 337 (defun calc-imaginary ()
 338   (interactive)
 339   (calc-slow-wrapper
 340    (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))))
 341
 342 (defun calc-to-degrees (arg)
 343   (interactive "P")
 344   (calc-wrapper
 345    (calc-unary-op ">deg" 'calcFunc-deg arg)))
 346
 347 (defun calc-to-radians (arg)
 348   (interactive "P")
 349   (calc-wrapper
 350    (calc-unary-op ">rad" 'calcFunc-rad arg)))
 351
 352
 353 (defun calc-degrees-mode (arg)
 354   (interactive "p")
 355   (cond ((= arg 1)
 356          (calc-wrapper
 357           (calc-change-mode 'calc-angle-mode 'deg)
 358           (message "Angles measured in degrees")))
 359         ((= arg 2) (calc-radians-mode))
 360         ((= arg 3) (calc-hms-mode))
 361         (t (error "Prefix argument out of range"))))
 362
 363 (defun calc-radians-mode ()
 364   (interactive)
 365   (calc-wrapper
 366    (calc-change-mode 'calc-angle-mode 'rad)
 367    (message "Angles measured in radians")))
 368
 369
 370 ;;; Compute the integer square-root floor(sqrt(A)).  A > 0.  [I I] [Public]
 371 ;;; This method takes advantage of the fact that Newton's method starting
 372 ;;; with an overestimate always works, even using truncating integer division!
 373 (defun math-isqrt (a)
 374   (cond ((Math-zerop a) a)
 375         ((not (math-natnump a))
 376          (math-reject-arg a 'natnump))
 377         ((integerp a)
 378          (math-isqrt-small a))
 379         (t
 380          (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a))))))))
 381
 382 (defun calcFunc-isqrt (a)
 383   (if (math-realp a)
 384       (math-isqrt (math-floor a))
 385     (math-floor (math-sqrt a))))
 386
 387
 388 ;;; This returns (flag . result) where the flag is t if A is a perfect square.
 389 (defun math-isqrt-bignum (a)   ; [P.l L]
 390   (let ((len (length a)))
 391     (if (= (% len 2) 0)
 392         (let* ((top (nthcdr (- len 2) a)))
 393           (math-isqrt-bignum-iter
 394            a
 395            (math-scale-bignum-digit-size
 396             (math-bignum-big
 397              (1+ (math-isqrt-small
 398                   (+ (* (nth 1 top) math-bignum-digit-size) (car top)))))
 399             (1- (/ len 2)))))
 400       (let* ((top (nth (1- len) a)))
 401         (math-isqrt-bignum-iter
 402          a
 403          (math-scale-bignum-digit-size
 404           (list (1+ (math-isqrt-small top)))
 405           (/ len 2)))))))
 406
 407 (defun math-isqrt-bignum-iter (a guess)   ; [l L l]
 408   (math-working "isqrt" (cons 'bigpos guess))
 409   (let* ((q (math-div-bignum a guess))
 410          (s (math-add-bignum (car q) guess))
 411          (g2 (math-div2-bignum s))
 412          (comp (math-compare-bignum g2 guess)))
 413     (if (< comp 0)
 414         (math-isqrt-bignum-iter a g2)
 415       (cons (and (= comp 0)
 416                  (math-zerop-bignum (cdr q))
 417                  (= (% (car s) 2) 0))
 418             guess))))
 419
 420 (defun math-zerop-bignum (a)
 421   (and (eq (car a) 0)
 422        (progn
 423          (while (eq (car (setq a (cdr a))) 0))
 424          (null a))))
 425
 426 (defun math-scale-bignum-digit-size (a n)   ; [L L S]
 427   (while (> n 0)
 428     (setq a (cons 0 a)
 429           n (1- n)))
 430   a)
 431
 432 (defun math-isqrt-small (a)   ; A > 0.  [S S]
 433   (let ((g (cond ((>= a 1000000) 10000)
 434                  ((>= a 10000) 1000)
 435                  ((>= a 100) 100)
 436                  (t 10)))
 437         g2)
 438     (while (< (setq g2 (/ (+ g (/ a g)) 2)) g)
 439       (setq g g2))
 440     g))
 441
 442
 443
 444
 445 ;;; Compute the square root of a number.
 446 ;;; [T N] if possible, else [F N] if possible, else [C N].  [Public]
 447 (defun math-sqrt (a)
 448   (or
 449    (and (Math-zerop a) a)
 450    (and (math-known-nonposp a)
 451         (math-imaginary (math-sqrt (math-neg a))))
 452    (and (integerp a)
 453         (let ((sqrt (math-isqrt-small a)))
 454           (if (= (* sqrt sqrt) a)
 455               sqrt
 456             (if calc-symbolic-mode
 457                 (list 'calcFunc-sqrt a)
 458               (math-sqrt-float (math-float a) (math-float sqrt))))))
 459    (and (eq (car-safe a) 'bigpos)
 460         (let* ((res (math-isqrt-bignum (cdr a)))
 461                (sqrt (math-normalize (cons 'bigpos (cdr res)))))
 462           (if (car res)
 463               sqrt
 464             (if calc-symbolic-mode
 465                 (list 'calcFunc-sqrt a)
 466               (math-sqrt-float (math-float a) (math-float sqrt))))))
 467    (and (eq (car-safe a) 'frac)
 468         (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a)))))
 469                (num-sqrt (math-normalize (cons 'bigpos (cdr num-res))))
 470                (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a)))))
 471                (den-sqrt (math-normalize (cons 'bigpos (cdr den-res)))))
 472           (if (and (car num-res) (car den-res))
 473               (list 'frac num-sqrt den-sqrt)
 474             (if calc-symbolic-mode
 475                 (if (or (car num-res) (car den-res))
 476                     (math-div (if (car num-res)
 477                                   num-sqrt (list 'calcFunc-sqrt (nth 1 a)))
 478                               (if (car den-res)
 479                                   den-sqrt (list 'calcFunc-sqrt (nth 2 a))))
 480                   (list 'calcFunc-sqrt a))
 481               (math-sqrt-float (math-float a)
 482                                (math-div (math-float num-sqrt) den-sqrt))))))
 483    (and (eq (car-safe a) 'float)
 484         (if calc-symbolic-mode
 485             (if (= (% (nth 2 a) 2) 0)
 486                 (let ((res (math-isqrt-bignum
 487                             (cdr (Math-bignum-test (nth 1 a))))))
 488                   (if (car res)
 489                       (math-make-float (math-normalize
 490                                         (cons 'bigpos (cdr res)))
 491                                        (/ (nth 2 a) 2))
 492                     (signal 'inexact-result nil)))
 493               (signal 'inexact-result nil))
 494           (math-sqrt-float a)))
 495    (and (eq (car-safe a) 'cplx)
 496         (math-with-extra-prec 2
 497           (let* ((d (math-abs a))
 498                  (imag (math-sqrt (math-mul (math-sub d (nth 1 a))
 499                                             '(float 5 -1)))))
 500             (list 'cplx
 501                   (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1)))
 502                   (if (math-negp (nth 2 a)) (math-neg imag) imag)))))
 503    (and (eq (car-safe a) 'polar)
 504         (list 'polar
 505               (math-sqrt (nth 1 a))
 506               (math-mul (nth 2 a) '(float 5 -1))))
 507    (and (eq (car-safe a) 'sdev)
 508         (let ((sqrt (math-sqrt (nth 1 a))))
 509           (math-make-sdev sqrt
 510                           (math-div (nth 2 a) (math-mul sqrt 2)))))
 511    (and (eq (car-safe a) 'intv)
 512         (not (math-negp (nth 2 a)))
 513         (math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a))))
 514    (and (eq (car-safe a) '*)
 515         (or (math-known-nonnegp (nth 1 a))
 516             (math-known-nonnegp (nth 2 a)))
 517         (math-mul (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
 518    (and (eq (car-safe a) '/)
 519         (or (and (math-known-nonnegp (nth 2 a))
 520                  (math-div (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
 521             (and (math-known-nonnegp (nth 1 a))
 522                  (not (math-equal-int (nth 1 a) 1))
 523                  (math-mul (math-sqrt (nth 1 a))
 524                            (math-sqrt (math-div 1 (nth 2 a)))))))
 525    (and (eq (car-safe a) '^)
 526         (math-known-evenp (nth 2 a))
 527         (math-known-realp (nth 1 a))
 528         (math-abs (math-pow (nth 1 a) (math-div (nth 2 a) 2))))
 529    (let ((inf (math-infinitep a)))
 530      (and inf
 531           (math-mul (math-sqrt (math-infinite-dir a inf)) inf)))
 532    (progn
 533      (calc-record-why 'numberp a)
 534      (list 'calcFunc-sqrt a))))
 535 (defalias 'calcFunc-sqrt 'math-sqrt)
 536
 537 (defun math-infinite-dir (a &optional inf)
 538   (or inf (setq inf (math-infinitep a)))
 539   (math-normalize (math-expr-subst a inf 1)))
 540
 541 (defun math-sqrt-float (a &optional guess)   ; [F F F]
 542   (if calc-symbolic-mode
 543       (signal 'inexact-result nil)
 544     (math-with-extra-prec 1 (math-sqrt-raw a guess))))
 545
 546 (defun math-sqrt-raw (a &optional guess)   ; [F F F]
 547   (if (not (Math-posp a))
 548       (math-sqrt a)
 549     (cond
 550      ((math-use-emacs-fn 'sqrt a))
 551      (t
 552       (if (null guess)
 553           (let ((ldiff (- (math-numdigs (nth 1 a)) 6)))
 554             (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff)))
 555             (setq guess (math-make-float (math-isqrt-small
 556                                           (math-scale-int (nth 1 a) (- ldiff)))
 557                                          (/ (+ (nth 2 a) ldiff) 2)))))
 558       (math-sqrt-float-iter a guess)))))
 559
 560 (defun math-sqrt-float-iter (a guess)   ; [F F F]
 561   (math-working "sqrt" guess)
 562   (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess))
 563                             '(float 5 -1))))
 564      (if (math-nearly-equal-float g2 guess)
 565          g2
 566        (math-sqrt-float-iter a g2))))
 567
 568 ;;; True if A and B differ only in the last digit of precision.  [P F F]
 569 (defun math-nearly-equal-float (a b)
 570   (let ((ediff (- (nth 2 a) (nth 2 b))))
 571     (cond ((= ediff 0)   ;; Expanded out for speed
 572            (setq ediff (math-add (Math-integer-neg (nth 1 a)) (nth 1 b)))
 573            (or (eq ediff 0)
 574                (and (not (consp ediff))
 575                     (< ediff 10)
 576                     (> ediff -10)
 577                     (= (math-numdigs (nth 1 a)) calc-internal-prec))))
 578           ((= ediff 1)
 579            (setq ediff (math-add (Math-integer-neg (nth 1 b))
 580                                  (math-scale-int (nth 1 a) 1)))
 581            (and (not (consp ediff))
 582                 (< ediff 10)
 583                 (> ediff -10)
 584                 (= (math-numdigs (nth 1 b)) calc-internal-prec)))
 585           ((= ediff -1)
 586            (setq ediff (math-add (Math-integer-neg (nth 1 a))
 587                                  (math-scale-int (nth 1 b) 1)))
 588            (and (not (consp ediff))
 589                 (< ediff 10)
 590                 (> ediff -10)
 591                 (= (math-numdigs (nth 1 a)) calc-internal-prec))))))
 592
 593 (defun math-nearly-equal (a b)   ;  [P N N] [Public]
 594   (setq a (math-float a))
 595   (setq b (math-float b))
 596   (if (eq (car a) 'polar) (setq a (math-complex a)))
 597   (if (eq (car b) 'polar) (setq b (math-complex b)))
 598   (if (eq (car a) 'cplx)
 599       (if (eq (car b) 'cplx)
 600           (and (or (math-nearly-equal-float (nth 1 a) (nth 1 b))
 601                    (and (math-nearly-zerop-float (nth 1 a) (nth 2 a))
 602                         (math-nearly-zerop-float (nth 1 b) (nth 2 b))))
 603                (or (math-nearly-equal-float (nth 2 a) (nth 2 b))
 604                    (and (math-nearly-zerop-float (nth 2 a) (nth 1 a))
 605                         (math-nearly-zerop-float (nth 2 b) (nth 1 b)))))
 606         (and (math-nearly-equal-float (nth 1 a) b)
 607              (math-nearly-zerop-float (nth 2 a) b)))
 608       (if (eq (car b) 'cplx)
 609           (and (math-nearly-equal-float a (nth 1 b))
 610                (math-nearly-zerop-float a (nth 2 b)))
 611         (math-nearly-equal-float a b))))
 612
 613 ;;; True if A is nearly zero compared to B.  [P F F]
 614 (defun math-nearly-zerop-float (a b)
 615   (or (eq (nth 1 a) 0)
 616       (<= (+ (math-numdigs (nth 1 a)) (nth 2 a))
 617           (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec)))))
 618
 619 (defun math-nearly-zerop (a b)   ; [P N R] [Public]
 620   (setq a (math-float a))
 621   (setq b (math-float b))
 622   (if (eq (car a) 'cplx)
 623       (and (math-nearly-zerop-float (nth 1 a) b)
 624            (math-nearly-zerop-float (nth 2 a) b))
 625     (if (eq (car a) 'polar)
 626         (math-nearly-zerop-float (nth 1 a) b)
 627       (math-nearly-zerop-float a b))))
 628
 629 ;;; This implementation could be improved, accuracy-wise.
 630 (defun math-hypot (a b)
 631   (cond ((Math-zerop a) (math-abs b))
 632         ((Math-zerop b) (math-abs a))
 633         ((not (Math-scalarp a))
 634          (if (math-infinitep a)
 635              (if (math-infinitep b)
 636                  (if (equal a b)
 637                      a
 638                    '(var nan var-nan))
 639                a)
 640            (calc-record-why 'scalarp a)
 641            (list 'calcFunc-hypot a b)))
 642         ((not (Math-scalarp b))
 643          (if (math-infinitep b)
 644              b
 645            (calc-record-why 'scalarp b)
 646            (list 'calcFunc-hypot a b)))
 647         ((and (Math-numberp a) (Math-numberp b))
 648          (math-with-extra-prec 1
 649            (math-sqrt (math-add (calcFunc-abssqr a) (calcFunc-abssqr b)))))
 650         ((eq (car-safe a) 'hms)
 651          (if (eq (car-safe b) 'hms)   ; this helps sdev's of hms forms
 652              (math-to-hms (math-hypot (math-from-hms a 'deg)
 653                                       (math-from-hms b 'deg)))
 654            (math-to-hms (math-hypot (math-from-hms a 'deg) b))))
 655         ((eq (car-safe b) 'hms)
 656          (math-to-hms (math-hypot a (math-from-hms b 'deg))))
 657         (t nil)))
 658 (defalias 'calcFunc-hypot 'math-hypot)
 659
 660 (defun calcFunc-sqr (x)
 661   (math-pow x 2))
 662
 663
 664
 665 (defun math-nth-root (a n)
 666   (cond ((= n 2) (math-sqrt a))
 667         ((Math-zerop a) a)
 668         ((Math-negp a) nil)
 669         ((Math-integerp a)
 670          (let ((root (math-nth-root-integer a n)))
 671            (if (car root)
 672                (cdr root)
 673              (and (not calc-symbolic-mode)
 674                   (math-nth-root-float (math-float a) n
 675                                        (math-float (cdr root)))))))
 676         ((eq (car-safe a) 'frac)
 677          (let* ((num-root (math-nth-root-integer (nth 1 a) n))
 678                 (den-root (math-nth-root-integer (nth 2 a) n)))
 679            (if (and (car num-root) (car den-root))
 680                (list 'frac (cdr num-root) (cdr den-root))
 681              (and (not calc-symbolic-mode)
 682                   (math-nth-root-float
 683                    (math-float a) n
 684                    (math-div-float (math-float (cdr num-root))
 685                                    (math-float (cdr den-root))))))))
 686         ((eq (car-safe a) 'float)
 687          (and (not calc-symbolic-mode)
 688               (math-nth-root-float a n)))
 689         ((eq (car-safe a) 'polar)
 690          (let ((root (math-nth-root (nth 1 a) n)))
 691            (and root (list 'polar root (math-div (nth 2 a) n)))))
 692         (t nil)))
 693
 694 ;; The variables math-nrf-n, math-nrf-nf and math-nrf-nfm1 are local
 695 ;; to math-nth-root-float, but are used by math-nth-root-float-iter,
 696 ;; which is called by math-nth-root-float.
 697 (defvar math-nrf-n)
 698 (defvar math-nrf-nf)
 699 (defvar math-nrf-nfm1)
 700
 701 (defun math-nth-root-float (a math-nrf-n &optional guess)
 702   (math-inexact-result)
 703   (math-with-extra-prec 1
 704     (let ((math-nrf-nf (math-float math-nrf-n))
 705           (math-nrf-nfm1 (math-float (1- math-nrf-n))))
 706       (math-nth-root-float-iter a (or guess
 707                                       (math-make-float
 708                                        1 (/ (+ (math-numdigs (nth 1 a))
 709                                                (nth 2 a)
 710                                                (/ math-nrf-n 2))
 711                                             math-nrf-n)))))))
 712
 713 (defun math-nth-root-float-iter (a guess)
 714   (math-working "root" guess)
 715   (let ((g2 (math-div-float (math-add-float (math-mul math-nrf-nfm1 guess)
 716                                             (math-div-float
 717                                              a (math-ipow guess (1- math-nrf-n))))
 718                             math-nrf-nf)))
 719     (if (math-nearly-equal-float g2 guess)
 720         g2
 721       (math-nth-root-float-iter a g2))))
 722
 723 ;; The variable math-nri-n is local to math-nth-root-integer, but
 724 ;; is used by math-nth-root-int-iter, which is called by
 725 ;; math-nth-root-int.
 726 (defvar math-nri-n)
 727
 728 (defun math-nth-root-integer (a math-nri-n &optional guess)   ; [I I S]
 729   (math-nth-root-int-iter a (or guess
 730                                 (math-scale-int 1 (/ (+ (math-numdigs a)
 731                                                         (1- math-nri-n))
 732                                                      math-nri-n)))))
 733
 734 (defun math-nth-root-int-iter (a guess)
 735   (math-working "root" guess)
 736   (let* ((q (math-idivmod a (math-ipow guess (1- math-nri-n))))
 737          (s (math-add (car q) (math-mul (1- math-nri-n) guess)))
 738          (g2 (math-idivmod s math-nri-n)))
 739     (if (Math-natnum-lessp (car g2) guess)
 740         (math-nth-root-int-iter a (car g2))
 741       (cons (and (equal (car g2) guess)
 742                  (eq (cdr q) 0)
 743                  (eq (cdr g2) 0))
 744             guess))))
 745
 746 (defun calcFunc-nroot (x n)
 747   (calcFunc-pow x (if (integerp n)
 748                       (math-make-frac 1 n)
 749                     (math-div 1 n))))
 750
 751
 752
 753
 754 ;;;; Transcendental functions.
 755
 756 ;;; All of these functions are defined on the complex plane.
 757 ;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
 758
 759 ;;; Most functions increase calc-internal-prec by 2 digits, then round
 760 ;;; down afterward.  "-raw" functions use the current precision, require
 761 ;;; their arguments to be in float (or complex float) format, and always
 762 ;;; work in radians (where applicable).
 763
 764 (defun math-to-radians (a)   ; [N N]
 765   (cond ((eq (car-safe a) 'hms)
 766          (math-from-hms a 'rad))
 767         ((memq calc-angle-mode '(deg hms))
 768          (math-mul a (math-pi-over-180)))
 769         (t a)))
 770
 771 (defun math-from-radians (a)   ; [N N]
 772   (cond ((eq calc-angle-mode 'deg)
 773          (if (math-constp a)
 774              (math-div a (math-pi-over-180))
 775            (list 'calcFunc-deg a)))
 776         ((eq calc-angle-mode 'hms)
 777          (math-to-hms a 'rad))
 778         (t a)))
 779
 780 (defun math-to-radians-2 (a)   ; [N N]
 781   (cond ((eq (car-safe a) 'hms)
 782          (math-from-hms a 'rad))
 783         ((memq calc-angle-mode '(deg hms))
 784          (if calc-symbolic-mode
 785              (math-div (math-mul a '(var pi var-pi)) 180)
 786            (math-mul a (math-pi-over-180))))
 787         (t a)))
 788
 789 (defun math-from-radians-2 (a)   ; [N N]
 790   (cond ((memq calc-angle-mode '(deg hms))
 791          (if calc-symbolic-mode
 792              (math-div (math-mul 180 a) '(var pi var-pi))
 793            (math-div a (math-pi-over-180))))
 794         (t a)))
 795
 796
 797
 798 ;;; Sine, cosine, and tangent.
 799
 800 (defun calcFunc-sin (x)   ; [N N] [Public]
 801   (cond ((and (integerp x)
 802               (if (eq calc-angle-mode 'deg)
 803                   (= (% x 90) 0)
 804                 (= x 0)))
 805          (aref [0 1 0 -1] (math-mod (/ x 90) 4)))
 806         ((Math-scalarp x)
 807          (math-with-extra-prec 2
 808            (math-sin-raw (math-to-radians (math-float x)))))
 809         ((eq (car x) 'sdev)
 810          (if (math-constp x)
 811              (math-with-extra-prec 2
 812                (let* ((xx (math-to-radians (math-float (nth 1 x))))
 813                       (xs (math-to-radians (math-float (nth 2 x))))
 814                       (sc (math-sin-cos-raw xx)))
 815                  (math-make-sdev (car sc) (math-mul xs (cdr sc)))))
 816            (math-make-sdev (calcFunc-sin (nth 1 x))
 817                            (math-mul (nth 2 x) (calcFunc-cos (nth 1 x))))))
 818         ((and (eq (car x) 'intv) (math-intv-constp x))
 819          (calcFunc-cos (math-sub x (math-quarter-circle nil))))
 820         ((equal x '(var nan var-nan))
 821          x)
 822         (t (calc-record-why 'scalarp x)
 823            (list 'calcFunc-sin x))))
 824
 825 (defun calcFunc-cos (x)   ; [N N] [Public]
 826   (cond ((and (integerp x)
 827               (if (eq calc-angle-mode 'deg)
 828                   (= (% x 90) 0)
 829                 (= x 0)))
 830          (aref [1 0 -1 0] (math-mod (/ x 90) 4)))
 831         ((Math-scalarp x)
 832          (math-with-extra-prec 2
 833            (math-cos-raw (math-to-radians (math-float x)))))
 834         ((eq (car x) 'sdev)
 835          (if (math-constp x)
 836              (math-with-extra-prec 2
 837                (let* ((xx (math-to-radians (math-float (nth 1 x))))
 838                       (xs (math-to-radians (math-float (nth 2 x))))
 839                       (sc (math-sin-cos-raw xx)))
 840                  (math-make-sdev (cdr sc) (math-mul xs (car sc)))))
 841            (math-make-sdev (calcFunc-cos (nth 1 x))
 842                            (math-mul (nth 2 x) (calcFunc-sin (nth 1 x))))))
 843         ((and (eq (car x) 'intv) (math-intv-constp x))
 844          (math-with-extra-prec 2
 845            (let* ((xx (math-to-radians (math-float x)))
 846                   (na (math-floor (math-div (nth 2 xx) (math-pi))))
 847                   (nb (math-floor (math-div (nth 3 xx) (math-pi))))
 848                   (span (math-sub nb na)))
 849              (if (memq span '(0 1))
 850                  (let ((int (math-sort-intv (nth 1 x)
 851                                             (math-cos-raw (nth 2 xx))
 852                                             (math-cos-raw (nth 3 xx)))))
 853                    (if (eq span 1)
 854                        (if (math-evenp na)
 855                            (math-make-intv (logior (nth 1 x) 2)
 856                                            -1
 857                                            (nth 3 int))
 858                          (math-make-intv (logior (nth 1 x) 1)
 859                                          (nth 2 int)
 860                                          1))
 861                      int))
 862                (list 'intv 3 -1 1)))))
 863         ((equal x '(var nan var-nan))
 864          x)
 865         (t (calc-record-why 'scalarp x)
 866            (list 'calcFunc-cos x))))
 867
 868 (defun calcFunc-sincos (x)   ; [V N] [Public]
 869   (if (Math-scalarp x)
 870       (math-with-extra-prec 2
 871         (let ((sc (math-sin-cos-raw (math-to-radians (math-float x)))))
 872           (list 'vec (cdr sc) (car sc))))    ; the vector [cos, sin]
 873     (list 'vec (calcFunc-sin x) (calcFunc-cos x))))
 874
 875 (defun calcFunc-tan (x)   ; [N N] [Public]
 876   (cond ((and (integerp x)
 877               (if (eq calc-angle-mode 'deg)
 878                   (= (% x 180) 0)
 879                 (= x 0)))
 880          0)
 881         ((Math-scalarp x)
 882          (math-with-extra-prec 2
 883            (math-tan-raw (math-to-radians (math-float x)))))
 884         ((eq (car x) 'sdev)
 885          (if (math-constp x)
 886              (math-with-extra-prec 2
 887                (let* ((xx (math-to-radians (math-float (nth 1 x))))
 888                       (xs (math-to-radians (math-float (nth 2 x))))
 889                       (sc (math-sin-cos-raw xx)))
 890                  (if (and (math-zerop (cdr sc)) (not calc-infinite-mode))
 891                      (progn
 892                        (calc-record-why "*Division by zero")
 893                        (list 'calcFunc-tan x))
 894                    (math-make-sdev (math-div-float (car sc) (cdr sc))
 895                                    (math-div-float xs (math-sqr (cdr sc)))))))
 896            (math-make-sdev (calcFunc-tan (nth 1 x))
 897                            (math-div (nth 2 x)
 898                                      (math-sqr (calcFunc-cos (nth 1 x)))))))
 899         ((and (eq (car x) 'intv) (math-intv-constp x))
 900          (or (math-with-extra-prec 2
 901                (let* ((xx (math-to-radians (math-float x)))
 902                       (na (math-floor (math-div (math-sub (nth 2 xx)
 903                                                           (math-pi-over-2))
 904                                                 (math-pi))))
 905                       (nb (math-floor (math-div (math-sub (nth 3 xx)
 906                                                           (math-pi-over-2))
 907                                                 (math-pi)))))
 908                  (and (equal na nb)
 909                       (math-sort-intv (nth 1 x)
 910                                       (math-tan-raw (nth 2 xx))
 911                                       (math-tan-raw (nth 3 xx))))))
 912              '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
 913         ((equal x '(var nan var-nan))
 914          x)
 915         (t (calc-record-why 'scalarp x)
 916            (list 'calcFunc-tan x))))
 917
 918 (defun calcFunc-sec (x)
 919   (cond ((and (integerp x)
 920               (eq calc-angle-mode 'deg)
 921               (= (% x 180) 0))
 922          (if (= (% x 360) 0)
 923              1
 924            -1))
 925         ((and (integerp x)
 926               (eq calc-angle-mode 'rad)
 927               (= x 0))
 928          1)
 929         ((Math-scalarp x)
 930          (math-with-extra-prec 2
 931            (math-sec-raw (math-to-radians (math-float x)))))
 932         ((eq (car x) 'sdev)
 933          (if (math-constp x)
 934              (math-with-extra-prec 2
 935                (let* ((xx (math-to-radians (math-float (nth 1 x))))
 936                       (xs (math-to-radians (math-float (nth 2 x))))
 937                       (sc (math-sin-cos-raw xx)))
 938                  (if (and (math-zerop (cdr sc))
 939                           (not calc-infinite-mode))
 940                      (progn
 941                        (calc-record-why "*Division by zero")
 942                        (list 'calcFunc-sec x))
 943                    (math-make-sdev (math-div-float '(float 1 0) (cdr sc))
 944                                    (math-div-float
 945                                     (math-mul xs (car sc))
 946                                     (math-sqr (cdr sc)))))))
 947            (math-make-sdev (calcFunc-sec (nth 1 x))
 948                            (math-div
 949                             (math-mul (nth 2 x)
 950                                       (calcFunc-sin (nth 1 x)))
 951                             (math-sqr (calcFunc-cos (nth 1 x)))))))
 952         ((and (eq (car x) 'intv)
 953               (math-intv-constp x))
 954          (math-with-extra-prec 2
 955            (let* ((xx (math-to-radians (math-float x)))
 956                   (na (math-floor (math-div (math-sub (nth 2 xx)
 957                                                       (math-pi-over-2))
 958                                             (math-pi))))
 959                   (nb (math-floor (math-div (math-sub (nth 3 xx)
 960                                                       (math-pi-over-2))
 961                                             (math-pi))))
 962                   (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2))))
 963                   (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2))))
 964                   (span (math-sub nbb naa)))
 965              (if (not (equal na nb))
 966                  '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
 967                (let ((int (math-sort-intv (nth 1 x)
 968                                           (math-sec-raw (nth 2 xx))
 969                                           (math-sec-raw (nth 3 xx)))))
 970                  (if (eq span 1)
 971                      (if (math-evenp (math-div (math-add naa 1) 2))
 972                          (math-make-intv (logior (nth 1 int) 2)
 973                                          1
 974                                          (nth 3 int))
 975                        (math-make-intv (logior (nth 1 int) 1)
 976                                        (nth 2 int)
 977                                        -1))
 978                    int))))))
 979         ((equal x '(var nan var-nan))
 980          x)
 981         (t (calc-record-why 'scalarp x)
 982            (list 'calcFunc-sec x))))
 983
 984 (defun calcFunc-csc (x)
 985   (cond ((and (integerp x)
 986               (eq calc-angle-mode 'deg)
 987               (= (% (- x 90) 180) 0))
 988          (if (= (% (- x 90) 360) 0)
 989              1
 990            -1))
 991         ((Math-scalarp x)
 992          (math-with-extra-prec 2
 993            (math-csc-raw (math-to-radians (math-float x)))))
 994         ((eq (car x) 'sdev)
 995          (if (math-constp x)
 996              (math-with-extra-prec 2
 997                (let* ((xx (math-to-radians (math-float (nth 1 x))))
 998                       (xs (math-to-radians (math-float (nth 2 x))))
 999                       (sc (math-sin-cos-raw xx)))
1000                  (if (and (math-zerop (car sc))
1001                           (not calc-infinite-mode))
1002                      (progn
1003                        (calc-record-why "*Division by zero")
1004                        (list 'calcFunc-csc x))
1005                    (math-make-sdev (math-div-float '(float 1 0) (car sc))
1006                                    (math-div-float
1007                                     (math-mul xs (cdr sc))
1008                                     (math-sqr (car sc)))))))
1009            (math-make-sdev (calcFunc-csc (nth 1 x))
1010                            (math-div
1011                             (math-mul (nth 2 x)
1012                                       (calcFunc-cos (nth 1 x)))
1013                             (math-sqr (calcFunc-sin (nth 1 x)))))))
1014         ((and (eq (car x) 'intv)
1015               (math-intv-constp x))
1016          (math-with-extra-prec 2
1017            (let* ((xx (math-to-radians (math-float x)))
1018                   (na (math-floor (math-div (nth 2 xx) (math-pi))))
1019                   (nb (math-floor (math-div (nth 3 xx) (math-pi))))
1020                   (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2))))
1021                   (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2))))
1022                   (span (math-sub nbb naa)))
1023              (if (not (equal na nb))
1024                  '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1025                (let ((int (math-sort-intv (nth 1 x)
1026                                           (math-csc-raw (nth 2 xx))
1027                                           (math-csc-raw (nth 3 xx)))))
1028                  (if (eq span 1)
1029                      (if (math-evenp (math-div naa 2))
1030                          (math-make-intv (logior (nth 1 int) 2)
1031                                          1
1032                                          (nth 3 int))
1033                        (math-make-intv (logior (nth 1 int) 1)
1034                                        (nth 2 int)
1035                                        -1))
1036                    int))))))
1037         ((equal x '(var nan var-nan))
1038          x)
1039         (t (calc-record-why 'scalarp x)
1040            (list 'calcFunc-csc x))))
1041
1042 (defun calcFunc-cot (x)   ; [N N] [Public]
1043   (cond ((and (integerp x)
1044               (if (eq calc-angle-mode 'deg)
1045                   (= (% (- x 90) 180) 0)
1046                 (= x 0)))
1047          0)
1048         ((Math-scalarp x)
1049          (math-with-extra-prec 2
1050            (math-cot-raw (math-to-radians (math-float x)))))
1051         ((eq (car x) 'sdev)
1052          (if (math-constp x)
1053              (math-with-extra-prec 2
1054                (let* ((xx (math-to-radians (math-float (nth 1 x))))
1055                       (xs (math-to-radians (math-float (nth 2 x))))
1056                       (sc (math-sin-cos-raw xx)))
1057                  (if (and (math-zerop (car sc)) (not calc-infinite-mode))
1058                      (progn
1059                        (calc-record-why "*Division by zero")
1060                        (list 'calcFunc-cot x))
1061                    (math-make-sdev (math-div-float (cdr sc) (car sc))
1062                                    (math-div-float xs (math-sqr (car sc)))))))
1063            (math-make-sdev (calcFunc-cot (nth 1 x))
1064                            (math-div (nth 2 x)
1065                                      (math-sqr (calcFunc-sin (nth 1 x)))))))
1066         ((and (eq (car x) 'intv) (math-intv-constp x))
1067          (or (math-with-extra-prec 2
1068                (let* ((xx (math-to-radians (math-float x)))
1069                       (na (math-floor (math-div (nth 2 xx) (math-pi))))
1070                       (nb (math-floor (math-div (nth 3 xx) (math-pi)))))
1071                  (and (equal na nb)
1072                       (math-sort-intv (nth 1 x)
1073                                       (math-cot-raw (nth 2 xx))
1074                                       (math-cot-raw (nth 3 xx))))))
1075              '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
1076         ((equal x '(var nan var-nan))
1077          x)
1078         (t (calc-record-why 'scalarp x)
1079            (list 'calcFunc-cot x))))
1080
1081 (defun math-sin-raw (x &optional orgx)   ; [N N]
1082   (cond ((eq (car x) 'cplx)
1083          (let* ((expx (math-exp-raw (nth 2 x)))
1084                 (expmx (math-div-float '(float 1 0) expx))
1085                 (sc (math-sin-cos-raw (nth 1 x))))
1086            (list 'cplx
1087                  (math-mul-float (car sc)
1088                                  (math-mul-float (math-add-float expx expmx)
1089                                                  '(float 5 -1)))
1090                  (math-mul-float (cdr sc)
1091                                  (math-mul-float (math-sub-float expx expmx)
1092                                                  '(float 5 -1))))))
1093         ((eq (car x) 'polar)
1094          (math-polar (math-sin-raw (math-complex x))))
1095         ((Math-integer-negp (nth 1 x))
1096          (math-neg-float (math-sin-raw (math-neg-float x) (if orgx orgx x))))
1097         ((math-lessp-float '(float 7 0) x)  ; avoid inf loops due to roundoff
1098          (math-sin-raw (math-mod x (math-two-pi)) (if orgx orgx x)))
1099         (t (math-sin-raw-2 x (if orgx orgx x)))))
1100
1101 (defun math-cos-raw (x)   ; [N N]
1102   (if (eq (car-safe x) 'polar)
1103       (math-polar (math-cos-raw (math-complex x)))
1104     (math-sin-raw (math-sub (math-pi-over-2) x) x)))
1105
1106 (defun math-sec-raw (x)   ; [N N]
1107   (cond ((eq (car x) 'cplx)
1108          (let* ((x (math-mul x '(float 1 0)))
1109                 (expx (math-exp-raw (nth 2 x)))
1110                 (expmx (math-div-float '(float 1 0) expx))
1111                 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
1112                 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
1113                 (sc (math-sin-cos-raw (nth 1 x)))
1114                 (d (math-add-float
1115                     (math-mul-float (math-sqr (car sc))
1116                                     (math-sqr sh))
1117                     (math-mul-float (math-sqr (cdr sc))
1118                                     (math-sqr ch)))))
1119            (and (not (eq (nth 1 d) 0))
1120                 (list 'cplx
1121                       (math-div-float (math-mul-float (cdr sc) ch) d)
1122                       (math-div-float (math-mul-float (car sc) sh) d)))))
1123         ((eq (car x) 'polar)
1124          (math-polar (math-sec-raw (math-complex x))))
1125         (t
1126          (let ((cs (math-cos-raw x)))
1127            (if (eq cs 0)
1128                (math-div 1 0)
1129              (math-div-float '(float 1 0) cs))))))
1130
1131 (defun math-csc-raw (x)   ; [N N]
1132   (cond ((eq (car x) 'cplx)
1133          (let* ((x (math-mul x '(float 1 0)))
1134                 (expx (math-exp-raw (nth 2 x)))
1135                 (expmx (math-div-float '(float 1 0) expx))
1136                 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
1137                 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
1138                 (sc (math-sin-cos-raw (nth 1 x)))
1139                 (d (math-add-float
1140                     (math-mul-float (math-sqr (car sc))
1141                                     (math-sqr ch))
1142                     (math-mul-float (math-sqr (cdr sc))
1143                                     (math-sqr sh)))))
1144            (and (not (eq (nth 1 d) 0))
1145                 (list 'cplx
1146                       (math-div-float (math-mul-float (car sc) ch) d)
1147                       (math-div-float (math-mul-float (cdr sc) sh) d)))))
1148         ((eq (car x) 'polar)
1149          (math-polar (math-csc-raw (math-complex x))))
1150         (t
1151          (let ((sn (math-sin-raw x)))
1152            (if (eq sn 0)
1153                (math-div 1 0)
1154              (math-div-float '(float 1 0) sn))))))
1155
1156 (defun math-cot-raw (x)   ; [N N]
1157   (cond ((eq (car x) 'cplx)
1158          (let* ((x (math-mul x '(float 1 0)))
1159                 (expx (math-exp-raw (nth 2 x)))
1160                 (expmx (math-div-float '(float 1 0) expx))
1161                 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
1162                 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
1163                 (sc (math-sin-cos-raw (nth 1 x)))
1164                 (d (math-add-float
1165                     (math-sqr (car sc))
1166                     (math-sqr sh))))
1167            (and (not (eq (nth 1 d) 0))
1168                 (list 'cplx
1169                       (math-div-float
1170                        (math-mul-float (car sc) (cdr sc))
1171                        d)
1172                       (math-neg
1173                        (math-div-float
1174                         (math-mul-float sh ch)
1175                         d))))))
1176         ((eq (car x) 'polar)
1177          (math-polar (math-cot-raw (math-complex x))))
1178         (t
1179          (let ((sc (math-sin-cos-raw x)))
1180            (if (eq (nth 1 (car sc)) 0)
1181                (math-div (cdr sc) 0)
1182              (math-div-float (cdr sc) (car sc)))))))
1183
1184
1185 ;;; This could use a smarter method:  Reduce x as in math-sin-raw, then
1186 ;;;   compute either sin(x) or cos(x), whichever is smaller, and compute
1187 ;;;   the other using the identity sin(x)^2 + cos(x)^2 = 1.
1188 (defun math-sin-cos-raw (x)   ; [F.F F]  (result is (sin x . cos x))
1189   (cons (math-sin-raw x) (math-cos-raw x)))
1190
1191 (defun math-tan-raw (x)   ; [N N]
1192   (cond ((eq (car x) 'cplx)
1193          (let* ((x (math-mul x '(float 2 0)))
1194                 (expx (math-exp-raw (nth 2 x)))
1195                 (expmx (math-div-float '(float 1 0) expx))
1196                 (sc (math-sin-cos-raw (nth 1 x)))
1197                 (d (math-add-float (cdr sc)
1198                                    (math-mul-float (math-add-float expx expmx)
1199                                                    '(float 5 -1)))))
1200            (and (not (eq (nth 1 d) 0))
1201                 (list 'cplx
1202                       (math-div-float (car sc) d)
1203                       (math-div-float (math-mul-float (math-sub-float expx
1204                                                                       expmx)
1205                                                       '(float 5 -1)) d)))))
1206         ((eq (car x) 'polar)
1207          (math-polar (math-tan-raw (math-complex x))))
1208         (t
1209          (let ((sc (math-sin-cos-raw x)))
1210            (if (eq (nth 1 (cdr sc)) 0)
1211                (math-div (car sc) 0)
1212              (math-div-float (car sc) (cdr sc)))))))
1213
1214 (defun math-sin-raw-2 (x orgx)   ; This avoids poss of inf recursion.  [F F]
1215   (let ((xmpo2 (math-sub-float (math-pi-over-2) x)))
1216     (cond ((Math-integer-negp (nth 1 xmpo2))
1217            (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi))
1218                                            orgx)))
1219           ((math-lessp-float (math-pi-over-4) x)
1220            (math-cos-raw-2 xmpo2 orgx))
1221           ((math-lessp-float x (math-neg (math-pi-over-4)))
1222            (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx)))
1223           ((math-with-extra-prec -1 (math-nearly-zerop-float x orgx))
1224            '(float 0 0))
1225           ((math-use-emacs-fn 'sin x))
1226           (calc-symbolic-mode (signal 'inexact-result nil))
1227           (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x)))))))
1228
1229 (defun math-cos-raw-2 (x orgx)   ; [F F]
1230   (cond ((math-with-extra-prec -1 (math-nearly-zerop-float x orgx))
1231          '(float 1 0))
1232         ((math-use-emacs-fn 'cos x))
1233         (calc-symbolic-mode (signal 'inexact-result nil))
1234         (t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
1235              (math-sin-series
1236               (math-add-float '(float 1 0)
1237                               (math-mul-float xnegsqr '(float 5 -1)))
1238               24 5 xnegsqr xnegsqr)))))
1239
1240 (defun math-sin-series (sum nfac n x xnegsqr)
1241   (math-working "sin" sum)
1242   (let* ((nextx (math-mul-float x xnegsqr))
1243          (nextsum (math-add-float sum (math-div-float nextx
1244                                                       (math-float nfac)))))
1245     (if (math-nearly-equal-float sum nextsum)
1246         sum
1247       (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
1248                        (+ n 2) nextx xnegsqr))))
1249
1250
1251 ;;; Inverse sine, cosine, tangent.
1252
1253 (defun calcFunc-arcsin (x)   ; [N N] [Public]
1254   (cond ((eq x 0) 0)
1255         ((and (eq x 1) (eq calc-angle-mode 'deg)) 90)
1256         ((and (eq x -1) (eq calc-angle-mode 'deg)) -90)
1257         (calc-symbolic-mode (signal 'inexact-result nil))
1258         ((Math-numberp x)
1259          (math-with-extra-prec 2
1260            (math-from-radians (math-arcsin-raw (math-float x)))))
1261         ((eq (car x) 'sdev)
1262          (math-make-sdev (calcFunc-arcsin (nth 1 x))
1263                          (math-from-radians
1264                           (math-div (nth 2 x)
1265                                     (math-sqrt
1266                                      (math-sub 1 (math-sqr (nth 1 x))))))))
1267         ((eq (car x) 'intv)
1268          (math-sort-intv (nth 1 x)
1269                          (calcFunc-arcsin (nth 2 x))
1270                          (calcFunc-arcsin (nth 3 x))))
1271         ((equal x '(var nan var-nan))
1272          x)
1273         (t (calc-record-why 'numberp x)
1274            (list 'calcFunc-arcsin x))))
1275
1276 (defun calcFunc-arccos (x)   ; [N N] [Public]
1277   (cond ((eq x 1) 0)
1278         ((and (eq x 0) (eq calc-angle-mode 'deg)) 90)
1279         ((and (eq x -1) (eq calc-angle-mode 'deg)) 180)
1280         (calc-symbolic-mode (signal 'inexact-result nil))
1281         ((Math-numberp x)
1282          (math-with-extra-prec 2
1283            (math-from-radians (math-arccos-raw (math-float x)))))
1284         ((eq (car x) 'sdev)
1285          (math-make-sdev (calcFunc-arccos (nth 1 x))
1286                          (math-from-radians
1287                           (math-div (nth 2 x)
1288                                     (math-sqrt
1289                                      (math-sub 1 (math-sqr (nth 1 x))))))))
1290         ((eq (car x) 'intv)
1291          (math-sort-intv (nth 1 x)
1292                          (calcFunc-arccos (nth 2 x))
1293                          (calcFunc-arccos (nth 3 x))))
1294         ((equal x '(var nan var-nan))
1295          x)
1296         (t (calc-record-why 'numberp x)
1297            (list 'calcFunc-arccos x))))
1298
1299 (defun calcFunc-arctan (x)   ; [N N] [Public]
1300   (cond ((eq x 0) 0)
1301         ((and (eq x 1) (eq calc-angle-mode 'deg)) 45)
1302         ((and (eq x -1) (eq calc-angle-mode 'deg)) -45)
1303         ((Math-numberp x)
1304          (math-with-extra-prec 2
1305            (math-from-radians (math-arctan-raw (math-float x)))))
1306         ((eq (car x) 'sdev)
1307          (math-make-sdev (calcFunc-arctan (nth 1 x))
1308                          (math-from-radians
1309                           (math-div (nth 2 x)
1310                                     (math-add 1 (math-sqr (nth 1 x)))))))
1311         ((eq (car x) 'intv)
1312          (math-sort-intv (nth 1 x)
1313                          (calcFunc-arctan (nth 2 x))
1314                          (calcFunc-arctan (nth 3 x))))
1315         ((equal x '(var inf var-inf))
1316          (math-quarter-circle t))
1317         ((equal x '(neg (var inf var-inf)))
1318          (math-neg (math-quarter-circle t)))
1319         ((equal x '(var nan var-nan))
1320          x)
1321         (t (calc-record-why 'numberp x)
1322            (list 'calcFunc-arctan x))))
1323
1324 (defun math-arcsin-raw (x)   ; [N N]
1325   (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x)))))
1326     (if (or (memq (car x) '(cplx polar))
1327             (memq (car a) '(cplx polar)))
1328         (math-with-extra-prec 2   ; use extra precision for difficult case
1329           (math-mul '(cplx 0 -1)
1330                     (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a))))
1331       (math-arctan2-raw x a))))
1332
1333 (defun math-arccos-raw (x)   ; [N N]
1334   (math-sub (math-pi-over-2) (math-arcsin-raw x)))
1335
1336 (defun math-arctan-raw (x)   ; [N N]
1337   (cond ((memq (car x) '(cplx polar))
1338          (math-with-extra-prec 2   ; extra-extra
1339            (math-div (math-sub
1340                       (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x)))
1341                       (math-ln-raw (math-add 1 (math-mul '(cplx 0 -1) x))))
1342                      '(cplx 0 2))))
1343         ((Math-integer-negp (nth 1 x))
1344          (math-neg-float (math-arctan-raw (math-neg-float x))))
1345         ((math-zerop x) x)
1346         ((math-use-emacs-fn 'atan x))
1347         (calc-symbolic-mode (signal 'inexact-result nil))
1348         ((math-equal-int x 1) (math-pi-over-4))
1349         ((math-equal-int x -1) (math-neg (math-pi-over-4)))
1350         ((math-lessp-float '(float 414214 -6) x)  ; if x > sqrt(2) - 1, reduce
1351          (if (math-lessp-float '(float 1 0) x)
1352              (math-sub-float (math-mul-float (math-pi) '(float 5 -1))
1353                              (math-arctan-raw (math-div-float '(float 1 0) x)))
1354            (math-sub-float (math-mul-float (math-pi) '(float 25 -2))
1355                            (math-arctan-raw (math-div-float
1356                                              (math-sub-float '(float 1 0) x)
1357                                              (math-add-float '(float 1 0)
1358                                                              x))))))
1359         (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x))))))
1360
1361 (defun math-arctan-series (sum n x xnegsqr)
1362   (math-working "arctan" sum)
1363   (let* ((nextx (math-mul-float x xnegsqr))
1364          (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1365     (if (math-nearly-equal-float sum nextsum)
1366         sum
1367       (math-arctan-series nextsum (+ n 2) nextx xnegsqr))))
1368
1369 (defun calcFunc-arctan2 (y x)   ; [F R R] [Public]
1370   (if (Math-anglep y)
1371       (if (Math-anglep x)
1372           (math-with-extra-prec 2
1373             (math-from-radians (math-arctan2-raw (math-float y)
1374                                                  (math-float x))))
1375         (calc-record-why 'anglep x)
1376         (list 'calcFunc-arctan2 y x))
1377     (if (and (or (math-infinitep x) (math-anglep x))
1378              (or (math-infinitep y) (math-anglep y)))
1379         (progn
1380           (if (math-posp x)
1381               (setq x 1)
1382             (if (math-negp x)
1383                 (setq x -1)
1384               (or (math-zerop x)
1385                   (setq x nil))))
1386           (if (math-posp y)
1387               (setq y 1)
1388             (if (math-negp y)
1389                 (setq y -1)
1390               (or (math-zerop y)
1391                   (setq y nil))))
1392           (if (and y x)
1393               (calcFunc-arctan2 y x)
1394             '(var nan var-nan)))
1395       (calc-record-why 'anglep y)
1396       (list 'calcFunc-arctan2 y x))))
1397
1398 (defun math-arctan2-raw (y x)   ; [F R R]
1399   (cond ((math-zerop y)
1400          (if (math-negp x) (math-pi)
1401            (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0)))
1402         ((math-zerop x)
1403          (if (math-posp y)
1404              (math-pi-over-2)
1405            (math-neg (math-pi-over-2))))
1406         ((math-posp x)
1407          (math-arctan-raw (math-div-float y x)))
1408         ((math-posp y)
1409          (math-add-float (math-arctan-raw (math-div-float y x))
1410                          (math-pi)))
1411         (t
1412          (math-sub-float (math-arctan-raw (math-div-float y x))
1413                          (math-pi)))))
1414
1415 (defun calcFunc-arcsincos (x)   ; [V N] [Public]
1416   (if (and (Math-vectorp x)
1417            (= (length x) 3))
1418       (calcFunc-arctan2 (nth 2 x) (nth 1 x))
1419     (math-reject-arg x "*Two-element vector expected")))
1420
1421
1422
1423 ;;; Exponential function.
1424
1425 (defun calcFunc-exp (x)   ; [N N] [Public]
1426   (cond ((eq x 0) 1)
1427         ((and (memq x '(1 -1)) calc-symbolic-mode)
1428          (if (eq x 1) '(var e var-e) (math-div 1 '(var e var-e))))
1429         ((Math-numberp x)
1430          (math-with-extra-prec 2 (math-exp-raw (math-float x))))
1431         ((eq (car-safe x) 'sdev)
1432          (let ((ex (calcFunc-exp (nth 1 x))))
1433            (math-make-sdev ex (math-mul (nth 2 x) ex))))
1434         ((eq (car-safe x) 'intv)
1435          (math-make-intv (nth 1 x) (calcFunc-exp (nth 2 x))
1436                          (calcFunc-exp (nth 3 x))))
1437         ((equal x '(var inf var-inf))
1438          x)
1439         ((equal x '(neg (var inf var-inf)))
1440          0)
1441         ((equal x '(var nan var-nan))
1442          x)
1443         (t (calc-record-why 'numberp x)
1444            (list 'calcFunc-exp x))))
1445
1446 (defun calcFunc-expm1 (x)   ; [N N] [Public]
1447   (cond ((eq x 0) 0)
1448         ((math-zerop x) '(float 0 0))
1449         (calc-symbolic-mode (signal 'inexact-result nil))
1450         ((Math-numberp x)
1451          (math-with-extra-prec 2
1452            (let ((x (math-float x)))
1453              (if (and (eq (car x) 'float)
1454                       (math-lessp-float x '(float 1 0))
1455                       (math-lessp-float '(float -1 0) x))
1456                  (math-exp-minus-1-raw x)
1457                (math-add (math-exp-raw x) -1)))))
1458         ((eq (car-safe x) 'sdev)
1459          (if (math-constp x)
1460              (let ((ex (calcFunc-expm1 (nth 1 x))))
1461                (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1))))
1462            (math-make-sdev (calcFunc-expm1 (nth 1 x))
1463                            (math-mul (nth 2 x) (calcFunc-exp (nth 1 x))))))
1464         ((eq (car-safe x) 'intv)
1465          (math-make-intv (nth 1 x)
1466                          (calcFunc-expm1 (nth 2 x))
1467                          (calcFunc-expm1 (nth 3 x))))
1468         ((equal x '(var inf var-inf))
1469          x)
1470         ((equal x '(neg (var inf var-inf)))
1471          -1)
1472         ((equal x '(var nan var-nan))
1473          x)
1474         (t (calc-record-why 'numberp x)
1475            (list 'calcFunc-expm1 x))))
1476
1477 (defun calcFunc-exp10 (x)   ; [N N] [Public]
1478   (if (eq x 0)
1479       1
1480     (math-pow '(float 1 1) x)))
1481
1482 (defun math-exp-raw (x)   ; [N N]
1483   (cond ((math-zerop x) '(float 1 0))
1484         (calc-symbolic-mode (signal 'inexact-result nil))
1485         ((eq (car x) 'cplx)
1486          (let ((expx (math-exp-raw (nth 1 x)))
1487                (sc (math-sin-cos-raw (nth 2 x))))
1488            (list 'cplx
1489                  (math-mul-float expx (cdr sc))
1490                  (math-mul-float expx (car sc)))))
1491         ((eq (car x) 'polar)
1492          (let ((xc (math-complex x)))
1493            (list 'polar
1494                  (math-exp-raw (nth 1 xc))
1495                  (math-from-radians (nth 2 xc)))))
1496         ((math-use-emacs-fn 'exp x))
1497         ((or (math-lessp-float '(float 5 -1) x)
1498              (math-lessp-float x '(float -5 -1)))
1499          (if (math-lessp-float '(float 921035 1) x)
1500              (math-overflow)
1501            (if (math-lessp-float x '(float -921035 1))
1502                (math-underflow)))
1503          (let* ((two-x (math-mul-float x '(float 2 0)))
1504                 (hint (math-scale-int (nth 1 two-x) (nth 2 two-x)))
1505                 (hfrac (math-sub-float x (math-mul-float (math-float hint)
1506                                                          '(float 5 -1)))))
1507            (math-mul-float (math-ipow (math-sqrt-e) hint)
1508                            (math-add-float '(float 1 0)
1509                                            (math-exp-minus-1-raw hfrac)))))
1510         (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x)))))
1511
1512 (defun math-exp-minus-1-raw (x)   ; [F F]
1513   (math-exp-series x 2 3 x x))
1514
1515 (defun math-exp-series (sum nfac n xpow x)
1516   (math-working "exp" sum)
1517   (let* ((nextx (math-mul-float xpow x))
1518          (nextsum (math-add-float sum (math-div-float nextx
1519                                                       (math-float nfac)))))
1520     (if (math-nearly-equal-float sum nextsum)
1521         sum
1522       (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x))))
1523
1524
1525
1526 ;;; Logarithms.
1527
1528 (defun calcFunc-ln (x)   ; [N N] [Public]
1529   (cond ((math-zerop x)
1530          (if calc-infinite-mode
1531              '(neg (var inf var-inf))
1532            (math-reject-arg x "*Logarithm of zero")))
1533         ((eq x 1) 0)
1534         ((Math-numberp x)
1535          (math-with-extra-prec 2 (math-ln-raw (math-float x))))
1536         ((eq (car-safe x) 'sdev)
1537          (math-make-sdev (calcFunc-ln (nth 1 x))
1538                          (math-div (nth 2 x) (nth 1 x))))
1539         ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1540                                           (Math-zerop (nth 2 x))
1541                                           (not (math-intv-constp x))))
1542          (let ((calc-infinite-mode t))
1543            (math-make-intv (nth 1 x) (calcFunc-ln (nth 2 x))
1544                            (calcFunc-ln (nth 3 x)))))
1545         ((equal x '(var e var-e))
1546          1)
1547         ((and (eq (car-safe x) '^)
1548               (equal (nth 1 x) '(var e var-e))
1549               (math-known-realp (nth 2 x)))
1550          (nth 2 x))
1551         ((math-infinitep x)
1552          (if (equal x '(var nan var-nan))
1553              x
1554            '(var inf var-inf)))
1555         (t (calc-record-why 'numberp x)
1556            (list 'calcFunc-ln x))))
1557
1558 (defun calcFunc-log10 (x)   ; [N N] [Public]
1559   (cond ((math-equal-int x 1)
1560          (if (math-floatp x) '(float 0 0) 0))
1561         ((and (Math-integerp x)
1562               (math-posp x)
1563               (let ((res (math-integer-log x 10)))
1564                 (and (car res)
1565                      (setq x (cdr res)))))
1566          x)
1567         ((and (eq (car-safe x) 'frac)
1568               (eq (nth 1 x) 1)
1569               (let ((res (math-integer-log (nth 2 x) 10)))
1570                 (and (car res)
1571                      (setq x (- (cdr res))))))
1572          x)
1573         ((math-zerop x)
1574          (if calc-infinite-mode
1575              '(neg (var inf var-inf))
1576            (math-reject-arg x "*Logarithm of zero")))
1577         (calc-symbolic-mode (signal 'inexact-result nil))
1578         ((Math-numberp x)
1579          (math-with-extra-prec 2
1580            (let ((xf (math-float x)))
1581              (if (eq (nth 1 xf) 0)
1582                  (math-reject-arg x "*Logarithm of zero"))
1583              (if (Math-integer-posp (nth 1 xf))
1584                  (if (eq (nth 1 xf) 1)    ; log10(1*10^n) = n
1585                      (math-float (nth 2 xf))
1586                    (let ((xdigs (1- (math-numdigs (nth 1 xf)))))
1587                      (math-add-float
1588                       (math-div-float (math-ln-raw-2
1589                                        (list 'float (nth 1 xf) (- xdigs)))
1590                                       (math-ln-10))
1591                       (math-float (+ (nth 2 xf) xdigs)))))
1592                (math-div (calcFunc-ln xf) (math-ln-10))))))
1593         ((eq (car-safe x) 'sdev)
1594          (math-make-sdev (calcFunc-log10 (nth 1 x))
1595                          (math-div (nth 2 x)
1596                                    (math-mul (nth 1 x) (math-ln-10)))))
1597         ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1598                                           (not (math-intv-constp x))))
1599          (math-make-intv (nth 1 x)
1600                          (calcFunc-log10 (nth 2 x))
1601                          (calcFunc-log10 (nth 3 x))))
1602         ((math-infinitep x)
1603          (if (equal x '(var nan var-nan))
1604              x
1605            '(var inf var-inf)))
1606         (t (calc-record-why 'numberp x)
1607            (list 'calcFunc-log10 x))))
1608
1609 (defun calcFunc-log (x &optional b)   ; [N N N] [Public]
1610   (cond ((or (null b) (equal b '(var e var-e)))
1611          (math-normalize (list 'calcFunc-ln x)))
1612         ((or (eq b 10) (equal b '(float 1 1)))
1613          (math-normalize (list 'calcFunc-log10 x)))
1614         ((math-zerop x)
1615          (if calc-infinite-mode
1616              (math-div (calcFunc-ln x) (calcFunc-ln b))
1617            (math-reject-arg x "*Logarithm of zero")))
1618         ((math-zerop b)
1619          (if calc-infinite-mode
1620              (math-div (calcFunc-ln x) (calcFunc-ln b))
1621            (math-reject-arg b "*Logarithm of zero")))
1622         ((math-equal-int b 1)
1623          (if calc-infinite-mode
1624              (math-div (calcFunc-ln x) 0)
1625            (math-reject-arg b "*Logarithm base one")))
1626         ((math-equal-int x 1)
1627          (if (math-floatp b) '(float 0 0) 0))
1628         ((and (Math-ratp x) (Math-ratp b)
1629               (math-posp x) (math-posp b)
1630               (let* ((sign 1) (inv nil)
1631                      (xx (if (Math-lessp 1 x)
1632                              x
1633                            (setq sign -1)
1634                            (math-div 1 x)))
1635                      (bb (if (Math-lessp 1 b)
1636                              b
1637                            (setq sign (- sign))
1638                            (math-div 1 b)))
1639                      (res (if (Math-lessp xx bb)
1640                               (setq inv (math-integer-log bb xx))
1641                             (math-integer-log xx bb))))
1642                 (and (car res)
1643                      (setq x (if inv
1644                                  (math-div 1 (* sign (cdr res)))
1645                                (* sign (cdr res)))))))
1646          x)
1647         (calc-symbolic-mode (signal 'inexact-result nil))
1648         ((and (Math-numberp x) (Math-numberp b))
1649          (math-with-extra-prec 2
1650            (math-div (math-ln-raw (math-float x))
1651                      (math-log-base-raw b))))
1652         ((and (eq (car-safe x) 'sdev)
1653               (Math-numberp b))
1654          (math-make-sdev (calcFunc-log (nth 1 x) b)
1655                          (math-div (nth 2 x)
1656                                    (math-mul (nth 1 x)
1657                                              (math-log-base-raw b)))))
1658         ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1659                                           (not (math-intv-constp x)))
1660               (math-realp b))
1661          (math-make-intv (nth 1 x)
1662                          (calcFunc-log (nth 2 x) b)
1663                          (calcFunc-log (nth 3 x) b)))
1664         ((or (eq (car-safe x) 'intv) (eq (car-safe b) 'intv))
1665          (math-div (calcFunc-ln x) (calcFunc-ln b)))
1666         ((or (math-infinitep x)
1667              (math-infinitep b))
1668          (math-div (calcFunc-ln x) (calcFunc-ln b)))
1669         (t (if (Math-numberp b)
1670                (calc-record-why 'numberp x)
1671              (calc-record-why 'numberp b))
1672            (list 'calcFunc-log x b))))
1673
1674 (defun calcFunc-alog (x &optional b)
1675   (cond ((or (null b) (equal b '(var e var-e)))
1676          (math-normalize (list 'calcFunc-exp x)))
1677         (t (math-pow b x))))
1678
1679 (defun calcFunc-ilog (x b)
1680   (if (and (math-natnump x) (not (eq x 0))
1681            (math-natnump b) (not (eq b 0)))
1682       (if (eq b 1)
1683           (math-reject-arg x "*Logarithm base one")
1684         (if (Math-natnum-lessp x b)
1685             0
1686           (cdr (math-integer-log x b))))
1687     (math-floor (calcFunc-log x b))))
1688
1689 (defun math-integer-log (x b)
1690   (let ((pows (list b))
1691         (pow (math-sqr b))
1692         next
1693         sum n)
1694     (while (not (Math-lessp x pow))
1695       (setq pows (cons pow pows)
1696             pow (math-sqr pow)))
1697     (setq n (lsh 1 (1- (length pows)))
1698           sum n
1699           pow (car pows))
1700     (while (and (setq pows (cdr pows))
1701                 (Math-lessp pow x))
1702       (setq n (/ n 2)
1703             next (math-mul pow (car pows)))
1704       (or (Math-lessp x next)
1705           (setq pow next
1706                 sum (+ sum n))))
1707     (cons (equal pow x) sum)))
1708
1709
1710 (defvar math-log-base-cache nil)
1711 (defun math-log-base-raw (b)   ; [N N]
1712   (if (not (and (equal (car math-log-base-cache) b)
1713                 (eq (nth 1 math-log-base-cache) calc-internal-prec)))
1714       (setq math-log-base-cache (list b calc-internal-prec
1715                                       (math-ln-raw (math-float b)))))
1716   (nth 2 math-log-base-cache))
1717
1718 (defun calcFunc-lnp1 (x)   ; [N N] [Public]
1719   (cond ((Math-equal-int x -1)
1720          (if calc-infinite-mode
1721              '(neg (var inf var-inf))
1722            (math-reject-arg x "*Logarithm of zero")))
1723         ((eq x 0) 0)
1724         ((math-zerop x) '(float 0 0))
1725         (calc-symbolic-mode (signal 'inexact-result nil))
1726         ((Math-numberp x)
1727          (math-with-extra-prec 2
1728            (let ((x (math-float x)))
1729              (if (and (eq (car x) 'float)
1730                       (math-lessp-float x '(float 5 -1))
1731                       (math-lessp-float '(float -5 -1) x))
1732                  (math-ln-plus-1-raw x)
1733                (math-ln-raw (math-add-float x '(float 1 0)))))))
1734         ((eq (car-safe x) 'sdev)
1735          (math-make-sdev (calcFunc-lnp1 (nth 1 x))
1736                          (math-div (nth 2 x) (math-add (nth 1 x) 1))))
1737         ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1738                                           (not (math-intv-constp x))))
1739          (math-make-intv (nth 1 x)
1740                          (calcFunc-lnp1 (nth 2 x))
1741                          (calcFunc-lnp1 (nth 3 x))))
1742         ((math-infinitep x)
1743          (if (equal x '(var nan var-nan))
1744              x
1745            '(var inf var-inf)))
1746         (t (calc-record-why 'numberp x)
1747            (list 'calcFunc-lnp1 x))))
1748
1749 (defun math-ln-raw (x)    ; [N N] --- must be float format!
1750   (cond ((eq (car-safe x) 'cplx)
1751          (list 'cplx
1752                (math-mul-float (math-ln-raw
1753                                 (math-add-float (math-sqr-float (nth 1 x))
1754                                                 (math-sqr-float (nth 2 x))))
1755                                '(float 5 -1))
1756                (math-arctan2-raw (nth 2 x) (nth 1 x))))
1757         ((eq (car x) 'polar)
1758          (math-polar (list 'cplx
1759                            (math-ln-raw (nth 1 x))
1760                            (math-to-radians (nth 2 x)))))
1761         ((Math-equal-int x 1)
1762          '(float 0 0))
1763         (calc-symbolic-mode (signal 'inexact-result nil))
1764         ((math-posp (nth 1 x))    ; positive and real
1765          (cond
1766           ((math-use-emacs-fn 'log x))
1767           (t
1768            (let ((xdigs (1- (math-numdigs (nth 1 x)))))
1769              (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs)))
1770                              (math-mul-float (math-float (+ (nth 2 x) xdigs))
1771                                              (math-ln-10)))))))
1772         ((math-zerop x)
1773          (math-reject-arg x "*Logarithm of zero"))
1774         ((eq calc-complex-mode 'polar)    ; negative and real
1775          (math-polar
1776           (list 'cplx   ; negative and real
1777                 (math-ln-raw (math-neg-float x))
1778                 (math-pi))))
1779         (t (list 'cplx   ; negative and real
1780                  (math-ln-raw (math-neg-float x))
1781                  (math-pi)))))
1782
1783 (defun math-ln-raw-2 (x)    ; [F F]
1784   (cond ((math-lessp-float '(float 14 -1) x)
1785          (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1)))
1786                          (math-ln-2)))
1787         (t    ; now .7 < x <= 1.4
1788          (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
1789                                         (math-add-float x '(float 1 0)))))))
1790
1791 (defun math-ln-raw-3 (x)   ; [F F]
1792   (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
1793                   '(float 2 0)))
1794
1795 ;;; Compute ln((1+x)/(1-x))
1796 (defun math-ln-raw-series (sum n x xsqr)
1797   (math-working "log" sum)
1798   (let* ((nextx (math-mul-float x xsqr))
1799          (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1800     (if (math-nearly-equal-float sum nextsum)
1801         sum
1802       (math-ln-raw-series nextsum (+ n 2) nextx xsqr))))
1803
1804 (defun math-ln-plus-1-raw (x)
1805   (math-lnp1-series x 2 x (math-neg x)))
1806
1807 (defun math-lnp1-series (sum n xpow x)
1808   (math-working "lnp1" sum)
1809   (let* ((nextx (math-mul-float xpow x))
1810          (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1811     (if (math-nearly-equal-float sum nextsum)
1812         sum
1813       (math-lnp1-series nextsum (1+ n) nextx x))))
1814
1815 (defconst math-approx-ln-10
1816   (math-read-number-simple "2.302585092994045684018")
1817   "An approximation for ln(10).")
1818
1819 (math-defcache math-ln-10 math-approx-ln-10
1820   (math-ln-raw-2 '(float 1 1)))
1821
1822 (defconst math-approx-ln-2
1823   (math-read-number-simple "0.693147180559945309417")
1824   "An approximation for ln(2).")
1825
1826 (math-defcache math-ln-2 math-approx-ln-2
1827   (math-ln-raw-3 (math-float '(frac 1 3))))
1828
1829
1830
1831 ;;; Hyperbolic functions.
1832
1833 (defun calcFunc-sinh (x)   ; [N N] [Public]
1834   (cond ((eq x 0) 0)
1835         (math-expand-formulas
1836          (math-normalize
1837           (list '/ (list '- (list 'calcFunc-exp x)
1838                          (list 'calcFunc-exp (list 'neg x))) 2)))
1839         ((Math-numberp x)
1840          (if calc-symbolic-mode (signal 'inexact-result nil))
1841          (math-with-extra-prec 2
1842            (let ((expx (math-exp-raw (math-float x))))
1843              (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1)))))
1844         ((eq (car-safe x) 'sdev)
1845          (math-make-sdev (calcFunc-sinh (nth 1 x))
1846                          (math-mul (nth 2 x) (calcFunc-cosh (nth 1 x)))))
1847         ((eq (car x) 'intv)
1848          (math-sort-intv (nth 1 x)
1849                          (calcFunc-sinh (nth 2 x))
1850                          (calcFunc-sinh (nth 3 x))))
1851         ((or (equal x '(var inf var-inf))
1852              (equal x '(neg (var inf var-inf)))
1853              (equal x '(var nan var-nan)))
1854          x)
1855         (t (calc-record-why 'numberp x)
1856            (list 'calcFunc-sinh x))))
1857 (put 'calcFunc-sinh 'math-expandable t)
1858
1859 (defun calcFunc-cosh (x)   ; [N N] [Public]
1860   (cond ((eq x 0) 1)
1861         (math-expand-formulas
1862          (math-normalize
1863           (list '/ (list '+ (list 'calcFunc-exp x)
1864                          (list 'calcFunc-exp (list 'neg x))) 2)))
1865         ((Math-numberp x)
1866          (if calc-symbolic-mode (signal 'inexact-result nil))
1867          (math-with-extra-prec 2
1868            (let ((expx (math-exp-raw (math-float x))))
1869              (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1)))))
1870         ((eq (car-safe x) 'sdev)
1871          (math-make-sdev (calcFunc-cosh (nth 1 x))
1872                          (math-mul (nth 2 x)
1873                                    (calcFunc-sinh (nth 1 x)))))
1874         ((and (eq (car x) 'intv) (math-intv-constp x))
1875          (setq x (math-abs x))
1876          (math-sort-intv (nth 1 x)
1877                          (calcFunc-cosh (nth 2 x))
1878                          (calcFunc-cosh (nth 3 x))))
1879         ((or (equal x '(var inf var-inf))
1880              (equal x '(neg (var inf var-inf)))
1881              (equal x '(var nan var-nan)))
1882          (math-abs x))
1883         (t (calc-record-why 'numberp x)
1884            (list 'calcFunc-cosh x))))
1885 (put 'calcFunc-cosh 'math-expandable t)
1886
1887 (defun calcFunc-tanh (x)   ; [N N] [Public]
1888   (cond ((eq x 0) 0)
1889         (math-expand-formulas
1890          (math-normalize
1891           (let ((expx (list 'calcFunc-exp x))
1892                 (expmx (list 'calcFunc-exp (list 'neg x))))
1893             (math-normalize
1894              (list '/ (list '- expx expmx) (list '+ expx expmx))))))
1895         ((Math-numberp x)
1896          (if calc-symbolic-mode (signal 'inexact-result nil))
1897          (math-with-extra-prec 2
1898            (let* ((expx (calcFunc-exp (math-float x)))
1899                   (expmx (math-div 1 expx)))
1900              (math-div (math-sub expx expmx)
1901                        (math-add expx expmx)))))
1902         ((eq (car-safe x) 'sdev)
1903          (math-make-sdev (calcFunc-tanh (nth 1 x))
1904                          (math-div (nth 2 x)
1905                                    (math-sqr (calcFunc-cosh (nth 1 x))))))
1906         ((eq (car x) 'intv)
1907          (math-sort-intv (nth 1 x)
1908                          (calcFunc-tanh (nth 2 x))
1909                          (calcFunc-tanh (nth 3 x))))
1910         ((equal x '(var inf var-inf))
1911          1)
1912         ((equal x '(neg (var inf var-inf)))
1913          -1)
1914         ((equal x '(var nan var-nan))
1915          x)
1916         (t (calc-record-why 'numberp x)
1917            (list 'calcFunc-tanh x))))
1918 (put 'calcFunc-tanh 'math-expandable t)
1919
1920 (defun calcFunc-sech (x)   ; [N N] [Public]
1921   (cond ((eq x 0) 1)
1922         (math-expand-formulas
1923          (math-normalize
1924           (list '/ 2 (list '+ (list 'calcFunc-exp x)
1925                            (list 'calcFunc-exp (list 'neg x))))))
1926         ((Math-numberp x)
1927          (if calc-symbolic-mode (signal 'inexact-result nil))
1928          (math-with-extra-prec 2
1929            (let ((expx (math-exp-raw (math-float x))))
1930              (math-div '(float 2 0) (math-add expx (math-div 1 expx))))))
1931         ((eq (car-safe x) 'sdev)
1932          (math-make-sdev (calcFunc-sech (nth 1 x))
1933                          (math-mul (nth 2 x)
1934                                    (math-mul (calcFunc-sech (nth 1 x))
1935                                              (calcFunc-tanh (nth 1 x))))))
1936         ((and (eq (car x) 'intv) (math-intv-constp x))
1937          (setq x (math-abs x))
1938          (math-sort-intv (nth 1 x)
1939                          (calcFunc-sech (nth 2 x))
1940                          (calcFunc-sech (nth 3 x))))
1941         ((or (equal x '(var inf var-inf))
1942              (equal x '(neg (var inf var-inf))))
1943          0)
1944         ((equal x '(var nan var-nan))
1945          x)
1946         (t (calc-record-why 'numberp x)
1947            (list 'calcFunc-sech x))))
1948 (put 'calcFunc-sech 'math-expandable t)
1949
1950 (defun calcFunc-csch (x)   ; [N N] [Public]
1951   (cond ((eq x 0) (math-div 1 0))
1952         (math-expand-formulas
1953          (math-normalize
1954           (list '/ 2 (list '- (list 'calcFunc-exp x)
1955                            (list 'calcFunc-exp (list 'neg x))))))
1956         ((Math-numberp x)
1957          (if calc-symbolic-mode (signal 'inexact-result nil))
1958          (math-with-extra-prec 2
1959            (let ((expx (math-exp-raw (math-float x))))
1960              (math-div '(float 2 0) (math-add expx (math-div -1 expx))))))
1961         ((eq (car-safe x) 'sdev)
1962          (math-make-sdev (calcFunc-csch (nth 1 x))
1963                          (math-mul (nth 2 x)
1964                                    (math-mul (calcFunc-csch (nth 1 x))
1965                                              (calcFunc-coth (nth 1 x))))))
1966         ((eq (car x) 'intv)
1967          (if (and (Math-negp (nth 2 x))
1968                   (Math-posp (nth 3 x)))
1969              '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1970            (math-sort-intv (nth 1 x)
1971                            (calcFunc-csch (nth 2 x))
1972                            (calcFunc-csch (nth 3 x)))))
1973         ((or (equal x '(var inf var-inf))
1974              (equal x '(neg (var inf var-inf))))
1975          0)
1976         ((equal x '(var nan var-nan))
1977          x)
1978         (t (calc-record-why 'numberp x)
1979            (list 'calcFunc-csch x))))
1980 (put 'calcFunc-csch 'math-expandable t)
1981
1982 (defun calcFunc-coth (x)   ; [N N] [Public]
1983   (cond ((eq x 0) (math-div 1 0))
1984         (math-expand-formulas
1985          (math-normalize
1986           (let ((expx (list 'calcFunc-exp x))
1987                 (expmx (list 'calcFunc-exp (list 'neg x))))
1988             (math-normalize
1989              (list '/ (list '+ expx expmx) (list '- expx expmx))))))
1990         ((Math-numberp x)
1991          (if calc-symbolic-mode (signal 'inexact-result nil))
1992          (math-with-extra-prec 2
1993            (let* ((expx (calcFunc-exp (math-float x)))
1994                   (expmx (math-div 1 expx)))
1995              (math-div (math-add expx expmx)
1996                        (math-sub expx expmx)))))
1997         ((eq (car-safe x) 'sdev)
1998          (math-make-sdev (calcFunc-coth (nth 1 x))
1999                          (math-div (nth 2 x)
2000                                    (math-sqr (calcFunc-sinh (nth 1 x))))))
2001         ((eq (car x) 'intv)
2002          (if (and (Math-negp (nth 2 x))
2003                   (Math-posp (nth 3 x)))
2004              '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
2005            (math-sort-intv (nth 1 x)
2006                            (calcFunc-coth (nth 2 x))
2007                            (calcFunc-coth (nth 3 x)))))
2008         ((equal x '(var inf var-inf))
2009          1)
2010         ((equal x '(neg (var inf var-inf)))
2011          -1)
2012         ((equal x '(var nan var-nan))
2013          x)
2014         (t (calc-record-why 'numberp x)
2015            (list 'calcFunc-coth x))))
2016 (put 'calcFunc-coth 'math-expandable t)
2017
2018 (defun calcFunc-arcsinh (x)   ; [N N] [Public]
2019   (cond ((eq x 0) 0)
2020         (math-expand-formulas
2021          (math-normalize
2022           (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
2023                                               (list '+ (list '^ x 2) 1))))))
2024         ((Math-numberp x)
2025          (if calc-symbolic-mode (signal 'inexact-result nil))
2026          (math-with-extra-prec 2
2027            (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x)
2028                                                              '(float 1 0)))))))
2029         ((eq (car-safe x) 'sdev)
2030          (math-make-sdev (calcFunc-arcsinh (nth 1 x))
2031                          (math-div (nth 2 x)
2032                                    (math-sqrt
2033                                     (math-add (math-sqr (nth 1 x)) 1)))))
2034         ((eq (car x) 'intv)
2035          (math-sort-intv (nth 1 x)
2036                          (calcFunc-arcsinh (nth 2 x))
2037                          (calcFunc-arcsinh (nth 3 x))))
2038         ((or (equal x '(var inf var-inf))
2039              (equal x '(neg (var inf var-inf)))
2040              (equal x '(var nan var-nan)))
2041          x)
2042         (t (calc-record-why 'numberp x)
2043            (list 'calcFunc-arcsinh x))))
2044 (put 'calcFunc-arcsinh 'math-expandable t)
2045
2046 (defun calcFunc-arccosh (x)   ; [N N] [Public]
2047   (cond ((eq x 1) 0)
2048         ((and (eq x -1) calc-symbolic-mode)
2049          '(var pi var-pi))
2050         ((and (eq x 0) calc-symbolic-mode)
2051          (math-div (math-mul '(var pi var-pi) '(var i var-i)) 2))
2052         (math-expand-formulas
2053          (math-normalize
2054           (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
2055                                               (list '- (list '^ x 2) 1))))))
2056         ((Math-numberp x)
2057          (if calc-symbolic-mode (signal 'inexact-result nil))
2058          (if (Math-equal-int x -1)
2059              (math-imaginary (math-pi))
2060            (math-with-extra-prec 2
2061              (if (or t    ; need to do this even in the real case!
2062                      (memq (car-safe x) '(cplx polar)))
2063                  (let ((xp1 (math-add 1 x)))  ; this gets the branch cuts right
2064                    (math-ln-raw
2065                     (math-add x (math-mul xp1
2066                                           (math-sqrt-raw
2067                                            (math-div (math-sub
2068                                                       x
2069                                                       '(float 1 0))
2070                                                      xp1))))))
2071                (math-ln-raw
2072                 (math-add x (math-sqrt-raw (math-add (math-sqr x)
2073                                                      '(float -1 0)))))))))
2074         ((eq (car-safe x) 'sdev)
2075          (math-make-sdev (calcFunc-arccosh (nth 1 x))
2076                          (math-div (nth 2 x)
2077                                    (math-sqrt
2078                                     (math-add (math-sqr (nth 1 x)) -1)))))
2079         ((eq (car x) 'intv)
2080          (math-sort-intv (nth 1 x)
2081                          (calcFunc-arccosh (nth 2 x))
2082                          (calcFunc-arccosh (nth 3 x))))
2083         ((or (equal x '(var inf var-inf))
2084              (equal x '(neg (var inf var-inf)))
2085              (equal x '(var nan var-nan)))
2086          x)
2087         (t (calc-record-why 'numberp x)
2088            (list 'calcFunc-arccosh x))))
2089 (put 'calcFunc-arccosh 'math-expandable t)
2090
2091 (defun calcFunc-arctanh (x)   ; [N N] [Public]
2092   (cond ((eq x 0) 0)
2093         ((and (Math-equal-int x 1) calc-infinite-mode)
2094          '(var inf var-inf))
2095         ((and (Math-equal-int x -1) calc-infinite-mode)
2096          '(neg (var inf var-inf)))
2097         (math-expand-formulas
2098          (list '/ (list '-
2099                         (list 'calcFunc-ln (list '+ 1 x))
2100                         (list 'calcFunc-ln (list '- 1 x))) 2))
2101         ((Math-numberp x)
2102          (if calc-symbolic-mode (signal 'inexact-result nil))
2103          (math-with-extra-prec 2
2104            (if (or (memq (car-safe x) '(cplx polar))
2105                    (Math-lessp 1 x))
2106                (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x))
2107                                    (math-ln-raw (math-sub '(float 1 0) x)))
2108                          '(float 5 -1))
2109              (if (and (math-equal-int x 1) calc-infinite-mode)
2110                  '(var inf var-inf)
2111                (if (and (math-equal-int x -1) calc-infinite-mode)
2112                    '(neg (var inf var-inf))
2113                  (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x)
2114                                                   (math-sub 1 x)))
2115                            '(float 5 -1)))))))
2116         ((eq (car-safe x) 'sdev)
2117          (math-make-sdev (calcFunc-arctanh (nth 1 x))
2118                          (math-div (nth 2 x)
2119                                    (math-sub 1 (math-sqr (nth 1 x))))))
2120         ((eq (car x) 'intv)
2121          (math-sort-intv (nth 1 x)
2122                          (calcFunc-arctanh (nth 2 x))
2123                          (calcFunc-arctanh (nth 3 x))))
2124         ((equal x '(var nan var-nan))
2125          x)
2126         (t (calc-record-why 'numberp x)
2127            (list 'calcFunc-arctanh x))))
2128 (put 'calcFunc-arctanh 'math-expandable t)
2129
2130
2131 ;;; Convert A from HMS or degrees to radians.
2132 (defun calcFunc-rad (a)   ; [R R] [Public]
2133   (cond ((or (Math-numberp a)
2134              (eq (car a) 'intv))
2135          (math-with-extra-prec 2
2136            (math-mul a (math-pi-over-180))))
2137         ((eq (car a) 'hms)
2138          (math-from-hms a 'rad))
2139         ((eq (car a) 'sdev)
2140          (math-make-sdev (calcFunc-rad (nth 1 a))
2141                          (calcFunc-rad (nth 2 a))))
2142         (math-expand-formulas
2143          (math-div (math-mul a '(var pi var-pi)) 180))
2144         ((math-infinitep a) a)
2145         (t (list 'calcFunc-rad a))))
2146 (put 'calcFunc-rad 'math-expandable t)
2147
2148 ;;; Convert A from HMS or radians to degrees.
2149 (defun calcFunc-deg (a)   ; [R R] [Public]
2150   (cond ((or (Math-numberp a)
2151              (eq (car a) 'intv))
2152          (math-with-extra-prec 2
2153            (math-div a (math-pi-over-180))))
2154         ((eq (car a) 'hms)
2155          (math-from-hms a 'deg))
2156         ((eq (car a) 'sdev)
2157          (math-make-sdev (calcFunc-deg (nth 1 a))
2158                          (calcFunc-deg (nth 2 a))))
2159         (math-expand-formulas
2160          (math-div (math-mul 180 a) '(var pi var-pi)))
2161         ((math-infinitep a) a)
2162         (t (list 'calcFunc-deg a))))
2163 (put 'calcFunc-deg 'math-expandable t)
2164
2165 (provide 'calc-math)
2166
2167 ;;; calc-math.el ends here