1 ;;; calc-math.el --- mathematical functions for Calc
3 ;; Copyright (C) 1990, 1991, 1992, 1993, 2001, 2002, 2003, 2004,
4 ;; 2005, 2006, 2007 Free Software Foundation, Inc.
6 ;; Author: David Gillespie <daveg@synaptics.com>
7 ;; Maintainer: Jay Belanger <jay.p.belanger@gmail.com>
9 ;; This file is part of GNU Emacs.
11 ;; GNU Emacs is free software; you can redistribute it and/or modify
12 ;; it under the terms of the GNU General Public License as published by
13 ;; the Free Software Foundation; either version 2, or (at your option)
16 ;; GNU Emacs is distributed in the hope that it will be useful,
17 ;; but WITHOUT ANY WARRANTY; without even the implied warranty of
18 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 ;; GNU General Public License for more details.
21 ;; You should have received a copy of the GNU General Public License
22 ;; along with GNU Emacs; see the file COPYING. If not, write to the
23 ;; Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
24 ;; Boston, MA 02110-1301, USA.
30 ;; This file is autoloaded from calc-ext.el.
35 (defun calc-sqrt (arg)
39 (calc-unary-op "^2" 'calcFunc-sqr arg)
40 (calc-unary-op "sqrt" 'calcFunc-sqrt arg))))
42 (defun calc-isqrt (arg)
46 (calc-unary-op "^2" 'calcFunc-sqr arg)
47 (calc-unary-op "isqt" 'calcFunc-isqrt arg))))
50 (defun calc-hypot (arg)
53 (calc-binary-op "hypt" 'calcFunc-hypot arg)))
60 (defun calc-log10 (arg)
62 (calc-hyperbolic-func)
69 (calc-binary-op "alog" 'calcFunc-alog arg)
70 (calc-binary-op "log" 'calcFunc-log arg))))
72 (defun calc-ilog (arg)
76 (calc-binary-op "alog" 'calcFunc-alog arg)
77 (calc-binary-op "ilog" 'calcFunc-ilog arg))))
79 (defun calc-lnp1 (arg)
87 (if (calc-is-hyperbolic)
89 (calc-unary-op "lg10" 'calcFunc-log10 arg)
90 (calc-unary-op "10^" 'calcFunc-exp10 arg))
92 (calc-unary-op "ln" 'calcFunc-ln arg)
93 (calc-unary-op "exp" 'calcFunc-exp arg)))))
95 (defun calc-expm1 (arg)
99 (calc-unary-op "ln+1" 'calcFunc-lnp1 arg)
100 (calc-unary-op "ex-1" 'calcFunc-expm1 arg))))
105 (if (calc-is-inverse)
106 (if (calc-is-hyperbolic)
107 (if calc-symbolic-mode
108 (calc-pop-push-record 0 "phi" '(var phi var-phi))
109 (calc-pop-push-record 0 "phi" (math-phi)))
110 (if calc-symbolic-mode
111 (calc-pop-push-record 0 "gmma" '(var gamma var-gamma))
112 (calc-pop-push-record 0 "gmma" (math-gamma-const))))
113 (if (calc-is-hyperbolic)
114 (if calc-symbolic-mode
115 (calc-pop-push-record 0 "e" '(var e var-e))
116 (calc-pop-push-record 0 "e" (math-e)))
117 (if calc-symbolic-mode
118 (calc-pop-push-record 0 "pi" '(var pi var-pi))
119 (calc-pop-push-record 0 "pi" (math-pi)))))))
121 (defun calc-sin (arg)
124 (if (calc-is-hyperbolic)
125 (if (calc-is-inverse)
126 (calc-unary-op "asnh" 'calcFunc-arcsinh arg)
127 (calc-unary-op "sinh" 'calcFunc-sinh arg))
128 (if (calc-is-inverse)
129 (calc-unary-op "asin" 'calcFunc-arcsin arg)
130 (calc-unary-op "sin" 'calcFunc-sin arg)))))
132 (defun calc-arcsin (arg)
137 (defun calc-sinh (arg)
139 (calc-hyperbolic-func)
142 (defun calc-arcsinh (arg)
145 (calc-hyperbolic-func)
148 (defun calc-sec (arg)
151 (if (calc-is-hyperbolic)
152 (calc-unary-op "sech" 'calcFunc-sech arg)
153 (calc-unary-op "sec" 'calcFunc-sec arg))))
155 (defun calc-sech (arg)
157 (calc-hyperbolic-func)
160 (defun calc-cos (arg)
163 (if (calc-is-hyperbolic)
164 (if (calc-is-inverse)
165 (calc-unary-op "acsh" 'calcFunc-arccosh arg)
166 (calc-unary-op "cosh" 'calcFunc-cosh arg))
167 (if (calc-is-inverse)
168 (calc-unary-op "acos" 'calcFunc-arccos arg)
169 (calc-unary-op "cos" 'calcFunc-cos arg)))))
171 (defun calc-arccos (arg)
176 (defun calc-cosh (arg)
178 (calc-hyperbolic-func)
181 (defun calc-arccosh (arg)
184 (calc-hyperbolic-func)
187 (defun calc-csc (arg)
190 (if (calc-is-hyperbolic)
191 (calc-unary-op "csch" 'calcFunc-csch arg)
192 (calc-unary-op "csc" 'calcFunc-csc arg))))
194 (defun calc-csch (arg)
196 (calc-hyperbolic-func)
199 (defun calc-sincos ()
202 (if (calc-is-inverse)
203 (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1)))
204 (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1))))))
206 (defun calc-tan (arg)
209 (if (calc-is-hyperbolic)
210 (if (calc-is-inverse)
211 (calc-unary-op "atnh" 'calcFunc-arctanh arg)
212 (calc-unary-op "tanh" 'calcFunc-tanh arg))
213 (if (calc-is-inverse)
214 (calc-unary-op "atan" 'calcFunc-arctan arg)
215 (calc-unary-op "tan" 'calcFunc-tan arg)))))
217 (defun calc-arctan (arg)
222 (defun calc-tanh (arg)
224 (calc-hyperbolic-func)
227 (defun calc-arctanh (arg)
230 (calc-hyperbolic-func)
233 (defun calc-cot (arg)
236 (if (calc-is-hyperbolic)
237 (calc-unary-op "coth" 'calcFunc-coth arg)
238 (calc-unary-op "cot" 'calcFunc-cot arg))))
240 (defun calc-coth (arg)
242 (calc-hyperbolic-func)
245 (defun calc-arctan2 ()
248 (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2)))))
250 (defun calc-conj (arg)
253 (calc-unary-op "conj" 'calcFunc-conj arg)))
255 (defun calc-imaginary ()
258 (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))))
260 (defun calc-to-degrees (arg)
263 (calc-unary-op ">deg" 'calcFunc-deg arg)))
265 (defun calc-to-radians (arg)
268 (calc-unary-op ">rad" 'calcFunc-rad arg)))
271 (defun calc-degrees-mode (arg)
275 (calc-change-mode 'calc-angle-mode 'deg)
276 (message "Angles measured in degrees")))
277 ((= arg 2) (calc-radians-mode))
278 ((= arg 3) (calc-hms-mode))
279 (t (error "Prefix argument out of range"))))
281 (defun calc-radians-mode ()
284 (calc-change-mode 'calc-angle-mode 'rad)
285 (message "Angles measured in radians")))
288 ;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public]
289 ;;; This method takes advantage of the fact that Newton's method starting
290 ;;; with an overestimate always works, even using truncating integer division!
291 (defun math-isqrt (a)
292 (cond ((Math-zerop a) a)
293 ((not (math-natnump a))
294 (math-reject-arg a 'natnump))
296 (math-isqrt-small a))
298 (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a))))))))
300 (defun calcFunc-isqrt (a)
302 (math-isqrt (math-floor a))
303 (math-floor (math-sqrt a))))
306 ;;; This returns (flag . result) where the flag is t if A is a perfect square.
307 (defun math-isqrt-bignum (a) ; [P.l L]
308 (let ((len (length a)))
310 (let* ((top (nthcdr (- len 2) a)))
311 (math-isqrt-bignum-iter
313 (math-scale-bignum-digit-size
315 (1+ (math-isqrt-small
316 (+ (* (nth 1 top) math-bignum-digit-size) (car top)))))
318 (let* ((top (nth (1- len) a)))
319 (math-isqrt-bignum-iter
321 (math-scale-bignum-digit-size
322 (list (1+ (math-isqrt-small top)))
325 (defun math-isqrt-bignum-iter (a guess) ; [l L l]
326 (math-working "isqrt" (cons 'bigpos guess))
327 (let* ((q (math-div-bignum a guess))
328 (s (math-add-bignum (car q) guess))
329 (g2 (math-div2-bignum s))
330 (comp (math-compare-bignum g2 guess)))
332 (math-isqrt-bignum-iter a g2)
333 (cons (and (= comp 0)
334 (math-zerop-bignum (cdr q))
338 (defun math-zerop-bignum (a)
341 (while (eq (car (setq a (cdr a))) 0))
344 (defun math-scale-bignum-digit-size (a n) ; [L L S]
350 (defun math-isqrt-small (a) ; A > 0. [S S]
351 (let ((g (cond ((>= a 1000000) 10000)
356 (while (< (setq g2 (/ (+ g (/ a g)) 2)) g)
363 ;;; Compute the square root of a number.
364 ;;; [T N] if possible, else [F N] if possible, else [C N]. [Public]
367 (and (Math-zerop a) a)
368 (and (math-known-nonposp a)
369 (math-imaginary (math-sqrt (math-neg a))))
371 (let ((sqrt (math-isqrt-small a)))
372 (if (= (* sqrt sqrt) a)
374 (if calc-symbolic-mode
375 (list 'calcFunc-sqrt a)
376 (math-sqrt-float (math-float a) (math-float sqrt))))))
377 (and (eq (car-safe a) 'bigpos)
378 (let* ((res (math-isqrt-bignum (cdr a)))
379 (sqrt (math-normalize (cons 'bigpos (cdr res)))))
382 (if calc-symbolic-mode
383 (list 'calcFunc-sqrt a)
384 (math-sqrt-float (math-float a) (math-float sqrt))))))
385 (and (eq (car-safe a) 'frac)
386 (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a)))))
387 (num-sqrt (math-normalize (cons 'bigpos (cdr num-res))))
388 (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a)))))
389 (den-sqrt (math-normalize (cons 'bigpos (cdr den-res)))))
390 (if (and (car num-res) (car den-res))
391 (list 'frac num-sqrt den-sqrt)
392 (if calc-symbolic-mode
393 (if (or (car num-res) (car den-res))
394 (math-div (if (car num-res)
395 num-sqrt (list 'calcFunc-sqrt (nth 1 a)))
397 den-sqrt (list 'calcFunc-sqrt (nth 2 a))))
398 (list 'calcFunc-sqrt a))
399 (math-sqrt-float (math-float a)
400 (math-div (math-float num-sqrt) den-sqrt))))))
401 (and (eq (car-safe a) 'float)
402 (if calc-symbolic-mode
403 (if (= (% (nth 2 a) 2) 0)
404 (let ((res (math-isqrt-bignum
405 (cdr (Math-bignum-test (nth 1 a))))))
407 (math-make-float (math-normalize
408 (cons 'bigpos (cdr res)))
410 (signal 'inexact-result nil)))
411 (signal 'inexact-result nil))
412 (math-sqrt-float a)))
413 (and (eq (car-safe a) 'cplx)
414 (math-with-extra-prec 2
415 (let* ((d (math-abs a))
416 (imag (math-sqrt (math-mul (math-sub d (nth 1 a))
419 (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1)))
420 (if (math-negp (nth 2 a)) (math-neg imag) imag)))))
421 (and (eq (car-safe a) 'polar)
423 (math-sqrt (nth 1 a))
424 (math-mul (nth 2 a) '(float 5 -1))))
425 (and (eq (car-safe a) 'sdev)
426 (let ((sqrt (math-sqrt (nth 1 a))))
428 (math-div (nth 2 a) (math-mul sqrt 2)))))
429 (and (eq (car-safe a) 'intv)
430 (not (math-negp (nth 2 a)))
431 (math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a))))
432 (and (eq (car-safe a) '*)
433 (or (math-known-nonnegp (nth 1 a))
434 (math-known-nonnegp (nth 2 a)))
435 (math-mul (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
436 (and (eq (car-safe a) '/)
437 (or (and (math-known-nonnegp (nth 2 a))
438 (math-div (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
439 (and (math-known-nonnegp (nth 1 a))
440 (not (math-equal-int (nth 1 a) 1))
441 (math-mul (math-sqrt (nth 1 a))
442 (math-sqrt (math-div 1 (nth 2 a)))))))
443 (and (eq (car-safe a) '^)
444 (math-known-evenp (nth 2 a))
445 (math-known-realp (nth 1 a))
446 (math-abs (math-pow (nth 1 a) (math-div (nth 2 a) 2))))
447 (let ((inf (math-infinitep a)))
449 (math-mul (math-sqrt (math-infinite-dir a inf)) inf)))
451 (calc-record-why 'numberp a)
452 (list 'calcFunc-sqrt a))))
453 (defalias 'calcFunc-sqrt 'math-sqrt)
455 (defun math-infinite-dir (a &optional inf)
456 (or inf (setq inf (math-infinitep a)))
457 (math-normalize (math-expr-subst a inf 1)))
459 (defun math-sqrt-float (a &optional guess) ; [F F F]
460 (if calc-symbolic-mode
461 (signal 'inexact-result nil)
462 (math-with-extra-prec 1 (math-sqrt-raw a guess))))
464 (defun math-sqrt-raw (a &optional guess) ; [F F F]
465 (if (not (Math-posp a))
468 (let ((ldiff (- (math-numdigs (nth 1 a)) 6)))
469 (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff)))
470 (setq guess (math-make-float (math-isqrt-small
471 (math-scale-int (nth 1 a) (- ldiff)))
472 (/ (+ (nth 2 a) ldiff) 2)))))
473 (math-sqrt-float-iter a guess)))
475 (defun math-sqrt-float-iter (a guess) ; [F F F]
476 (math-working "sqrt" guess)
477 (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess))
479 (if (math-nearly-equal-float g2 guess)
481 (math-sqrt-float-iter a g2))))
483 ;;; True if A and B differ only in the last digit of precision. [P F F]
484 (defun math-nearly-equal-float (a b)
485 (let ((ediff (- (nth 2 a) (nth 2 b))))
486 (cond ((= ediff 0) ;; Expanded out for speed
487 (setq ediff (math-add (Math-integer-neg (nth 1 a)) (nth 1 b)))
489 (and (not (consp ediff))
492 (= (math-numdigs (nth 1 a)) calc-internal-prec))))
494 (setq ediff (math-add (Math-integer-neg (nth 1 b))
495 (math-scale-int (nth 1 a) 1)))
496 (and (not (consp ediff))
499 (= (math-numdigs (nth 1 b)) calc-internal-prec)))
501 (setq ediff (math-add (Math-integer-neg (nth 1 a))
502 (math-scale-int (nth 1 b) 1)))
503 (and (not (consp ediff))
506 (= (math-numdigs (nth 1 a)) calc-internal-prec))))))
508 (defun math-nearly-equal (a b) ; [P N N] [Public]
509 (setq a (math-float a))
510 (setq b (math-float b))
511 (if (eq (car a) 'polar) (setq a (math-complex a)))
512 (if (eq (car b) 'polar) (setq b (math-complex b)))
513 (if (eq (car a) 'cplx)
514 (if (eq (car b) 'cplx)
515 (and (or (math-nearly-equal-float (nth 1 a) (nth 1 b))
516 (and (math-nearly-zerop-float (nth 1 a) (nth 2 a))
517 (math-nearly-zerop-float (nth 1 b) (nth 2 b))))
518 (or (math-nearly-equal-float (nth 2 a) (nth 2 b))
519 (and (math-nearly-zerop-float (nth 2 a) (nth 1 a))
520 (math-nearly-zerop-float (nth 2 b) (nth 1 b)))))
521 (and (math-nearly-equal-float (nth 1 a) b)
522 (math-nearly-zerop-float (nth 2 a) b)))
523 (if (eq (car b) 'cplx)
524 (and (math-nearly-equal-float a (nth 1 b))
525 (math-nearly-zerop-float a (nth 2 b)))
526 (math-nearly-equal-float a b))))
528 ;;; True if A is nearly zero compared to B. [P F F]
529 (defun math-nearly-zerop-float (a b)
531 (<= (+ (math-numdigs (nth 1 a)) (nth 2 a))
532 (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec)))))
534 (defun math-nearly-zerop (a b) ; [P N R] [Public]
535 (setq a (math-float a))
536 (setq b (math-float b))
537 (if (eq (car a) 'cplx)
538 (and (math-nearly-zerop-float (nth 1 a) b)
539 (math-nearly-zerop-float (nth 2 a) b))
540 (if (eq (car a) 'polar)
541 (math-nearly-zerop-float (nth 1 a) b)
542 (math-nearly-zerop-float a b))))
544 ;;; This implementation could be improved, accuracy-wise.
545 (defun math-hypot (a b)
546 (cond ((Math-zerop a) (math-abs b))
547 ((Math-zerop b) (math-abs a))
548 ((not (Math-scalarp a))
549 (if (math-infinitep a)
550 (if (math-infinitep b)
555 (calc-record-why 'scalarp a)
556 (list 'calcFunc-hypot a b)))
557 ((not (Math-scalarp b))
558 (if (math-infinitep b)
560 (calc-record-why 'scalarp b)
561 (list 'calcFunc-hypot a b)))
562 ((and (Math-numberp a) (Math-numberp b))
563 (math-with-extra-prec 1
564 (math-sqrt (math-add (calcFunc-abssqr a) (calcFunc-abssqr b)))))
565 ((eq (car-safe a) 'hms)
566 (if (eq (car-safe b) 'hms) ; this helps sdev's of hms forms
567 (math-to-hms (math-hypot (math-from-hms a 'deg)
568 (math-from-hms b 'deg)))
569 (math-to-hms (math-hypot (math-from-hms a 'deg) b))))
570 ((eq (car-safe b) 'hms)
571 (math-to-hms (math-hypot a (math-from-hms b 'deg))))
573 (defalias 'calcFunc-hypot 'math-hypot)
575 (defun calcFunc-sqr (x)
580 (defun math-nth-root (a n)
581 (cond ((= n 2) (math-sqrt a))
585 (let ((root (math-nth-root-integer a n)))
588 (and (not calc-symbolic-mode)
589 (math-nth-root-float (math-float a) n
590 (math-float (cdr root)))))))
591 ((eq (car-safe a) 'frac)
592 (let* ((num-root (math-nth-root-integer (nth 1 a) n))
593 (den-root (math-nth-root-integer (nth 2 a) n)))
594 (if (and (car num-root) (car den-root))
595 (list 'frac (cdr num-root) (cdr den-root))
596 (and (not calc-symbolic-mode)
599 (math-div-float (math-float (cdr num-root))
600 (math-float (cdr den-root))))))))
601 ((eq (car-safe a) 'float)
602 (and (not calc-symbolic-mode)
603 (math-nth-root-float a n)))
604 ((eq (car-safe a) 'polar)
605 (let ((root (math-nth-root (nth 1 a) n)))
606 (and root (list 'polar root (math-div (nth 2 a) n)))))
609 ;; The variables math-nrf-n, math-nrf-nf and math-nrf-nfm1 are local
610 ;; to math-nth-root-float, but are used by math-nth-root-float-iter,
611 ;; which is called by math-nth-root-float.
614 (defvar math-nrf-nfm1)
616 (defun math-nth-root-float (a math-nrf-n &optional guess)
617 (math-inexact-result)
618 (math-with-extra-prec 1
619 (let ((math-nrf-nf (math-float math-nrf-n))
620 (math-nrf-nfm1 (math-float (1- math-nrf-n))))
621 (math-nth-root-float-iter a (or guess
623 1 (/ (+ (math-numdigs (nth 1 a))
628 (defun math-nth-root-float-iter (a guess)
629 (math-working "root" guess)
630 (let ((g2 (math-div-float (math-add-float (math-mul math-nrf-nfm1 guess)
632 a (math-ipow guess (1- math-nrf-n))))
634 (if (math-nearly-equal-float g2 guess)
636 (math-nth-root-float-iter a g2))))
638 ;; The variable math-nri-n is local to math-nth-root-integer, but
639 ;; is used by math-nth-root-int-iter, which is called by
640 ;; math-nth-root-int.
643 (defun math-nth-root-integer (a math-nri-n &optional guess) ; [I I S]
644 (math-nth-root-int-iter a (or guess
645 (math-scale-int 1 (/ (+ (math-numdigs a)
649 (defun math-nth-root-int-iter (a guess)
650 (math-working "root" guess)
651 (let* ((q (math-idivmod a (math-ipow guess (1- math-nri-n))))
652 (s (math-add (car q) (math-mul (1- math-nri-n) guess)))
653 (g2 (math-idivmod s math-nri-n)))
654 (if (Math-natnum-lessp (car g2) guess)
655 (math-nth-root-int-iter a (car g2))
656 (cons (and (equal (car g2) guess)
661 (defun calcFunc-nroot (x n)
662 (calcFunc-pow x (if (integerp n)
669 ;;;; Transcendental functions.
671 ;;; All of these functions are defined on the complex plane.
672 ;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
674 ;;; Most functions increase calc-internal-prec by 2 digits, then round
675 ;;; down afterward. "-raw" functions use the current precision, require
676 ;;; their arguments to be in float (or complex float) format, and always
677 ;;; work in radians (where applicable).
679 (defun math-to-radians (a) ; [N N]
680 (cond ((eq (car-safe a) 'hms)
681 (math-from-hms a 'rad))
682 ((memq calc-angle-mode '(deg hms))
683 (math-mul a (math-pi-over-180)))
686 (defun math-from-radians (a) ; [N N]
687 (cond ((eq calc-angle-mode 'deg)
689 (math-div a (math-pi-over-180))
690 (list 'calcFunc-deg a)))
691 ((eq calc-angle-mode 'hms)
692 (math-to-hms a 'rad))
695 (defun math-to-radians-2 (a) ; [N N]
696 (cond ((eq (car-safe a) 'hms)
697 (math-from-hms a 'rad))
698 ((memq calc-angle-mode '(deg hms))
699 (if calc-symbolic-mode
700 (math-div (math-mul a '(var pi var-pi)) 180)
701 (math-mul a (math-pi-over-180))))
704 (defun math-from-radians-2 (a) ; [N N]
705 (cond ((memq calc-angle-mode '(deg hms))
706 (if calc-symbolic-mode
707 (math-div (math-mul 180 a) '(var pi var-pi))
708 (math-div a (math-pi-over-180))))
713 ;;; Sine, cosine, and tangent.
715 (defun calcFunc-sin (x) ; [N N] [Public]
716 (cond ((and (integerp x)
717 (if (eq calc-angle-mode 'deg)
720 (aref [0 1 0 -1] (math-mod (/ x 90) 4)))
722 (math-with-extra-prec 2
723 (math-sin-raw (math-to-radians (math-float x)))))
726 (math-with-extra-prec 2
727 (let* ((xx (math-to-radians (math-float (nth 1 x))))
728 (xs (math-to-radians (math-float (nth 2 x))))
729 (sc (math-sin-cos-raw xx)))
730 (math-make-sdev (car sc) (math-mul xs (cdr sc)))))
731 (math-make-sdev (calcFunc-sin (nth 1 x))
732 (math-mul (nth 2 x) (calcFunc-cos (nth 1 x))))))
733 ((and (eq (car x) 'intv) (math-intv-constp x))
734 (calcFunc-cos (math-sub x (math-quarter-circle nil))))
735 ((equal x '(var nan var-nan))
737 (t (calc-record-why 'scalarp x)
738 (list 'calcFunc-sin x))))
740 (defun calcFunc-cos (x) ; [N N] [Public]
741 (cond ((and (integerp x)
742 (if (eq calc-angle-mode 'deg)
745 (aref [1 0 -1 0] (math-mod (/ x 90) 4)))
747 (math-with-extra-prec 2
748 (math-cos-raw (math-to-radians (math-float x)))))
751 (math-with-extra-prec 2
752 (let* ((xx (math-to-radians (math-float (nth 1 x))))
753 (xs (math-to-radians (math-float (nth 2 x))))
754 (sc (math-sin-cos-raw xx)))
755 (math-make-sdev (cdr sc) (math-mul xs (car sc)))))
756 (math-make-sdev (calcFunc-cos (nth 1 x))
757 (math-mul (nth 2 x) (calcFunc-sin (nth 1 x))))))
758 ((and (eq (car x) 'intv) (math-intv-constp x))
759 (math-with-extra-prec 2
760 (let* ((xx (math-to-radians (math-float x)))
761 (na (math-floor (math-div (nth 2 xx) (math-pi))))
762 (nb (math-floor (math-div (nth 3 xx) (math-pi))))
763 (span (math-sub nb na)))
764 (if (memq span '(0 1))
765 (let ((int (math-sort-intv (nth 1 x)
766 (math-cos-raw (nth 2 xx))
767 (math-cos-raw (nth 3 xx)))))
770 (math-make-intv (logior (nth 1 x) 2)
773 (math-make-intv (logior (nth 1 x) 1)
777 (list 'intv 3 -1 1)))))
778 ((equal x '(var nan var-nan))
780 (t (calc-record-why 'scalarp x)
781 (list 'calcFunc-cos x))))
783 (defun calcFunc-sincos (x) ; [V N] [Public]
785 (math-with-extra-prec 2
786 (let ((sc (math-sin-cos-raw (math-to-radians (math-float x)))))
787 (list 'vec (cdr sc) (car sc)))) ; the vector [cos, sin]
788 (list 'vec (calcFunc-sin x) (calcFunc-cos x))))
790 (defun calcFunc-tan (x) ; [N N] [Public]
791 (cond ((and (integerp x)
792 (if (eq calc-angle-mode 'deg)
797 (math-with-extra-prec 2
798 (math-tan-raw (math-to-radians (math-float x)))))
801 (math-with-extra-prec 2
802 (let* ((xx (math-to-radians (math-float (nth 1 x))))
803 (xs (math-to-radians (math-float (nth 2 x))))
804 (sc (math-sin-cos-raw xx)))
805 (if (and (math-zerop (cdr sc)) (not calc-infinite-mode))
807 (calc-record-why "*Division by zero")
808 (list 'calcFunc-tan x))
809 (math-make-sdev (math-div-float (car sc) (cdr sc))
810 (math-div-float xs (math-sqr (cdr sc)))))))
811 (math-make-sdev (calcFunc-tan (nth 1 x))
813 (math-sqr (calcFunc-cos (nth 1 x)))))))
814 ((and (eq (car x) 'intv) (math-intv-constp x))
815 (or (math-with-extra-prec 2
816 (let* ((xx (math-to-radians (math-float x)))
817 (na (math-floor (math-div (math-sub (nth 2 xx)
820 (nb (math-floor (math-div (math-sub (nth 3 xx)
824 (math-sort-intv (nth 1 x)
825 (math-tan-raw (nth 2 xx))
826 (math-tan-raw (nth 3 xx))))))
827 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
828 ((equal x '(var nan var-nan))
830 (t (calc-record-why 'scalarp x)
831 (list 'calcFunc-tan x))))
833 (defun calcFunc-sec (x)
834 (cond ((and (integerp x)
835 (eq calc-angle-mode 'deg)
841 (eq calc-angle-mode 'rad)
845 (math-with-extra-prec 2
846 (math-sec-raw (math-to-radians (math-float x)))))
849 (math-with-extra-prec 2
850 (let* ((xx (math-to-radians (math-float (nth 1 x))))
851 (xs (math-to-radians (math-float (nth 2 x))))
852 (sc (math-sin-cos-raw xx)))
853 (if (and (math-zerop (cdr sc))
854 (not calc-infinite-mode))
856 (calc-record-why "*Division by zero")
857 (list 'calcFunc-sec x))
858 (math-make-sdev (math-div-float '(float 1 0) (cdr sc))
860 (math-mul xs (car sc))
861 (math-sqr (cdr sc)))))))
862 (math-make-sdev (calcFunc-sec (nth 1 x))
865 (calcFunc-sin (nth 1 x)))
866 (math-sqr (calcFunc-cos (nth 1 x)))))))
867 ((and (eq (car x) 'intv)
868 (math-intv-constp x))
869 (math-with-extra-prec 2
870 (let* ((xx (math-to-radians (math-float x)))
871 (na (math-floor (math-div (math-sub (nth 2 xx)
874 (nb (math-floor (math-div (math-sub (nth 3 xx)
877 (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2))))
878 (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2))))
879 (span (math-sub nbb naa)))
880 (if (not (equal na nb))
881 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
882 (let ((int (math-sort-intv (nth 1 x)
883 (math-sec-raw (nth 2 xx))
884 (math-sec-raw (nth 3 xx)))))
886 (if (math-evenp (math-div (math-add naa 1) 2))
887 (math-make-intv (logior (nth 1 int) 2)
890 (math-make-intv (logior (nth 1 int) 1)
894 ((equal x '(var nan var-nan))
896 (t (calc-record-why 'scalarp x)
897 (list 'calcFunc-sec x))))
899 (defun calcFunc-csc (x)
900 (cond ((and (integerp x)
901 (eq calc-angle-mode 'deg)
902 (= (% (- x 90) 180) 0))
903 (if (= (% (- x 90) 360) 0)
907 (math-with-extra-prec 2
908 (math-csc-raw (math-to-radians (math-float x)))))
911 (math-with-extra-prec 2
912 (let* ((xx (math-to-radians (math-float (nth 1 x))))
913 (xs (math-to-radians (math-float (nth 2 x))))
914 (sc (math-sin-cos-raw xx)))
915 (if (and (math-zerop (car sc))
916 (not calc-infinite-mode))
918 (calc-record-why "*Division by zero")
919 (list 'calcFunc-csc x))
920 (math-make-sdev (math-div-float '(float 1 0) (car sc))
922 (math-mul xs (cdr sc))
923 (math-sqr (car sc)))))))
924 (math-make-sdev (calcFunc-csc (nth 1 x))
927 (calcFunc-cos (nth 1 x)))
928 (math-sqr (calcFunc-sin (nth 1 x)))))))
929 ((and (eq (car x) 'intv)
930 (math-intv-constp x))
931 (math-with-extra-prec 2
932 (let* ((xx (math-to-radians (math-float x)))
933 (na (math-floor (math-div (nth 2 xx) (math-pi))))
934 (nb (math-floor (math-div (nth 3 xx) (math-pi))))
935 (naa (math-floor (math-div (nth 2 xx) (math-pi-over-2))))
936 (nbb (math-floor (math-div (nth 3 xx) (math-pi-over-2))))
937 (span (math-sub nbb naa)))
938 (if (not (equal na nb))
939 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
940 (let ((int (math-sort-intv (nth 1 x)
941 (math-csc-raw (nth 2 xx))
942 (math-csc-raw (nth 3 xx)))))
944 (if (math-evenp (math-div naa 2))
945 (math-make-intv (logior (nth 1 int) 2)
948 (math-make-intv (logior (nth 1 int) 1)
952 ((equal x '(var nan var-nan))
954 (t (calc-record-why 'scalarp x)
955 (list 'calcFunc-csc x))))
957 (defun calcFunc-cot (x) ; [N N] [Public]
958 (cond ((and (integerp x)
959 (if (eq calc-angle-mode 'deg)
960 (= (% (- x 90) 180) 0)
964 (math-with-extra-prec 2
965 (math-cot-raw (math-to-radians (math-float x)))))
968 (math-with-extra-prec 2
969 (let* ((xx (math-to-radians (math-float (nth 1 x))))
970 (xs (math-to-radians (math-float (nth 2 x))))
971 (sc (math-sin-cos-raw xx)))
972 (if (and (math-zerop (car sc)) (not calc-infinite-mode))
974 (calc-record-why "*Division by zero")
975 (list 'calcFunc-cot x))
976 (math-make-sdev (math-div-float (cdr sc) (car sc))
977 (math-div-float xs (math-sqr (car sc)))))))
978 (math-make-sdev (calcFunc-cot (nth 1 x))
980 (math-sqr (calcFunc-sin (nth 1 x)))))))
981 ((and (eq (car x) 'intv) (math-intv-constp x))
982 (or (math-with-extra-prec 2
983 (let* ((xx (math-to-radians (math-float x)))
984 (na (math-floor (math-div (nth 2 xx) (math-pi))))
985 (nb (math-floor (math-div (nth 3 xx) (math-pi)))))
987 (math-sort-intv (nth 1 x)
988 (math-cot-raw (nth 2 xx))
989 (math-cot-raw (nth 3 xx))))))
990 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
991 ((equal x '(var nan var-nan))
993 (t (calc-record-why 'scalarp x)
994 (list 'calcFunc-cot x))))
996 (defun math-sin-raw (x) ; [N N]
997 (cond ((eq (car x) 'cplx)
998 (let* ((expx (math-exp-raw (nth 2 x)))
999 (expmx (math-div-float '(float 1 0) expx))
1000 (sc (math-sin-cos-raw (nth 1 x))))
1002 (math-mul-float (car sc)
1003 (math-mul-float (math-add-float expx expmx)
1005 (math-mul-float (cdr sc)
1006 (math-mul-float (math-sub-float expx expmx)
1008 ((eq (car x) 'polar)
1009 (math-polar (math-sin-raw (math-complex x))))
1010 ((Math-integer-negp (nth 1 x))
1011 (math-neg-float (math-sin-raw (math-neg-float x))))
1012 ((math-lessp-float '(float 7 0) x) ; avoid inf loops due to roundoff
1013 (math-sin-raw (math-mod x (math-two-pi))))
1014 (t (math-sin-raw-2 x x))))
1016 (defun math-cos-raw (x) ; [N N]
1017 (if (eq (car-safe x) 'polar)
1018 (math-polar (math-cos-raw (math-complex x)))
1019 (math-sin-raw (math-sub (math-pi-over-2) x))))
1021 (defun math-sec-raw (x) ; [N N]
1022 (cond ((eq (car x) 'cplx)
1023 (let* ((x (math-mul x '(float 1 0)))
1024 (expx (math-exp-raw (nth 2 x)))
1025 (expmx (math-div-float '(float 1 0) expx))
1026 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
1027 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
1028 (sc (math-sin-cos-raw (nth 1 x)))
1030 (math-mul-float (math-sqr (car sc))
1032 (math-mul-float (math-sqr (cdr sc))
1034 (and (not (eq (nth 1 d) 0))
1036 (math-div-float (math-mul-float (cdr sc) ch) d)
1037 (math-div-float (math-mul-float (car sc) sh) d)))))
1038 ((eq (car x) 'polar)
1039 (math-polar (math-sec-raw (math-complex x))))
1041 (let ((cs (math-cos-raw x)))
1044 (math-div-float '(float 1 0) cs))))))
1046 (defun math-csc-raw (x) ; [N N]
1047 (cond ((eq (car x) 'cplx)
1048 (let* ((x (math-mul x '(float 1 0)))
1049 (expx (math-exp-raw (nth 2 x)))
1050 (expmx (math-div-float '(float 1 0) expx))
1051 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
1052 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
1053 (sc (math-sin-cos-raw (nth 1 x)))
1055 (math-mul-float (math-sqr (car sc))
1057 (math-mul-float (math-sqr (cdr sc))
1059 (and (not (eq (nth 1 d) 0))
1061 (math-div-float (math-mul-float (car sc) ch) d)
1062 (math-div-float (math-mul-float (cdr sc) sh) d)))))
1063 ((eq (car x) 'polar)
1064 (math-polar (math-csc-raw (math-complex x))))
1066 (let ((sn (math-sin-raw x)))
1069 (math-div-float '(float 1 0) sn))))))
1071 (defun math-cot-raw (x) ; [N N]
1072 (cond ((eq (car x) 'cplx)
1073 (let* ((x (math-mul x '(float 1 0)))
1074 (expx (math-exp-raw (nth 2 x)))
1075 (expmx (math-div-float '(float 1 0) expx))
1076 (sh (math-mul-float (math-sub-float expx expmx) '(float 5 -1)))
1077 (ch (math-mul-float (math-add-float expx expmx) '(float 5 -1)))
1078 (sc (math-sin-cos-raw (nth 1 x)))
1082 (and (not (eq (nth 1 d) 0))
1085 (math-mul-float (car sc) (cdr sc))
1089 (math-mul-float sh ch)
1091 ((eq (car x) 'polar)
1092 (math-polar (math-cot-raw (math-complex x))))
1094 (let ((sc (math-sin-cos-raw x)))
1095 (if (eq (nth 1 (car sc)) 0)
1096 (math-div (cdr sc) 0)
1097 (math-div-float (cdr sc) (car sc)))))))
1100 ;;; This could use a smarter method: Reduce x as in math-sin-raw, then
1101 ;;; compute either sin(x) or cos(x), whichever is smaller, and compute
1102 ;;; the other using the identity sin(x)^2 + cos(x)^2 = 1.
1103 (defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x))
1104 (cons (math-sin-raw x) (math-cos-raw x)))
1106 (defun math-tan-raw (x) ; [N N]
1107 (cond ((eq (car x) 'cplx)
1108 (let* ((x (math-mul x '(float 2 0)))
1109 (expx (math-exp-raw (nth 2 x)))
1110 (expmx (math-div-float '(float 1 0) expx))
1111 (sc (math-sin-cos-raw (nth 1 x)))
1112 (d (math-add-float (cdr sc)
1113 (math-mul-float (math-add-float expx expmx)
1115 (and (not (eq (nth 1 d) 0))
1117 (math-div-float (car sc) d)
1118 (math-div-float (math-mul-float (math-sub-float expx
1120 '(float 5 -1)) d)))))
1121 ((eq (car x) 'polar)
1122 (math-polar (math-tan-raw (math-complex x))))
1124 (let ((sc (math-sin-cos-raw x)))
1125 (if (eq (nth 1 (cdr sc)) 0)
1126 (math-div (car sc) 0)
1127 (math-div-float (car sc) (cdr sc)))))))
1129 (defun math-sin-raw-2 (x orgx) ; This avoids poss of inf recursion. [F F]
1130 (let ((xmpo2 (math-sub-float (math-pi-over-2) x)))
1131 (cond ((Math-integer-negp (nth 1 xmpo2))
1132 (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi))
1134 ((math-lessp-float (math-pi-over-4) x)
1135 (math-cos-raw-2 xmpo2 orgx))
1136 ((math-lessp-float x (math-neg (math-pi-over-4)))
1137 (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx)))
1138 ((math-nearly-zerop-float x orgx) '(float 0 0))
1139 (calc-symbolic-mode (signal 'inexact-result nil))
1140 (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x)))))))
1142 (defun math-cos-raw-2 (x orgx) ; [F F]
1143 (cond ((math-nearly-zerop-float x orgx) '(float 1 0))
1144 (calc-symbolic-mode (signal 'inexact-result nil))
1145 (t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
1147 (math-add-float '(float 1 0)
1148 (math-mul-float xnegsqr '(float 5 -1)))
1149 24 5 xnegsqr xnegsqr)))))
1151 (defun math-sin-series (sum nfac n x xnegsqr)
1152 (math-working "sin" sum)
1153 (let* ((nextx (math-mul-float x xnegsqr))
1154 (nextsum (math-add-float sum (math-div-float nextx
1155 (math-float nfac)))))
1156 (if (math-nearly-equal-float sum nextsum)
1158 (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
1159 (+ n 2) nextx xnegsqr))))
1162 ;;; Inverse sine, cosine, tangent.
1164 (defun calcFunc-arcsin (x) ; [N N] [Public]
1166 ((and (eq x 1) (eq calc-angle-mode 'deg)) 90)
1167 ((and (eq x -1) (eq calc-angle-mode 'deg)) -90)
1168 (calc-symbolic-mode (signal 'inexact-result nil))
1170 (math-with-extra-prec 2
1171 (math-from-radians (math-arcsin-raw (math-float x)))))
1173 (math-make-sdev (calcFunc-arcsin (nth 1 x))
1177 (math-sub 1 (math-sqr (nth 1 x))))))))
1179 (math-sort-intv (nth 1 x)
1180 (calcFunc-arcsin (nth 2 x))
1181 (calcFunc-arcsin (nth 3 x))))
1182 ((equal x '(var nan var-nan))
1184 (t (calc-record-why 'numberp x)
1185 (list 'calcFunc-arcsin x))))
1187 (defun calcFunc-arccos (x) ; [N N] [Public]
1189 ((and (eq x 0) (eq calc-angle-mode 'deg)) 90)
1190 ((and (eq x -1) (eq calc-angle-mode 'deg)) 180)
1191 (calc-symbolic-mode (signal 'inexact-result nil))
1193 (math-with-extra-prec 2
1194 (math-from-radians (math-arccos-raw (math-float x)))))
1196 (math-make-sdev (calcFunc-arccos (nth 1 x))
1200 (math-sub 1 (math-sqr (nth 1 x))))))))
1202 (math-sort-intv (nth 1 x)
1203 (calcFunc-arccos (nth 2 x))
1204 (calcFunc-arccos (nth 3 x))))
1205 ((equal x '(var nan var-nan))
1207 (t (calc-record-why 'numberp x)
1208 (list 'calcFunc-arccos x))))
1210 (defun calcFunc-arctan (x) ; [N N] [Public]
1212 ((and (eq x 1) (eq calc-angle-mode 'deg)) 45)
1213 ((and (eq x -1) (eq calc-angle-mode 'deg)) -45)
1215 (math-with-extra-prec 2
1216 (math-from-radians (math-arctan-raw (math-float x)))))
1218 (math-make-sdev (calcFunc-arctan (nth 1 x))
1221 (math-add 1 (math-sqr (nth 1 x)))))))
1223 (math-sort-intv (nth 1 x)
1224 (calcFunc-arctan (nth 2 x))
1225 (calcFunc-arctan (nth 3 x))))
1226 ((equal x '(var inf var-inf))
1227 (math-quarter-circle t))
1228 ((equal x '(neg (var inf var-inf)))
1229 (math-neg (math-quarter-circle t)))
1230 ((equal x '(var nan var-nan))
1232 (t (calc-record-why 'numberp x)
1233 (list 'calcFunc-arctan x))))
1235 (defun math-arcsin-raw (x) ; [N N]
1236 (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x)))))
1237 (if (or (memq (car x) '(cplx polar))
1238 (memq (car a) '(cplx polar)))
1239 (math-with-extra-prec 2 ; use extra precision for difficult case
1240 (math-mul '(cplx 0 -1)
1241 (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a))))
1242 (math-arctan2-raw x a))))
1244 (defun math-arccos-raw (x) ; [N N]
1245 (math-sub (math-pi-over-2) (math-arcsin-raw x)))
1247 (defun math-arctan-raw (x) ; [N N]
1248 (cond ((memq (car x) '(cplx polar))
1249 (math-with-extra-prec 2 ; extra-extra
1251 (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x)))
1252 (math-ln-raw (math-add 1 (math-mul '(cplx 0 -1) x))))
1254 ((Math-integer-negp (nth 1 x))
1255 (math-neg-float (math-arctan-raw (math-neg-float x))))
1257 (calc-symbolic-mode (signal 'inexact-result nil))
1258 ((math-equal-int x 1) (math-pi-over-4))
1259 ((math-equal-int x -1) (math-neg (math-pi-over-4)))
1260 ((math-lessp-float '(float 414214 -6) x) ; if x > sqrt(2) - 1, reduce
1261 (if (math-lessp-float '(float 1 0) x)
1262 (math-sub-float (math-mul-float (math-pi) '(float 5 -1))
1263 (math-arctan-raw (math-div-float '(float 1 0) x)))
1264 (math-sub-float (math-mul-float (math-pi) '(float 25 -2))
1265 (math-arctan-raw (math-div-float
1266 (math-sub-float '(float 1 0) x)
1267 (math-add-float '(float 1 0)
1269 (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x))))))
1271 (defun math-arctan-series (sum n x xnegsqr)
1272 (math-working "arctan" sum)
1273 (let* ((nextx (math-mul-float x xnegsqr))
1274 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1275 (if (math-nearly-equal-float sum nextsum)
1277 (math-arctan-series nextsum (+ n 2) nextx xnegsqr))))
1279 (defun calcFunc-arctan2 (y x) ; [F R R] [Public]
1282 (math-with-extra-prec 2
1283 (math-from-radians (math-arctan2-raw (math-float y)
1285 (calc-record-why 'anglep x)
1286 (list 'calcFunc-arctan2 y x))
1287 (if (and (or (math-infinitep x) (math-anglep x))
1288 (or (math-infinitep y) (math-anglep y)))
1303 (calcFunc-arctan2 y x)
1304 '(var nan var-nan)))
1305 (calc-record-why 'anglep y)
1306 (list 'calcFunc-arctan2 y x))))
1308 (defun math-arctan2-raw (y x) ; [F R R]
1309 (cond ((math-zerop y)
1310 (if (math-negp x) (math-pi)
1311 (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0)))
1315 (math-neg (math-pi-over-2))))
1317 (math-arctan-raw (math-div-float y x)))
1319 (math-add-float (math-arctan-raw (math-div-float y x))
1322 (math-sub-float (math-arctan-raw (math-div-float y x))
1325 (defun calcFunc-arcsincos (x) ; [V N] [Public]
1326 (if (and (Math-vectorp x)
1328 (calcFunc-arctan2 (nth 2 x) (nth 1 x))
1329 (math-reject-arg x "*Two-element vector expected")))
1333 ;;; Exponential function.
1335 (defun calcFunc-exp (x) ; [N N] [Public]
1337 ((and (memq x '(1 -1)) calc-symbolic-mode)
1338 (if (eq x 1) '(var e var-e) (math-div 1 '(var e var-e))))
1340 (math-with-extra-prec 2 (math-exp-raw (math-float x))))
1341 ((eq (car-safe x) 'sdev)
1342 (let ((ex (calcFunc-exp (nth 1 x))))
1343 (math-make-sdev ex (math-mul (nth 2 x) ex))))
1344 ((eq (car-safe x) 'intv)
1345 (math-make-intv (nth 1 x) (calcFunc-exp (nth 2 x))
1346 (calcFunc-exp (nth 3 x))))
1347 ((equal x '(var inf var-inf))
1349 ((equal x '(neg (var inf var-inf)))
1351 ((equal x '(var nan var-nan))
1353 (t (calc-record-why 'numberp x)
1354 (list 'calcFunc-exp x))))
1356 (defun calcFunc-expm1 (x) ; [N N] [Public]
1358 ((math-zerop x) '(float 0 0))
1359 (calc-symbolic-mode (signal 'inexact-result nil))
1361 (math-with-extra-prec 2
1362 (let ((x (math-float x)))
1363 (if (and (eq (car x) 'float)
1364 (math-lessp-float x '(float 1 0))
1365 (math-lessp-float '(float -1 0) x))
1366 (math-exp-minus-1-raw x)
1367 (math-add (math-exp-raw x) -1)))))
1368 ((eq (car-safe x) 'sdev)
1370 (let ((ex (calcFunc-expm1 (nth 1 x))))
1371 (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1))))
1372 (math-make-sdev (calcFunc-expm1 (nth 1 x))
1373 (math-mul (nth 2 x) (calcFunc-exp (nth 1 x))))))
1374 ((eq (car-safe x) 'intv)
1375 (math-make-intv (nth 1 x)
1376 (calcFunc-expm1 (nth 2 x))
1377 (calcFunc-expm1 (nth 3 x))))
1378 ((equal x '(var inf var-inf))
1380 ((equal x '(neg (var inf var-inf)))
1382 ((equal x '(var nan var-nan))
1384 (t (calc-record-why 'numberp x)
1385 (list 'calcFunc-expm1 x))))
1387 (defun calcFunc-exp10 (x) ; [N N] [Public]
1390 (math-pow '(float 1 1) x)))
1392 (defun math-exp-raw (x) ; [N N]
1393 (cond ((math-zerop x) '(float 1 0))
1394 (calc-symbolic-mode (signal 'inexact-result nil))
1396 (let ((expx (math-exp-raw (nth 1 x)))
1397 (sc (math-sin-cos-raw (nth 2 x))))
1399 (math-mul-float expx (cdr sc))
1400 (math-mul-float expx (car sc)))))
1401 ((eq (car x) 'polar)
1402 (let ((xc (math-complex x)))
1404 (math-exp-raw (nth 1 xc))
1405 (math-from-radians (nth 2 xc)))))
1406 ((or (math-lessp-float '(float 5 -1) x)
1407 (math-lessp-float x '(float -5 -1)))
1408 (if (math-lessp-float '(float 921035 1) x)
1410 (if (math-lessp-float x '(float -921035 1))
1412 (let* ((two-x (math-mul-float x '(float 2 0)))
1413 (hint (math-scale-int (nth 1 two-x) (nth 2 two-x)))
1414 (hfrac (math-sub-float x (math-mul-float (math-float hint)
1416 (math-mul-float (math-ipow (math-sqrt-e) hint)
1417 (math-add-float '(float 1 0)
1418 (math-exp-minus-1-raw hfrac)))))
1419 (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x)))))
1421 (defun math-exp-minus-1-raw (x) ; [F F]
1422 (math-exp-series x 2 3 x x))
1424 (defun math-exp-series (sum nfac n xpow x)
1425 (math-working "exp" sum)
1426 (let* ((nextx (math-mul-float xpow x))
1427 (nextsum (math-add-float sum (math-div-float nextx
1428 (math-float nfac)))))
1429 (if (math-nearly-equal-float sum nextsum)
1431 (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x))))
1437 (defun calcFunc-ln (x) ; [N N] [Public]
1438 (cond ((math-zerop x)
1439 (if calc-infinite-mode
1440 '(neg (var inf var-inf))
1441 (math-reject-arg x "*Logarithm of zero")))
1444 (math-with-extra-prec 2 (math-ln-raw (math-float x))))
1445 ((eq (car-safe x) 'sdev)
1446 (math-make-sdev (calcFunc-ln (nth 1 x))
1447 (math-div (nth 2 x) (nth 1 x))))
1448 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1449 (Math-zerop (nth 2 x))
1450 (not (math-intv-constp x))))
1451 (let ((calc-infinite-mode t))
1452 (math-make-intv (nth 1 x) (calcFunc-ln (nth 2 x))
1453 (calcFunc-ln (nth 3 x)))))
1454 ((equal x '(var e var-e))
1456 ((and (eq (car-safe x) '^)
1457 (equal (nth 1 x) '(var e var-e))
1458 (math-known-realp (nth 2 x)))
1461 (if (equal x '(var nan var-nan))
1463 '(var inf var-inf)))
1464 (t (calc-record-why 'numberp x)
1465 (list 'calcFunc-ln x))))
1467 (defun calcFunc-log10 (x) ; [N N] [Public]
1468 (cond ((math-equal-int x 1)
1469 (if (math-floatp x) '(float 0 0) 0))
1470 ((and (Math-integerp x)
1472 (let ((res (math-integer-log x 10)))
1474 (setq x (cdr res)))))
1476 ((and (eq (car-safe x) 'frac)
1478 (let ((res (math-integer-log (nth 2 x) 10)))
1480 (setq x (- (cdr res))))))
1483 (if calc-infinite-mode
1484 '(neg (var inf var-inf))
1485 (math-reject-arg x "*Logarithm of zero")))
1486 (calc-symbolic-mode (signal 'inexact-result nil))
1488 (math-with-extra-prec 2
1489 (let ((xf (math-float x)))
1490 (if (eq (nth 1 xf) 0)
1491 (math-reject-arg x "*Logarithm of zero"))
1492 (if (Math-integer-posp (nth 1 xf))
1493 (if (eq (nth 1 xf) 1) ; log10(1*10^n) = n
1494 (math-float (nth 2 xf))
1495 (let ((xdigs (1- (math-numdigs (nth 1 xf)))))
1497 (math-div-float (math-ln-raw-2
1498 (list 'float (nth 1 xf) (- xdigs)))
1500 (math-float (+ (nth 2 xf) xdigs)))))
1501 (math-div (calcFunc-ln xf) (math-ln-10))))))
1502 ((eq (car-safe x) 'sdev)
1503 (math-make-sdev (calcFunc-log10 (nth 1 x))
1505 (math-mul (nth 1 x) (math-ln-10)))))
1506 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1507 (not (math-intv-constp x))))
1508 (math-make-intv (nth 1 x)
1509 (calcFunc-log10 (nth 2 x))
1510 (calcFunc-log10 (nth 3 x))))
1512 (if (equal x '(var nan var-nan))
1514 '(var inf var-inf)))
1515 (t (calc-record-why 'numberp x)
1516 (list 'calcFunc-log10 x))))
1518 (defun calcFunc-log (x &optional b) ; [N N N] [Public]
1519 (cond ((or (null b) (equal b '(var e var-e)))
1520 (math-normalize (list 'calcFunc-ln x)))
1521 ((or (eq b 10) (equal b '(float 1 1)))
1522 (math-normalize (list 'calcFunc-log10 x)))
1524 (if calc-infinite-mode
1525 (math-div (calcFunc-ln x) (calcFunc-ln b))
1526 (math-reject-arg x "*Logarithm of zero")))
1528 (if calc-infinite-mode
1529 (math-div (calcFunc-ln x) (calcFunc-ln b))
1530 (math-reject-arg b "*Logarithm of zero")))
1531 ((math-equal-int b 1)
1532 (if calc-infinite-mode
1533 (math-div (calcFunc-ln x) 0)
1534 (math-reject-arg b "*Logarithm base one")))
1535 ((math-equal-int x 1)
1536 (if (math-floatp b) '(float 0 0) 0))
1537 ((and (Math-ratp x) (Math-ratp b)
1538 (math-posp x) (math-posp b)
1539 (let* ((sign 1) (inv nil)
1540 (xx (if (Math-lessp 1 x)
1544 (bb (if (Math-lessp 1 b)
1546 (setq sign (- sign))
1548 (res (if (Math-lessp xx bb)
1549 (setq inv (math-integer-log bb xx))
1550 (math-integer-log xx bb))))
1553 (math-div 1 (* sign (cdr res)))
1554 (* sign (cdr res)))))))
1556 (calc-symbolic-mode (signal 'inexact-result nil))
1557 ((and (Math-numberp x) (Math-numberp b))
1558 (math-with-extra-prec 2
1559 (math-div (math-ln-raw (math-float x))
1560 (math-log-base-raw b))))
1561 ((and (eq (car-safe x) 'sdev)
1563 (math-make-sdev (calcFunc-log (nth 1 x) b)
1566 (math-log-base-raw b)))))
1567 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1568 (not (math-intv-constp x)))
1570 (math-make-intv (nth 1 x)
1571 (calcFunc-log (nth 2 x) b)
1572 (calcFunc-log (nth 3 x) b)))
1573 ((or (eq (car-safe x) 'intv) (eq (car-safe b) 'intv))
1574 (math-div (calcFunc-ln x) (calcFunc-ln b)))
1575 ((or (math-infinitep x)
1577 (math-div (calcFunc-ln x) (calcFunc-ln b)))
1578 (t (if (Math-numberp b)
1579 (calc-record-why 'numberp x)
1580 (calc-record-why 'numberp b))
1581 (list 'calcFunc-log x b))))
1583 (defun calcFunc-alog (x &optional b)
1584 (cond ((or (null b) (equal b '(var e var-e)))
1585 (math-normalize (list 'calcFunc-exp x)))
1586 (t (math-pow b x))))
1588 (defun calcFunc-ilog (x b)
1589 (if (and (math-natnump x) (not (eq x 0))
1590 (math-natnump b) (not (eq b 0)))
1592 (math-reject-arg x "*Logarithm base one")
1593 (if (Math-natnum-lessp x b)
1595 (cdr (math-integer-log x b))))
1596 (math-floor (calcFunc-log x b))))
1598 (defun math-integer-log (x b)
1599 (let ((pows (list b))
1603 (while (not (Math-lessp x pow))
1604 (setq pows (cons pow pows)
1605 pow (math-sqr pow)))
1606 (setq n (lsh 1 (1- (length pows)))
1609 (while (and (setq pows (cdr pows))
1612 next (math-mul pow (car pows)))
1613 (or (Math-lessp x next)
1616 (cons (equal pow x) sum)))
1619 (defvar math-log-base-cache nil)
1620 (defun math-log-base-raw (b) ; [N N]
1621 (if (not (and (equal (car math-log-base-cache) b)
1622 (eq (nth 1 math-log-base-cache) calc-internal-prec)))
1623 (setq math-log-base-cache (list b calc-internal-prec
1624 (math-ln-raw (math-float b)))))
1625 (nth 2 math-log-base-cache))
1627 (defun calcFunc-lnp1 (x) ; [N N] [Public]
1628 (cond ((Math-equal-int x -1)
1629 (if calc-infinite-mode
1630 '(neg (var inf var-inf))
1631 (math-reject-arg x "*Logarithm of zero")))
1633 ((math-zerop x) '(float 0 0))
1634 (calc-symbolic-mode (signal 'inexact-result nil))
1636 (math-with-extra-prec 2
1637 (let ((x (math-float x)))
1638 (if (and (eq (car x) 'float)
1639 (math-lessp-float x '(float 5 -1))
1640 (math-lessp-float '(float -5 -1) x))
1641 (math-ln-plus-1-raw x)
1642 (math-ln-raw (math-add-float x '(float 1 0)))))))
1643 ((eq (car-safe x) 'sdev)
1644 (math-make-sdev (calcFunc-lnp1 (nth 1 x))
1645 (math-div (nth 2 x) (math-add (nth 1 x) 1))))
1646 ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
1647 (not (math-intv-constp x))))
1648 (math-make-intv (nth 1 x)
1649 (calcFunc-lnp1 (nth 2 x))
1650 (calcFunc-lnp1 (nth 3 x))))
1652 (if (equal x '(var nan var-nan))
1654 '(var inf var-inf)))
1655 (t (calc-record-why 'numberp x)
1656 (list 'calcFunc-lnp1 x))))
1658 (defun math-ln-raw (x) ; [N N] --- must be float format!
1659 (cond ((eq (car-safe x) 'cplx)
1661 (math-mul-float (math-ln-raw
1662 (math-add-float (math-sqr-float (nth 1 x))
1663 (math-sqr-float (nth 2 x))))
1665 (math-arctan2-raw (nth 2 x) (nth 1 x))))
1666 ((eq (car x) 'polar)
1667 (math-polar (list 'cplx
1668 (math-ln-raw (nth 1 x))
1669 (math-to-radians (nth 2 x)))))
1670 ((Math-equal-int x 1)
1672 (calc-symbolic-mode (signal 'inexact-result nil))
1673 ((math-posp (nth 1 x)) ; positive and real
1674 (let ((xdigs (1- (math-numdigs (nth 1 x)))))
1675 (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs)))
1676 (math-mul-float (math-float (+ (nth 2 x) xdigs))
1679 (math-reject-arg x "*Logarithm of zero"))
1680 ((eq calc-complex-mode 'polar) ; negative and real
1682 (list 'cplx ; negative and real
1683 (math-ln-raw (math-neg-float x))
1685 (t (list 'cplx ; negative and real
1686 (math-ln-raw (math-neg-float x))
1689 (defun math-ln-raw-2 (x) ; [F F]
1690 (cond ((math-lessp-float '(float 14 -1) x)
1691 (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1)))
1693 (t ; now .7 < x <= 1.4
1694 (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
1695 (math-add-float x '(float 1 0)))))))
1697 (defun math-ln-raw-3 (x) ; [F F]
1698 (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
1701 ;;; Compute ln((1+x)/(1-x))
1702 (defun math-ln-raw-series (sum n x xsqr)
1703 (math-working "log" sum)
1704 (let* ((nextx (math-mul-float x xsqr))
1705 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1706 (if (math-nearly-equal-float sum nextsum)
1708 (math-ln-raw-series nextsum (+ n 2) nextx xsqr))))
1710 (defun math-ln-plus-1-raw (x)
1711 (math-lnp1-series x 2 x (math-neg x)))
1713 (defun math-lnp1-series (sum n xpow x)
1714 (math-working "lnp1" sum)
1715 (let* ((nextx (math-mul-float xpow x))
1716 (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
1717 (if (math-nearly-equal-float sum nextsum)
1719 (math-lnp1-series nextsum (1+ n) nextx x))))
1721 (defconst math-approx-ln-10
1723 (math-read-number-simple "2.302585092994045684018"))
1724 "An approximation for ln(10).")
1726 (math-defcache math-ln-10 math-approx-ln-10
1727 (math-ln-raw-2 '(float 1 1)))
1729 (defconst math-approx-ln-2
1731 (math-read-number-simple "0.693147180559945309417"))
1732 "An approximation for ln(2).")
1734 (math-defcache math-ln-2 math-approx-ln-2
1735 (math-ln-raw-3 (math-float '(frac 1 3))))
1739 ;;; Hyperbolic functions.
1741 (defun calcFunc-sinh (x) ; [N N] [Public]
1743 (math-expand-formulas
1745 (list '/ (list '- (list 'calcFunc-exp x)
1746 (list 'calcFunc-exp (list 'neg x))) 2)))
1748 (if calc-symbolic-mode (signal 'inexact-result nil))
1749 (math-with-extra-prec 2
1750 (let ((expx (math-exp-raw (math-float x))))
1751 (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1)))))
1752 ((eq (car-safe x) 'sdev)
1753 (math-make-sdev (calcFunc-sinh (nth 1 x))
1754 (math-mul (nth 2 x) (calcFunc-cosh (nth 1 x)))))
1756 (math-sort-intv (nth 1 x)
1757 (calcFunc-sinh (nth 2 x))
1758 (calcFunc-sinh (nth 3 x))))
1759 ((or (equal x '(var inf var-inf))
1760 (equal x '(neg (var inf var-inf)))
1761 (equal x '(var nan var-nan)))
1763 (t (calc-record-why 'numberp x)
1764 (list 'calcFunc-sinh x))))
1765 (put 'calcFunc-sinh 'math-expandable t)
1767 (defun calcFunc-cosh (x) ; [N N] [Public]
1769 (math-expand-formulas
1771 (list '/ (list '+ (list 'calcFunc-exp x)
1772 (list 'calcFunc-exp (list 'neg x))) 2)))
1774 (if calc-symbolic-mode (signal 'inexact-result nil))
1775 (math-with-extra-prec 2
1776 (let ((expx (math-exp-raw (math-float x))))
1777 (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1)))))
1778 ((eq (car-safe x) 'sdev)
1779 (math-make-sdev (calcFunc-cosh (nth 1 x))
1781 (calcFunc-sinh (nth 1 x)))))
1782 ((and (eq (car x) 'intv) (math-intv-constp x))
1783 (setq x (math-abs x))
1784 (math-sort-intv (nth 1 x)
1785 (calcFunc-cosh (nth 2 x))
1786 (calcFunc-cosh (nth 3 x))))
1787 ((or (equal x '(var inf var-inf))
1788 (equal x '(neg (var inf var-inf)))
1789 (equal x '(var nan var-nan)))
1791 (t (calc-record-why 'numberp x)
1792 (list 'calcFunc-cosh x))))
1793 (put 'calcFunc-cosh 'math-expandable t)
1795 (defun calcFunc-tanh (x) ; [N N] [Public]
1797 (math-expand-formulas
1799 (let ((expx (list 'calcFunc-exp x))
1800 (expmx (list 'calcFunc-exp (list 'neg x))))
1802 (list '/ (list '- expx expmx) (list '+ expx expmx))))))
1804 (if calc-symbolic-mode (signal 'inexact-result nil))
1805 (math-with-extra-prec 2
1806 (let* ((expx (calcFunc-exp (math-float x)))
1807 (expmx (math-div 1 expx)))
1808 (math-div (math-sub expx expmx)
1809 (math-add expx expmx)))))
1810 ((eq (car-safe x) 'sdev)
1811 (math-make-sdev (calcFunc-tanh (nth 1 x))
1813 (math-sqr (calcFunc-cosh (nth 1 x))))))
1815 (math-sort-intv (nth 1 x)
1816 (calcFunc-tanh (nth 2 x))
1817 (calcFunc-tanh (nth 3 x))))
1818 ((equal x '(var inf var-inf))
1820 ((equal x '(neg (var inf var-inf)))
1822 ((equal x '(var nan var-nan))
1824 (t (calc-record-why 'numberp x)
1825 (list 'calcFunc-tanh x))))
1826 (put 'calcFunc-tanh 'math-expandable t)
1828 (defun calcFunc-sech (x) ; [N N] [Public]
1830 (math-expand-formulas
1832 (list '/ 2 (list '+ (list 'calcFunc-exp x)
1833 (list 'calcFunc-exp (list 'neg x))))))
1835 (if calc-symbolic-mode (signal 'inexact-result nil))
1836 (math-with-extra-prec 2
1837 (let ((expx (math-exp-raw (math-float x))))
1838 (math-div '(float 2 0) (math-add expx (math-div 1 expx))))))
1839 ((eq (car-safe x) 'sdev)
1840 (math-make-sdev (calcFunc-sech (nth 1 x))
1842 (math-mul (calcFunc-sech (nth 1 x))
1843 (calcFunc-tanh (nth 1 x))))))
1844 ((and (eq (car x) 'intv) (math-intv-constp x))
1845 (setq x (math-abs x))
1846 (math-sort-intv (nth 1 x)
1847 (calcFunc-sech (nth 2 x))
1848 (calcFunc-sech (nth 3 x))))
1849 ((or (equal x '(var inf var-inf))
1850 (equal x '(neg (var inf var-inf))))
1852 ((equal x '(var nan var-nan))
1854 (t (calc-record-why 'numberp x)
1855 (list 'calcFunc-sech x))))
1856 (put 'calcFunc-sech 'math-expandable t)
1858 (defun calcFunc-csch (x) ; [N N] [Public]
1859 (cond ((eq x 0) (math-div 1 0))
1860 (math-expand-formulas
1862 (list '/ 2 (list '- (list 'calcFunc-exp x)
1863 (list 'calcFunc-exp (list 'neg x))))))
1865 (if calc-symbolic-mode (signal 'inexact-result nil))
1866 (math-with-extra-prec 2
1867 (let ((expx (math-exp-raw (math-float x))))
1868 (math-div '(float 2 0) (math-add expx (math-div -1 expx))))))
1869 ((eq (car-safe x) 'sdev)
1870 (math-make-sdev (calcFunc-csch (nth 1 x))
1872 (math-mul (calcFunc-csch (nth 1 x))
1873 (calcFunc-coth (nth 1 x))))))
1875 (if (and (Math-negp (nth 2 x))
1876 (Math-posp (nth 3 x)))
1877 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1878 (math-sort-intv (nth 1 x)
1879 (calcFunc-csch (nth 2 x))
1880 (calcFunc-csch (nth 3 x)))))
1881 ((or (equal x '(var inf var-inf))
1882 (equal x '(neg (var inf var-inf))))
1884 ((equal x '(var nan var-nan))
1886 (t (calc-record-why 'numberp x)
1887 (list 'calcFunc-csch x))))
1888 (put 'calcFunc-csch 'math-expandable t)
1890 (defun calcFunc-coth (x) ; [N N] [Public]
1891 (cond ((eq x 0) (math-div 1 0))
1892 (math-expand-formulas
1894 (let ((expx (list 'calcFunc-exp x))
1895 (expmx (list 'calcFunc-exp (list 'neg x))))
1897 (list '/ (list '+ expx expmx) (list '- expx expmx))))))
1899 (if calc-symbolic-mode (signal 'inexact-result nil))
1900 (math-with-extra-prec 2
1901 (let* ((expx (calcFunc-exp (math-float x)))
1902 (expmx (math-div 1 expx)))
1903 (math-div (math-add expx expmx)
1904 (math-sub expx expmx)))))
1905 ((eq (car-safe x) 'sdev)
1906 (math-make-sdev (calcFunc-coth (nth 1 x))
1908 (math-sqr (calcFunc-sinh (nth 1 x))))))
1910 (if (and (Math-negp (nth 2 x))
1911 (Math-posp (nth 3 x)))
1912 '(intv 3 (neg (var inf var-inf)) (var inf var-inf))
1913 (math-sort-intv (nth 1 x)
1914 (calcFunc-coth (nth 2 x))
1915 (calcFunc-coth (nth 3 x)))))
1916 ((equal x '(var inf var-inf))
1918 ((equal x '(neg (var inf var-inf)))
1920 ((equal x '(var nan var-nan))
1922 (t (calc-record-why 'numberp x)
1923 (list 'calcFunc-coth x))))
1924 (put 'calcFunc-coth 'math-expandable t)
1926 (defun calcFunc-arcsinh (x) ; [N N] [Public]
1928 (math-expand-formulas
1930 (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
1931 (list '+ (list '^ x 2) 1))))))
1933 (if calc-symbolic-mode (signal 'inexact-result nil))
1934 (math-with-extra-prec 2
1935 (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x)
1937 ((eq (car-safe x) 'sdev)
1938 (math-make-sdev (calcFunc-arcsinh (nth 1 x))
1941 (math-add (math-sqr (nth 1 x)) 1)))))
1943 (math-sort-intv (nth 1 x)
1944 (calcFunc-arcsinh (nth 2 x))
1945 (calcFunc-arcsinh (nth 3 x))))
1946 ((or (equal x '(var inf var-inf))
1947 (equal x '(neg (var inf var-inf)))
1948 (equal x '(var nan var-nan)))
1950 (t (calc-record-why 'numberp x)
1951 (list 'calcFunc-arcsinh x))))
1952 (put 'calcFunc-arcsinh 'math-expandable t)
1954 (defun calcFunc-arccosh (x) ; [N N] [Public]
1956 ((and (eq x -1) calc-symbolic-mode)
1958 ((and (eq x 0) calc-symbolic-mode)
1959 (math-div (math-mul '(var pi var-pi) '(var i var-i)) 2))
1960 (math-expand-formulas
1962 (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
1963 (list '- (list '^ x 2) 1))))))
1965 (if calc-symbolic-mode (signal 'inexact-result nil))
1966 (if (Math-equal-int x -1)
1967 (math-imaginary (math-pi))
1968 (math-with-extra-prec 2
1969 (if (or t ; need to do this even in the real case!
1970 (memq (car-safe x) '(cplx polar)))
1971 (let ((xp1 (math-add 1 x))) ; this gets the branch cuts right
1973 (math-add x (math-mul xp1
1980 (math-add x (math-sqrt-raw (math-add (math-sqr x)
1981 '(float -1 0)))))))))
1982 ((eq (car-safe x) 'sdev)
1983 (math-make-sdev (calcFunc-arccosh (nth 1 x))
1986 (math-add (math-sqr (nth 1 x)) -1)))))
1988 (math-sort-intv (nth 1 x)
1989 (calcFunc-arccosh (nth 2 x))
1990 (calcFunc-arccosh (nth 3 x))))
1991 ((or (equal x '(var inf var-inf))
1992 (equal x '(neg (var inf var-inf)))
1993 (equal x '(var nan var-nan)))
1995 (t (calc-record-why 'numberp x)
1996 (list 'calcFunc-arccosh x))))
1997 (put 'calcFunc-arccosh 'math-expandable t)
1999 (defun calcFunc-arctanh (x) ; [N N] [Public]
2001 ((and (Math-equal-int x 1) calc-infinite-mode)
2003 ((and (Math-equal-int x -1) calc-infinite-mode)
2004 '(neg (var inf var-inf)))
2005 (math-expand-formulas
2007 (list 'calcFunc-ln (list '+ 1 x))
2008 (list 'calcFunc-ln (list '- 1 x))) 2))
2010 (if calc-symbolic-mode (signal 'inexact-result nil))
2011 (math-with-extra-prec 2
2012 (if (or (memq (car-safe x) '(cplx polar))
2014 (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x))
2015 (math-ln-raw (math-sub '(float 1 0) x)))
2017 (if (and (math-equal-int x 1) calc-infinite-mode)
2019 (if (and (math-equal-int x -1) calc-infinite-mode)
2020 '(neg (var inf var-inf))
2021 (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x)
2024 ((eq (car-safe x) 'sdev)
2025 (math-make-sdev (calcFunc-arctanh (nth 1 x))
2027 (math-sub 1 (math-sqr (nth 1 x))))))
2029 (math-sort-intv (nth 1 x)
2030 (calcFunc-arctanh (nth 2 x))
2031 (calcFunc-arctanh (nth 3 x))))
2032 ((equal x '(var nan var-nan))
2034 (t (calc-record-why 'numberp x)
2035 (list 'calcFunc-arctanh x))))
2036 (put 'calcFunc-arctanh 'math-expandable t)
2039 ;;; Convert A from HMS or degrees to radians.
2040 (defun calcFunc-rad (a) ; [R R] [Public]
2041 (cond ((or (Math-numberp a)
2043 (math-with-extra-prec 2
2044 (math-mul a (math-pi-over-180))))
2046 (math-from-hms a 'rad))
2048 (math-make-sdev (calcFunc-rad (nth 1 a))
2049 (calcFunc-rad (nth 2 a))))
2050 (math-expand-formulas
2051 (math-div (math-mul a '(var pi var-pi)) 180))
2052 ((math-infinitep a) a)
2053 (t (list 'calcFunc-rad a))))
2054 (put 'calcFunc-rad 'math-expandable t)
2056 ;;; Convert A from HMS or radians to degrees.
2057 (defun calcFunc-deg (a) ; [R R] [Public]
2058 (cond ((or (Math-numberp a)
2060 (math-with-extra-prec 2
2061 (math-div a (math-pi-over-180))))
2063 (math-from-hms a 'deg))
2065 (math-make-sdev (calcFunc-deg (nth 1 a))
2066 (calcFunc-deg (nth 2 a))))
2067 (math-expand-formulas
2068 (math-div (math-mul 180 a) '(var pi var-pi)))
2069 ((math-infinitep a) a)
2070 (t (list 'calcFunc-deg a))))
2071 (put 'calcFunc-deg 'math-expandable t)
2073 (provide 'calc-math)
2075 ;;; arch-tag: c7367e8e-d0b8-4f70-8577-2fb3f31dbb4c
2076 ;;; calc-math.el ends here