diff options
author | Eli Zaretskii <eliz@gnu.org> | 2001-11-06 18:59:06 +0000 |
---|---|---|
committer | Eli Zaretskii <eliz@gnu.org> | 2001-11-06 18:59:06 +0000 |
commit | 136211a997eb94f7dc6f97219052317116e114da (patch) | |
tree | 014fd8ffa0fa5c5d81869ec26426fb262471ee23 /lisp/calc/calc-math.el | |
parent | 0ffbbdeb4464b5b0d63e83fe3f8e91674248d84d (diff) | |
download | emacs-136211a997eb94f7dc6f97219052317116e114da.tar.gz |
Initial import of Calc 2.02f.
Diffstat (limited to 'lisp/calc/calc-math.el')
-rw-r--r-- | lisp/calc/calc-math.el | 1783 |
1 files changed, 1783 insertions, 0 deletions
diff --git a/lisp/calc/calc-math.el b/lisp/calc/calc-math.el new file mode 100644 index 00000000000..c7b841851e1 --- /dev/null +++ b/lisp/calc/calc-math.el @@ -0,0 +1,1783 @@ +;; Calculator for GNU Emacs, part II [calc-math.el] +;; Copyright (C) 1990, 1991, 1992, 1993 Free Software Foundation, Inc. +;; Written by Dave Gillespie, daveg@synaptics.com. + +;; This file is part of GNU Emacs. + +;; GNU Emacs is distributed in the hope that it will be useful, +;; but WITHOUT ANY WARRANTY. No author or distributor +;; accepts responsibility to anyone for the consequences of using it +;; or for whether it serves any particular purpose or works at all, +;; unless he says so in writing. Refer to the GNU Emacs General Public +;; License for full details. + +;; Everyone is granted permission to copy, modify and redistribute +;; GNU Emacs, but only under the conditions described in the +;; GNU Emacs General Public License. A copy of this license is +;; supposed to have been given to you along with GNU Emacs so you +;; can know your rights and responsibilities. It should be in a +;; file named COPYING. Among other things, the copyright notice +;; and this notice must be preserved on all copies. + + + +;; This file is autoloaded from calc-ext.el. +(require 'calc-ext) + +(require 'calc-macs) + +(defun calc-Need-calc-math () nil) + + +(defun calc-sqrt (arg) + (interactive "P") + (calc-slow-wrapper + (if (calc-is-inverse) + (calc-unary-op "^2" 'calcFunc-sqr arg) + (calc-unary-op "sqrt" 'calcFunc-sqrt arg))) +) + +(defun calc-isqrt (arg) + (interactive "P") + (calc-slow-wrapper + (if (calc-is-inverse) + (calc-unary-op "^2" 'calcFunc-sqr arg) + (calc-unary-op "isqt" 'calcFunc-isqrt arg))) +) + + +(defun calc-hypot (arg) + (interactive "P") + (calc-slow-wrapper + (calc-binary-op "hypt" 'calcFunc-hypot arg)) +) + +(defun calc-ln (arg) + (interactive "P") + (calc-invert-func) + (calc-exp arg) +) + +(defun calc-log10 (arg) + (interactive "P") + (calc-hyperbolic-func) + (calc-ln arg) +) + +(defun calc-log (arg) + (interactive "P") + (calc-slow-wrapper + (if (calc-is-inverse) + (calc-binary-op "alog" 'calcFunc-alog arg) + (calc-binary-op "log" 'calcFunc-log arg))) +) + +(defun calc-ilog (arg) + (interactive "P") + (calc-slow-wrapper + (if (calc-is-inverse) + (calc-binary-op "alog" 'calcFunc-alog arg) + (calc-binary-op "ilog" 'calcFunc-ilog arg))) +) + +(defun calc-lnp1 (arg) + (interactive "P") + (calc-invert-func) + (calc-expm1 arg) +) + +(defun calc-exp (arg) + (interactive "P") + (calc-slow-wrapper + (if (calc-is-hyperbolic) + (if (calc-is-inverse) + (calc-unary-op "lg10" 'calcFunc-log10 arg) + (calc-unary-op "10^" 'calcFunc-exp10 arg)) + (if (calc-is-inverse) + (calc-unary-op "ln" 'calcFunc-ln arg) + (calc-unary-op "exp" 'calcFunc-exp arg)))) +) + +(defun calc-expm1 (arg) + (interactive "P") + (calc-slow-wrapper + (if (calc-is-inverse) + (calc-unary-op "ln+1" 'calcFunc-lnp1 arg) + (calc-unary-op "ex-1" 'calcFunc-expm1 arg))) +) + +(defun calc-pi () + (interactive) + (calc-slow-wrapper + (if (calc-is-inverse) + (if (calc-is-hyperbolic) + (if calc-symbolic-mode + (calc-pop-push-record 0 "phi" '(var phi var-phi)) + (calc-pop-push-record 0 "phi" (math-phi))) + (if calc-symbolic-mode + (calc-pop-push-record 0 "gmma" '(var gamma var-gamma)) + (calc-pop-push-record 0 "gmma" (math-gamma-const)))) + (if (calc-is-hyperbolic) + (if calc-symbolic-mode + (calc-pop-push-record 0 "e" '(var e var-e)) + (calc-pop-push-record 0 "e" (math-e))) + (if calc-symbolic-mode + (calc-pop-push-record 0 "pi" '(var pi var-pi)) + (calc-pop-push-record 0 "pi" (math-pi)))))) +) + +(defun calc-sin (arg) + (interactive "P") + (calc-slow-wrapper + (if (calc-is-hyperbolic) + (if (calc-is-inverse) + (calc-unary-op "asnh" 'calcFunc-arcsinh arg) + (calc-unary-op "sinh" 'calcFunc-sinh arg)) + (if (calc-is-inverse) + (calc-unary-op "asin" 'calcFunc-arcsin arg) + (calc-unary-op "sin" 'calcFunc-sin arg)))) +) + +(defun calc-arcsin (arg) + (interactive "P") + (calc-invert-func) + (calc-sin arg) +) + +(defun calc-sinh (arg) + (interactive "P") + (calc-hyperbolic-func) + (calc-sin arg) +) + +(defun calc-arcsinh (arg) + (interactive "P") + (calc-invert-func) + (calc-hyperbolic-func) + (calc-sin arg) +) + +(defun calc-cos (arg) + (interactive "P") + (calc-slow-wrapper + (if (calc-is-hyperbolic) + (if (calc-is-inverse) + (calc-unary-op "acsh" 'calcFunc-arccosh arg) + (calc-unary-op "cosh" 'calcFunc-cosh arg)) + (if (calc-is-inverse) + (calc-unary-op "acos" 'calcFunc-arccos arg) + (calc-unary-op "cos" 'calcFunc-cos arg)))) +) + +(defun calc-arccos (arg) + (interactive "P") + (calc-invert-func) + (calc-cos arg) +) + +(defun calc-cosh (arg) + (interactive "P") + (calc-hyperbolic-func) + (calc-cos arg) +) + +(defun calc-arccosh (arg) + (interactive "P") + (calc-invert-func) + (calc-hyperbolic-func) + (calc-cos arg) +) + +(defun calc-sincos () + (interactive) + (calc-slow-wrapper + (if (calc-is-inverse) + (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1))) + (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1))))) +) + +(defun calc-tan (arg) + (interactive "P") + (calc-slow-wrapper + (if (calc-is-hyperbolic) + (if (calc-is-inverse) + (calc-unary-op "atnh" 'calcFunc-arctanh arg) + (calc-unary-op "tanh" 'calcFunc-tanh arg)) + (if (calc-is-inverse) + (calc-unary-op "atan" 'calcFunc-arctan arg) + (calc-unary-op "tan" 'calcFunc-tan arg)))) +) + +(defun calc-arctan (arg) + (interactive "P") + (calc-invert-func) + (calc-tan arg) +) + +(defun calc-tanh (arg) + (interactive "P") + (calc-hyperbolic-func) + (calc-tan arg) +) + +(defun calc-arctanh (arg) + (interactive "P") + (calc-invert-func) + (calc-hyperbolic-func) + (calc-tan arg) +) + +(defun calc-arctan2 () + (interactive) + (calc-slow-wrapper + (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2)))) +) + +(defun calc-conj (arg) + (interactive "P") + (calc-wrapper + (calc-unary-op "conj" 'calcFunc-conj arg)) +) + +(defun calc-imaginary () + (interactive) + (calc-slow-wrapper + (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1)))) +) + + + +(defun calc-to-degrees (arg) + (interactive "P") + (calc-wrapper + (calc-unary-op ">deg" 'calcFunc-deg arg)) +) + +(defun calc-to-radians (arg) + (interactive "P") + (calc-wrapper + (calc-unary-op ">rad" 'calcFunc-rad arg)) +) + + +(defun calc-degrees-mode (arg) + (interactive "p") + (cond ((= arg 1) + (calc-wrapper + (calc-change-mode 'calc-angle-mode 'deg) + (message "Angles measured in degrees."))) + ((= arg 2) (calc-radians-mode)) + ((= arg 3) (calc-hms-mode)) + (t (error "Prefix argument out of range"))) +) + +(defun calc-radians-mode () + (interactive) + (calc-wrapper + (calc-change-mode 'calc-angle-mode 'rad) + (message "Angles measured in radians.")) +) + + +;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public] +;;; This method takes advantage of the fact that Newton's method starting +;;; with an overestimate always works, even using truncating integer division! +(defun math-isqrt (a) + (cond ((Math-zerop a) a) + ((not (math-natnump a)) + (math-reject-arg a 'natnump)) + ((integerp a) + (math-isqrt-small a)) + (t + (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a))))))) +) + +(defun calcFunc-isqrt (a) + (if (math-realp a) + (math-isqrt (math-floor a)) + (math-floor (math-sqrt a))) +) + + +;;; This returns (flag . result) where the flag is T if A is a perfect square. +(defun math-isqrt-bignum (a) ; [P.l L] + (let ((len (length a))) + (if (= (% len 2) 0) + (let* ((top (nthcdr (- len 2) a))) + (math-isqrt-bignum-iter + a + (math-scale-bignum-3 + (math-bignum-big + (1+ (math-isqrt-small + (+ (* (nth 1 top) 1000) (car top))))) + (1- (/ len 2))))) + (let* ((top (nth (1- len) a))) + (math-isqrt-bignum-iter + a + (math-scale-bignum-3 + (list (1+ (math-isqrt-small top))) + (/ len 2)))))) +) + +(defun math-isqrt-bignum-iter (a guess) ; [l L l] + (math-working "isqrt" (cons 'bigpos guess)) + (let* ((q (math-div-bignum a guess)) + (s (math-add-bignum (car q) guess)) + (g2 (math-div2-bignum s)) + (comp (math-compare-bignum g2 guess))) + (if (< comp 0) + (math-isqrt-bignum-iter a g2) + (cons (and (= comp 0) + (math-zerop-bignum (cdr q)) + (= (% (car s) 2) 0)) + guess))) +) + +(defun math-zerop-bignum (a) + (and (eq (car a) 0) + (progn + (while (eq (car (setq a (cdr a))) 0)) + (null a))) +) + +(defun math-scale-bignum-3 (a n) ; [L L S] + (while (> n 0) + (setq a (cons 0 a) + n (1- n))) + a +) + +(defun math-isqrt-small (a) ; A > 0. [S S] + (let ((g (cond ((>= a 10000) 1000) + ((>= a 100) 100) + (t 10))) + g2) + (while (< (setq g2 (/ (+ g (/ a g)) 2)) g) + (setq g g2)) + g) +) + + + + +;;; Compute the square root of a number. +;;; [T N] if possible, else [F N] if possible, else [C N]. [Public] +(defun math-sqrt (a) + (or + (and (Math-zerop a) a) + (and (math-known-nonposp a) + (math-imaginary (math-sqrt (math-neg a)))) + (and (integerp a) + (let ((sqrt (math-isqrt-small a))) + (if (= (* sqrt sqrt) a) + sqrt + (if calc-symbolic-mode + (list 'calcFunc-sqrt a) + (math-sqrt-float (math-float a) (math-float sqrt)))))) + (and (eq (car-safe a) 'bigpos) + (let* ((res (math-isqrt-bignum (cdr a))) + (sqrt (math-normalize (cons 'bigpos (cdr res))))) + (if (car res) + sqrt + (if calc-symbolic-mode + (list 'calcFunc-sqrt a) + (math-sqrt-float (math-float a) (math-float sqrt)))))) + (and (eq (car-safe a) 'frac) + (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a))))) + (num-sqrt (math-normalize (cons 'bigpos (cdr num-res)))) + (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a))))) + (den-sqrt (math-normalize (cons 'bigpos (cdr den-res))))) + (if (and (car num-res) (car den-res)) + (list 'frac num-sqrt den-sqrt) + (if calc-symbolic-mode + (if (or (car num-res) (car den-res)) + (math-div (if (car num-res) + num-sqrt (list 'calcFunc-sqrt (nth 1 a))) + (if (car den-res) + den-sqrt (list 'calcFunc-sqrt (nth 2 a)))) + (list 'calcFunc-sqrt a)) + (math-sqrt-float (math-float a) + (math-div (math-float num-sqrt) den-sqrt)))))) + (and (eq (car-safe a) 'float) + (if calc-symbolic-mode + (if (= (% (nth 2 a) 2) 0) + (let ((res (math-isqrt-bignum + (cdr (Math-bignum-test (nth 1 a)))))) + (if (car res) + (math-make-float (math-normalize + (cons 'bigpos (cdr res))) + (/ (nth 2 a) 2)) + (signal 'inexact-result nil))) + (signal 'inexact-result nil)) + (math-sqrt-float a))) + (and (eq (car-safe a) 'cplx) + (math-with-extra-prec 2 + (let* ((d (math-abs a)) + (imag (math-sqrt (math-mul (math-sub d (nth 1 a)) + '(float 5 -1))))) + (list 'cplx + (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1))) + (if (math-negp (nth 2 a)) (math-neg imag) imag))))) + (and (eq (car-safe a) 'polar) + (list 'polar + (math-sqrt (nth 1 a)) + (math-mul (nth 2 a) '(float 5 -1)))) + (and (eq (car-safe a) 'sdev) + (let ((sqrt (math-sqrt (nth 1 a)))) + (math-make-sdev sqrt + (math-div (nth 2 a) (math-mul sqrt 2))))) + (and (eq (car-safe a) 'intv) + (not (math-negp (nth 2 a))) + (math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a)))) + (and (eq (car-safe a) '*) + (or (math-known-nonnegp (nth 1 a)) + (math-known-nonnegp (nth 2 a))) + (math-mul (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a)))) + (and (eq (car-safe a) '/) + (or (and (math-known-nonnegp (nth 2 a)) + (math-div (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a)))) + (and (math-known-nonnegp (nth 1 a)) + (not (math-equal-int (nth 1 a) 1)) + (math-mul (math-sqrt (nth 1 a)) + (math-sqrt (math-div 1 (nth 2 a))))))) + (and (eq (car-safe a) '^) + (math-known-evenp (nth 2 a)) + (math-known-realp (nth 1 a)) + (math-abs (math-pow (nth 1 a) (math-div (nth 2 a) 2)))) + (let ((inf (math-infinitep a))) + (and inf + (math-mul (math-sqrt (math-infinite-dir a inf)) inf))) + (progn + (calc-record-why 'numberp a) + (list 'calcFunc-sqrt a))) +) +(fset 'calcFunc-sqrt (symbol-function 'math-sqrt)) + +(defun math-infinite-dir (a &optional inf) + (or inf (setq inf (math-infinitep a))) + (math-normalize (math-expr-subst a inf 1)) +) + +(defun math-sqrt-float (a &optional guess) ; [F F F] + (if calc-symbolic-mode + (signal 'inexact-result nil) + (math-with-extra-prec 1 (math-sqrt-raw a guess))) +) + +(defun math-sqrt-raw (a &optional guess) ; [F F F] + (if (not (Math-posp a)) + (math-sqrt a) + (if (null guess) + (let ((ldiff (- (math-numdigs (nth 1 a)) 6))) + (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff))) + (setq guess (math-make-float (math-isqrt-small + (math-scale-int (nth 1 a) (- ldiff))) + (/ (+ (nth 2 a) ldiff) 2))))) + (math-sqrt-float-iter a guess)) +) + +(defun math-sqrt-float-iter (a guess) ; [F F F] + (math-working "sqrt" guess) + (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess)) + '(float 5 -1)))) + (if (math-nearly-equal-float g2 guess) + g2 + (math-sqrt-float-iter a g2))) +) + +;;; True if A and B differ only in the last digit of precision. [P F F] +(defun math-nearly-equal-float (a b) + (let ((ediff (- (nth 2 a) (nth 2 b)))) + (cond ((= ediff 0) ;; Expanded out for speed + (setq ediff (math-add (Math-integer-neg (nth 1 a)) (nth 1 b))) + (or (eq ediff 0) + (and (not (consp ediff)) + (< ediff 10) + (> ediff -10) + (= (math-numdigs (nth 1 a)) calc-internal-prec)))) + ((= ediff 1) + (setq ediff (math-add (Math-integer-neg (nth 1 b)) + (math-scale-int (nth 1 a) 1))) + (and (not (consp ediff)) + (< ediff 10) + (> ediff -10) + (= (math-numdigs (nth 1 b)) calc-internal-prec))) + ((= ediff -1) + (setq ediff (math-add (Math-integer-neg (nth 1 a)) + (math-scale-int (nth 1 b) 1))) + (and (not (consp ediff)) + (< ediff 10) + (> ediff -10) + (= (math-numdigs (nth 1 a)) calc-internal-prec))))) +) + +(defun math-nearly-equal (a b) ; [P N N] [Public] + (setq a (math-float a)) + (setq b (math-float b)) + (if (eq (car a) 'polar) (setq a (math-complex a))) + (if (eq (car b) 'polar) (setq b (math-complex b))) + (if (eq (car a) 'cplx) + (if (eq (car b) 'cplx) + (and (or (math-nearly-equal-float (nth 1 a) (nth 1 b)) + (and (math-nearly-zerop-float (nth 1 a) (nth 2 a)) + (math-nearly-zerop-float (nth 1 b) (nth 2 b)))) + (or (math-nearly-equal-float (nth 2 a) (nth 2 b)) + (and (math-nearly-zerop-float (nth 2 a) (nth 1 a)) + (math-nearly-zerop-float (nth 2 b) (nth 1 b))))) + (and (math-nearly-equal-float (nth 1 a) b) + (math-nearly-zerop-float (nth 2 a) b))) + (if (eq (car b) 'cplx) + (and (math-nearly-equal-float a (nth 1 b)) + (math-nearly-zerop-float a (nth 2 b))) + (math-nearly-equal-float a b))) +) + +;;; True if A is nearly zero compared to B. [P F F] +(defun math-nearly-zerop-float (a b) + (or (eq (nth 1 a) 0) + (<= (+ (math-numdigs (nth 1 a)) (nth 2 a)) + (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec)))) +) + +(defun math-nearly-zerop (a b) ; [P N R] [Public] + (setq a (math-float a)) + (setq b (math-float b)) + (if (eq (car a) 'cplx) + (and (math-nearly-zerop-float (nth 1 a) b) + (math-nearly-zerop-float (nth 2 a) b)) + (if (eq (car a) 'polar) + (math-nearly-zerop-float (nth 1 a) b) + (math-nearly-zerop-float a b))) +) + +;;; This implementation could be improved, accuracy-wise. +(defun math-hypot (a b) + (cond ((Math-zerop a) (math-abs b)) + ((Math-zerop b) (math-abs a)) + ((not (Math-scalarp a)) + (if (math-infinitep a) + (if (math-infinitep b) + (if (equal a b) + a + '(var nan var-nan)) + a) + (calc-record-why 'scalarp a) + (list 'calcFunc-hypot a b))) + ((not (Math-scalarp b)) + (if (math-infinitep b) + b + (calc-record-why 'scalarp b) + (list 'calcFunc-hypot a b))) + ((and (Math-numberp a) (Math-numberp b)) + (math-with-extra-prec 1 + (math-sqrt (math-add (calcFunc-abssqr a) (calcFunc-abssqr b))))) + ((eq (car-safe a) 'hms) + (if (eq (car-safe b) 'hms) ; this helps sdev's of hms forms + (math-to-hms (math-hypot (math-from-hms a 'deg) + (math-from-hms b 'deg))) + (math-to-hms (math-hypot (math-from-hms a 'deg) b)))) + ((eq (car-safe b) 'hms) + (math-to-hms (math-hypot a (math-from-hms b 'deg)))) + (t nil)) +) +(fset 'calcFunc-hypot (symbol-function 'math-hypot)) + +(defun calcFunc-sqr (x) + (math-pow x 2) +) + + + +(defun math-nth-root (a n) + (cond ((= n 2) (math-sqrt a)) + ((Math-zerop a) a) + ((Math-negp a) nil) + ((Math-integerp a) + (let ((root (math-nth-root-integer a n))) + (if (car root) + (cdr root) + (and (not calc-symbolic-mode) + (math-nth-root-float (math-float a) n + (math-float (cdr root))))))) + ((eq (car-safe a) 'frac) + (let* ((num-root (math-nth-root-integer (nth 1 a) n)) + (den-root (math-nth-root-integer (nth 2 a) n))) + (if (and (car num-root) (car den-root)) + (list 'frac (cdr num-root) (cdr den-root)) + (and (not calc-symbolic-mode) + (math-nth-root-float + (math-float a) n + (math-div-float (math-float (cdr num-root)) + (math-float (cdr den-root)))))))) + ((eq (car-safe a) 'float) + (and (not calc-symbolic-mode) + (math-nth-root-float a n))) + ((eq (car-safe a) 'polar) + (let ((root (math-nth-root (nth 1 a) n))) + (and root (list 'polar root (math-div (nth 2 a) n))))) + (t nil)) +) + +(defun math-nth-root-float (a n &optional guess) + (math-inexact-result) + (math-with-extra-prec 1 + (let ((nf (math-float n)) + (nfm1 (math-float (1- n)))) + (math-nth-root-float-iter a (or guess + (math-make-float + 1 (/ (+ (math-numdigs (nth 1 a)) + (nth 2 a) + (/ n 2)) + n)))))) +) + +(defun math-nth-root-float-iter (a guess) ; uses "n", "nf", "nfm1" + (math-working "root" guess) + (let ((g2 (math-div-float (math-add-float (math-mul nfm1 guess) + (math-div-float + a (math-ipow guess (1- n)))) + nf))) + (if (math-nearly-equal-float g2 guess) + g2 + (math-nth-root-float-iter a g2))) +) + +(defun math-nth-root-integer (a n &optional guess) ; [I I S] + (math-nth-root-int-iter a (or guess + (math-scale-int 1 (/ (+ (math-numdigs a) + (1- n)) + n)))) +) + +(defun math-nth-root-int-iter (a guess) ; uses "n" + (math-working "root" guess) + (let* ((q (math-idivmod a (math-ipow guess (1- n)))) + (s (math-add (car q) (math-mul (1- n) guess))) + (g2 (math-idivmod s n))) + (if (Math-natnum-lessp (car g2) guess) + (math-nth-root-int-iter a (car g2)) + (cons (and (equal (car g2) guess) + (eq (cdr q) 0) + (eq (cdr g2) 0)) + guess))) +) + +(defun calcFunc-nroot (x n) + (calcFunc-pow x (if (integerp n) + (math-make-frac 1 n) + (math-div 1 n))) +) + + + + +;;;; Transcendental functions. + +;;; All of these functions are defined on the complex plane. +;;; (Branch cuts, etc. follow Steele's Common Lisp book.) + +;;; Most functions increase calc-internal-prec by 2 digits, then round +;;; down afterward. "-raw" functions use the current precision, require +;;; their arguments to be in float (or complex float) format, and always +;;; work in radians (where applicable). + +(defun math-to-radians (a) ; [N N] + (cond ((eq (car-safe a) 'hms) + (math-from-hms a 'rad)) + ((memq calc-angle-mode '(deg hms)) + (math-mul a (math-pi-over-180))) + (t a)) +) + +(defun math-from-radians (a) ; [N N] + (cond ((eq calc-angle-mode 'deg) + (if (math-constp a) + (math-div a (math-pi-over-180)) + (list 'calcFunc-deg a))) + ((eq calc-angle-mode 'hms) + (math-to-hms a 'rad)) + (t a)) +) + +(defun math-to-radians-2 (a) ; [N N] + (cond ((eq (car-safe a) 'hms) + (math-from-hms a 'rad)) + ((memq calc-angle-mode '(deg hms)) + (if calc-symbolic-mode + (math-div (math-mul a '(var pi var-pi)) 180) + (math-mul a (math-pi-over-180)))) + (t a)) +) + +(defun math-from-radians-2 (a) ; [N N] + (cond ((memq calc-angle-mode '(deg hms)) + (if calc-symbolic-mode + (math-div (math-mul 180 a) '(var pi var-pi)) + (math-div a (math-pi-over-180)))) + (t a)) +) + + + +;;; Sine, cosine, and tangent. + +(defun calcFunc-sin (x) ; [N N] [Public] + (cond ((and (integerp x) + (if (eq calc-angle-mode 'deg) + (= (% x 90) 0) + (= x 0))) + (aref [0 1 0 -1] (math-mod (/ x 90) 4))) + ((Math-scalarp x) + (math-with-extra-prec 2 + (math-sin-raw (math-to-radians (math-float x))))) + ((eq (car x) 'sdev) + (if (math-constp x) + (math-with-extra-prec 2 + (let* ((xx (math-to-radians (math-float (nth 1 x)))) + (xs (math-to-radians (math-float (nth 2 x)))) + (sc (math-sin-cos-raw xx))) + (math-make-sdev (car sc) (math-mul xs (cdr sc))))) + (math-make-sdev (calcFunc-sin (nth 1 x)) + (math-mul (nth 2 x) (calcFunc-cos (nth 1 x)))))) + ((and (eq (car x) 'intv) (math-intv-constp x)) + (calcFunc-cos (math-sub x (math-quarter-circle nil)))) + ((equal x '(var nan var-nan)) + x) + (t (calc-record-why 'scalarp x) + (list 'calcFunc-sin x))) +) + +(defun calcFunc-cos (x) ; [N N] [Public] + (cond ((and (integerp x) + (if (eq calc-angle-mode 'deg) + (= (% x 90) 0) + (= x 0))) + (aref [1 0 -1 0] (math-mod (/ x 90) 4))) + ((Math-scalarp x) + (math-with-extra-prec 2 + (math-cos-raw (math-to-radians (math-float x))))) + ((eq (car x) 'sdev) + (if (math-constp x) + (math-with-extra-prec 2 + (let* ((xx (math-to-radians (math-float (nth 1 x)))) + (xs (math-to-radians (math-float (nth 2 x)))) + (sc (math-sin-cos-raw xx))) + (math-make-sdev (cdr sc) (math-mul xs (car sc))))) + (math-make-sdev (calcFunc-cos (nth 1 x)) + (math-mul (nth 2 x) (calcFunc-sin (nth 1 x)))))) + ((and (eq (car x) 'intv) (math-intv-constp x)) + (math-with-extra-prec 2 + (let* ((xx (math-to-radians (math-float x))) + (na (math-floor (math-div (nth 2 xx) (math-pi)))) + (nb (math-floor (math-div (nth 3 xx) (math-pi)))) + (span (math-sub nb na))) + (if (memq span '(0 1)) + (let ((int (math-sort-intv (nth 1 x) + (math-cos-raw (nth 2 xx)) + (math-cos-raw (nth 3 xx))))) + (if (eq span 1) + (if (math-evenp na) + (math-make-intv (logior (nth 1 x) 2) + -1 + (nth 3 int)) + (math-make-intv (logior (nth 1 x) 1) + (nth 2 int) + 1)) + int)) + (list 'intv 3 -1 1))))) + ((equal x '(var nan var-nan)) + x) + (t (calc-record-why 'scalarp x) + (list 'calcFunc-cos x))) +) + +(defun calcFunc-sincos (x) ; [V N] [Public] + (if (Math-scalarp x) + (math-with-extra-prec 2 + (let ((sc (math-sin-cos-raw (math-to-radians (math-float x))))) + (list 'vec (cdr sc) (car sc)))) ; the vector [cos, sin] + (list 'vec (calcFunc-sin x) (calcFunc-cos x))) +) + +(defun calcFunc-tan (x) ; [N N] [Public] + (cond ((and (integerp x) + (if (eq calc-angle-mode 'deg) + (= (% x 180) 0) + (= x 0))) + 0) + ((Math-scalarp x) + (math-with-extra-prec 2 + (math-tan-raw (math-to-radians (math-float x))))) + ((eq (car x) 'sdev) + (if (math-constp x) + (math-with-extra-prec 2 + (let* ((xx (math-to-radians (math-float (nth 1 x)))) + (xs (math-to-radians (math-float (nth 2 x)))) + (sc (math-sin-cos-raw xx))) + (if (and (math-zerop (cdr sc)) (not calc-infinite-mode)) + (progn + (calc-record-why "*Division by zero") + (list 'calcFunc-tan x)) + (math-make-sdev (math-div-float (car sc) (cdr sc)) + (math-div-float xs (math-sqr (cdr sc))))))) + (math-make-sdev (calcFunc-tan (nth 1 x)) + (math-div (nth 2 x) + (math-sqr (calcFunc-cos (nth 1 x))))))) + ((and (eq (car x) 'intv) (math-intv-constp x)) + (or (math-with-extra-prec 2 + (let* ((xx (math-to-radians (math-float x))) + (na (math-floor (math-div (math-sub (nth 2 xx) + (math-pi-over-2)) + (math-pi)))) + (nb (math-floor (math-div (math-sub (nth 3 xx) + (math-pi-over-2)) + (math-pi))))) + (and (equal na nb) + (math-sort-intv (nth 1 x) + (math-tan-raw (nth 2 xx)) + (math-tan-raw (nth 3 xx)))))) + '(intv 3 (neg (var inf var-inf)) (var inf var-inf)))) + ((equal x '(var nan var-nan)) + x) + (t (calc-record-why 'scalarp x) + (list 'calcFunc-tan x))) +) + +(defun math-sin-raw (x) ; [N N] + (cond ((eq (car x) 'cplx) + (let* ((expx (math-exp-raw (nth 2 x))) + (expmx (math-div-float '(float 1 0) expx)) + (sc (math-sin-cos-raw (nth 1 x)))) + (list 'cplx + (math-mul-float (car sc) + (math-mul-float (math-add-float expx expmx) + '(float 5 -1))) + (math-mul-float (cdr sc) + (math-mul-float (math-sub-float expx expmx) + '(float 5 -1)))))) + ((eq (car x) 'polar) + (math-polar (math-sin-raw (math-complex x)))) + ((Math-integer-negp (nth 1 x)) + (math-neg-float (math-sin-raw (math-neg-float x)))) + ((math-lessp-float '(float 7 0) x) ; avoid inf loops due to roundoff + (math-sin-raw (math-mod x (math-two-pi)))) + (t (math-sin-raw-2 x x))) +) + +(defun math-cos-raw (x) ; [N N] + (if (eq (car-safe x) 'polar) + (math-polar (math-cos-raw (math-complex x))) + (math-sin-raw (math-sub (math-pi-over-2) x))) +) + +;;; This could use a smarter method: Reduce x as in math-sin-raw, then +;;; compute either sin(x) or cos(x), whichever is smaller, and compute +;;; the other using the identity sin(x)^2 + cos(x)^2 = 1. +(defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x)) + (cons (math-sin-raw x) (math-cos-raw x)) +) + +(defun math-tan-raw (x) ; [N N] + (cond ((eq (car x) 'cplx) + (let* ((x (math-mul x '(float 2 0))) + (expx (math-exp-raw (nth 2 x))) + (expmx (math-div-float '(float 1 0) expx)) + (sc (math-sin-cos-raw (nth 1 x))) + (d (math-add-float (cdr sc) + (math-mul-float (math-add-float expx expmx) + '(float 5 -1))))) + (and (not (eq (nth 1 d) 0)) + (list 'cplx + (math-div-float (car sc) d) + (math-div-float (math-mul-float (math-sub-float expx + expmx) + '(float 5 -1)) d))))) + ((eq (car x) 'polar) + (math-polar (math-tan-raw (math-complex x)))) + (t + (let ((sc (math-sin-cos-raw x))) + (if (eq (nth 1 (cdr sc)) 0) + (math-div (car sc) 0) + (math-div-float (car sc) (cdr sc)))))) +) + +(defun math-sin-raw-2 (x orgx) ; This avoids poss of inf recursion. [F F] + (let ((xmpo2 (math-sub-float (math-pi-over-2) x))) + (cond ((Math-integer-negp (nth 1 xmpo2)) + (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi)) + orgx))) + ((math-lessp-float (math-pi-over-4) x) + (math-cos-raw-2 xmpo2 orgx)) + ((math-lessp-float x (math-neg (math-pi-over-4))) + (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx))) + ((math-nearly-zerop-float x orgx) '(float 0 0)) + (calc-symbolic-mode (signal 'inexact-result nil)) + (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x)))))) +) + +(defun math-cos-raw-2 (x orgx) ; [F F] + (cond ((math-nearly-zerop-float x orgx) '(float 1 0)) + (calc-symbolic-mode (signal 'inexact-result nil)) + (t (let ((xnegsqr (math-neg-float (math-sqr-float x)))) + (math-sin-series + (math-add-float '(float 1 0) + (math-mul-float xnegsqr '(float 5 -1))) + 24 5 xnegsqr xnegsqr)))) +) + +(defun math-sin-series (sum nfac n x xnegsqr) + (math-working "sin" sum) + (let* ((nextx (math-mul-float x xnegsqr)) + (nextsum (math-add-float sum (math-div-float nextx + (math-float nfac))))) + (if (math-nearly-equal-float sum nextsum) + sum + (math-sin-series nextsum (math-mul nfac (* n (1+ n))) + (+ n 2) nextx xnegsqr))) +) + + +;;; Inverse sine, cosine, tangent. + +(defun calcFunc-arcsin (x) ; [N N] [Public] + (cond ((eq x 0) 0) + ((and (eq x 1) (eq calc-angle-mode 'deg)) 90) + ((and (eq x -1) (eq calc-angle-mode 'deg)) -90) + (calc-symbolic-mode (signal 'inexact-result nil)) + ((Math-numberp x) + (math-with-extra-prec 2 + (math-from-radians (math-arcsin-raw (math-float x))))) + ((eq (car x) 'sdev) + (math-make-sdev (calcFunc-arcsin (nth 1 x)) + (math-from-radians + (math-div (nth 2 x) + (math-sqrt + (math-sub 1 (math-sqr (nth 1 x)))))))) + ((eq (car x) 'intv) + (math-sort-intv (nth 1 x) + (calcFunc-arcsin (nth 2 x)) + (calcFunc-arcsin (nth 3 x)))) + ((equal x '(var nan var-nan)) + x) + (t (calc-record-why 'numberp x) + (list 'calcFunc-arcsin x))) +) + +(defun calcFunc-arccos (x) ; [N N] [Public] + (cond ((eq x 1) 0) + ((and (eq x 0) (eq calc-angle-mode 'deg)) 90) + ((and (eq x -1) (eq calc-angle-mode 'deg)) 180) + (calc-symbolic-mode (signal 'inexact-result nil)) + ((Math-numberp x) + (math-with-extra-prec 2 + (math-from-radians (math-arccos-raw (math-float x))))) + ((eq (car x) 'sdev) + (math-make-sdev (calcFunc-arccos (nth 1 x)) + (math-from-radians + (math-div (nth 2 x) + (math-sqrt + (math-sub 1 (math-sqr (nth 1 x)))))))) + ((eq (car x) 'intv) + (math-sort-intv (nth 1 x) + (calcFunc-arccos (nth 2 x)) + (calcFunc-arccos (nth 3 x)))) + ((equal x '(var nan var-nan)) + x) + (t (calc-record-why 'numberp x) + (list 'calcFunc-arccos x))) +) + +(defun calcFunc-arctan (x) ; [N N] [Public] + (cond ((eq x 0) 0) + ((and (eq x 1) (eq calc-angle-mode 'deg)) 45) + ((and (eq x -1) (eq calc-angle-mode 'deg)) -45) + ((Math-numberp x) + (math-with-extra-prec 2 + (math-from-radians (math-arctan-raw (math-float x))))) + ((eq (car x) 'sdev) + (math-make-sdev (calcFunc-arctan (nth 1 x)) + (math-from-radians + (math-div (nth 2 x) + (math-add 1 (math-sqr (nth 1 x))))))) + ((eq (car x) 'intv) + (math-sort-intv (nth 1 x) + (calcFunc-arctan (nth 2 x)) + (calcFunc-arctan (nth 3 x)))) + ((equal x '(var inf var-inf)) + (math-quarter-circle t)) + ((equal x '(neg (var inf var-inf))) + (math-neg (math-quarter-circle t))) + ((equal x '(var nan var-nan)) + x) + (t (calc-record-why 'numberp x) + (list 'calcFunc-arctan x))) +) + +(defun math-arcsin-raw (x) ; [N N] + (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x))))) + (if (or (memq (car x) '(cplx polar)) + (memq (car a) '(cplx polar))) + (math-with-extra-prec 2 ; use extra precision for difficult case + (math-mul '(cplx 0 -1) + (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a)))) + (math-arctan2-raw x a))) +) + +(defun math-arccos-raw (x) ; [N N] + (math-sub (math-pi-over-2) (math-arcsin-raw x)) +) + +(defun math-arctan-raw (x) ; [N N] + (cond ((memq (car x) '(cplx polar)) + (math-with-extra-prec 2 ; extra-extra + (math-div (math-sub + (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x))) + (math-ln-raw (math-add 1 (math-mul '(cplx 0 -1) x)))) + '(cplx 0 2)))) + ((Math-integer-negp (nth 1 x)) + (math-neg-float (math-arctan-raw (math-neg-float x)))) + ((math-zerop x) x) + (calc-symbolic-mode (signal 'inexact-result nil)) + ((math-equal-int x 1) (math-pi-over-4)) + ((math-equal-int x -1) (math-neg (math-pi-over-4))) + ((math-lessp-float '(float 414214 -6) x) ; if x > sqrt(2) - 1, reduce + (if (math-lessp-float '(float 1 0) x) + (math-sub-float (math-mul-float (math-pi) '(float 5 -1)) + (math-arctan-raw (math-div-float '(float 1 0) x))) + (math-sub-float (math-mul-float (math-pi) '(float 25 -2)) + (math-arctan-raw (math-div-float + (math-sub-float '(float 1 0) x) + (math-add-float '(float 1 0) + x)))))) + (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x))))) +) + +(defun math-arctan-series (sum n x xnegsqr) + (math-working "arctan" sum) + (let* ((nextx (math-mul-float x xnegsqr)) + (nextsum (math-add-float sum (math-div-float nextx (math-float n))))) + (if (math-nearly-equal-float sum nextsum) + sum + (math-arctan-series nextsum (+ n 2) nextx xnegsqr))) +) + +(defun calcFunc-arctan2 (y x) ; [F R R] [Public] + (if (Math-anglep y) + (if (Math-anglep x) + (math-with-extra-prec 2 + (math-from-radians (math-arctan2-raw (math-float y) + (math-float x)))) + (calc-record-why 'anglep x) + (list 'calcFunc-arctan2 y x)) + (if (and (or (math-infinitep x) (math-anglep x)) + (or (math-infinitep y) (math-anglep y))) + (progn + (if (math-posp x) + (setq x 1) + (if (math-negp x) + (setq x -1) + (or (math-zerop x) + (setq x nil)))) + (if (math-posp y) + (setq y 1) + (if (math-negp y) + (setq y -1) + (or (math-zerop y) + (setq y nil)))) + (if (and y x) + (calcFunc-arctan2 y x) + '(var nan var-nan))) + (calc-record-why 'anglep y) + (list 'calcFunc-arctan2 y x))) +) + +(defun math-arctan2-raw (y x) ; [F R R] + (cond ((math-zerop y) + (if (math-negp x) (math-pi) + (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0))) + ((math-zerop x) + (if (math-posp y) + (math-pi-over-2) + (math-neg (math-pi-over-2)))) + ((math-posp x) + (math-arctan-raw (math-div-float y x))) + ((math-posp y) + (math-add-float (math-arctan-raw (math-div-float y x)) + (math-pi))) + (t + (math-sub-float (math-arctan-raw (math-div-float y x)) + (math-pi)))) +) + +(defun calcFunc-arcsincos (x) ; [V N] [Public] + (if (and (Math-vectorp x) + (= (length x) 3)) + (calcFunc-arctan2 (nth 2 x) (nth 1 x)) + (math-reject-arg x "*Two-element vector expected")) +) + + + +;;; Exponential function. + +(defun calcFunc-exp (x) ; [N N] [Public] + (cond ((eq x 0) 1) + ((and (memq x '(1 -1)) calc-symbolic-mode) + (if (eq x 1) '(var e var-e) (math-div 1 '(var e var-e)))) + ((Math-numberp x) + (math-with-extra-prec 2 (math-exp-raw (math-float x)))) + ((eq (car-safe x) 'sdev) + (let ((ex (calcFunc-exp (nth 1 x)))) + (math-make-sdev ex (math-mul (nth 2 x) ex)))) + ((eq (car-safe x) 'intv) + (math-make-intv (nth 1 x) (calcFunc-exp (nth 2 x)) + (calcFunc-exp (nth 3 x)))) + ((equal x '(var inf var-inf)) + x) + ((equal x '(neg (var inf var-inf))) + 0) + ((equal x '(var nan var-nan)) + x) + (t (calc-record-why 'numberp x) + (list 'calcFunc-exp x))) +) + +(defun calcFunc-expm1 (x) ; [N N] [Public] + (cond ((eq x 0) 0) + ((math-zerop x) '(float 0 0)) + (calc-symbolic-mode (signal 'inexact-result nil)) + ((Math-numberp x) + (math-with-extra-prec 2 + (let ((x (math-float x))) + (if (and (eq (car x) 'float) + (math-lessp-float x '(float 1 0)) + (math-lessp-float '(float -1 0) x)) + (math-exp-minus-1-raw x) + (math-add (math-exp-raw x) -1))))) + ((eq (car-safe x) 'sdev) + (if (math-constp x) + (let ((ex (calcFunc-expm1 (nth 1 x)))) + (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1)))) + (math-make-sdev (calcFunc-expm1 (nth 1 x)) + (math-mul (nth 2 x) (calcFunc-exp (nth 1 x)))))) + ((eq (car-safe x) 'intv) + (math-make-intv (nth 1 x) + (calcFunc-expm1 (nth 2 x)) + (calcFunc-expm1 (nth 3 x)))) + ((equal x '(var inf var-inf)) + x) + ((equal x '(neg (var inf var-inf))) + -1) + ((equal x '(var nan var-nan)) + x) + (t (calc-record-why 'numberp x) + (list 'calcFunc-expm1 x))) +) + +(defun calcFunc-exp10 (x) ; [N N] [Public] + (if (eq x 0) + 1 + (math-pow '(float 1 1) x)) +) + +(defun math-exp-raw (x) ; [N N] + (cond ((math-zerop x) '(float 1 0)) + (calc-symbolic-mode (signal 'inexact-result nil)) + ((eq (car x) 'cplx) + (let ((expx (math-exp-raw (nth 1 x))) + (sc (math-sin-cos-raw (nth 2 x)))) + (list 'cplx + (math-mul-float expx (cdr sc)) + (math-mul-float expx (car sc))))) + ((eq (car x) 'polar) + (let ((xc (math-complex x))) + (list 'polar + (math-exp-raw (nth 1 xc)) + (math-from-radians (nth 2 xc))))) + ((or (math-lessp-float '(float 5 -1) x) + (math-lessp-float x '(float -5 -1))) + (if (math-lessp-float '(float 921035 1) x) + (math-overflow) + (if (math-lessp-float x '(float -921035 1)) + (math-underflow))) + (let* ((two-x (math-mul-float x '(float 2 0))) + (hint (math-scale-int (nth 1 two-x) (nth 2 two-x))) + (hfrac (math-sub-float x (math-mul-float (math-float hint) + '(float 5 -1))))) + (math-mul-float (math-ipow (math-sqrt-e) hint) + (math-add-float '(float 1 0) + (math-exp-minus-1-raw hfrac))))) + (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x)))) +) + +(defun math-exp-minus-1-raw (x) ; [F F] + (math-exp-series x 2 3 x x) +) + +(defun math-exp-series (sum nfac n xpow x) + (math-working "exp" sum) + (let* ((nextx (math-mul-float xpow x)) + (nextsum (math-add-float sum (math-div-float nextx + (math-float nfac))))) + (if (math-nearly-equal-float sum nextsum) + sum + (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x))) +) + + + +;;; Logarithms. + +(defun calcFunc-ln (x) ; [N N] [Public] + (cond ((math-zerop x) + (if calc-infinite-mode + '(neg (var inf var-inf)) + (math-reject-arg x "*Logarithm of zero"))) + ((eq x 1) 0) + ((Math-numberp x) + (math-with-extra-prec 2 (math-ln-raw (math-float x)))) + ((eq (car-safe x) 'sdev) + (math-make-sdev (calcFunc-ln (nth 1 x)) + (math-div (nth 2 x) (nth 1 x)))) + ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) + (Math-zerop (nth 2 x)) + (not (math-intv-constp x)))) + (let ((calc-infinite-mode t)) + (math-make-intv (nth 1 x) (calcFunc-ln (nth 2 x)) + (calcFunc-ln (nth 3 x))))) + ((equal x '(var e var-e)) + 1) + ((and (eq (car-safe x) '^) + (equal (nth 1 x) '(var e var-e)) + (math-known-realp (nth 2 x))) + (nth 2 x)) + ((math-infinitep x) + (if (equal x '(var nan var-nan)) + x + '(var inf var-inf))) + (t (calc-record-why 'numberp x) + (list 'calcFunc-ln x))) +) + +(defun calcFunc-log10 (x) ; [N N] [Public] + (cond ((math-equal-int x 1) + (if (math-floatp x) '(float 0 0) 0)) + ((and (Math-integerp x) + (math-posp x) + (let ((res (math-integer-log x 10))) + (and (car res) + (setq x (cdr res))))) + x) + ((and (eq (car-safe x) 'frac) + (eq (nth 1 x) 1) + (let ((res (math-integer-log (nth 2 x) 10))) + (and (car res) + (setq x (- (cdr res)))))) + x) + ((math-zerop x) + (if calc-infinite-mode + '(neg (var inf var-inf)) + (math-reject-arg x "*Logarithm of zero"))) + (calc-symbolic-mode (signal 'inexact-result nil)) + ((Math-numberp x) + (math-with-extra-prec 2 + (let ((xf (math-float x))) + (if (eq (nth 1 xf) 0) + (math-reject-arg x "*Logarithm of zero")) + (if (Math-integer-posp (nth 1 xf)) + (if (eq (nth 1 xf) 1) ; log10(1*10^n) = n + (math-float (nth 2 xf)) + (let ((xdigs (1- (math-numdigs (nth 1 xf))))) + (math-add-float + (math-div-float (math-ln-raw-2 + (list 'float (nth 1 xf) (- xdigs))) + (math-ln-10)) + (math-float (+ (nth 2 xf) xdigs))))) + (math-div (calcFunc-ln xf) (math-ln-10)))))) + ((eq (car-safe x) 'sdev) + (math-make-sdev (calcFunc-log10 (nth 1 x)) + (math-div (nth 2 x) + (math-mul (nth 1 x) (math-ln-10))))) + ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) + (not (math-intv-constp x)))) + (math-make-intv (nth 1 x) + (calcFunc-log10 (nth 2 x)) + (calcFunc-log10 (nth 3 x)))) + ((math-infinitep x) + (if (equal x '(var nan var-nan)) + x + '(var inf var-inf))) + (t (calc-record-why 'numberp x) + (list 'calcFunc-log10 x))) +) + +(defun calcFunc-log (x &optional b) ; [N N N] [Public] + (cond ((or (null b) (equal b '(var e var-e))) + (math-normalize (list 'calcFunc-ln x))) + ((or (eq b 10) (equal b '(float 1 1))) + (math-normalize (list 'calcFunc-log10 x))) + ((math-zerop x) + (if calc-infinite-mode + (math-div (calcFunc-ln x) (calcFunc-ln b)) + (math-reject-arg x "*Logarithm of zero"))) + ((math-zerop b) + (if calc-infinite-mode + (math-div (calcFunc-ln x) (calcFunc-ln b)) + (math-reject-arg b "*Logarithm of zero"))) + ((math-equal-int b 1) + (if calc-infinite-mode + (math-div (calcFunc-ln x) 0) + (math-reject-arg b "*Logarithm base one"))) + ((math-equal-int x 1) + (if (or (math-floatp a) (math-floatp b)) '(float 0 0) 0)) + ((and (Math-ratp x) (Math-ratp b) + (math-posp x) (math-posp b) + (let* ((sign 1) (inv nil) + (xx (if (Math-lessp 1 x) + x + (setq sign -1) + (math-div 1 x))) + (bb (if (Math-lessp 1 b) + b + (setq sign (- sign)) + (math-div 1 b))) + (res (if (Math-lessp xx bb) + (setq inv (math-integer-log bb xx)) + (math-integer-log xx bb)))) + (and (car res) + (setq x (if inv + (math-div 1 (* sign (cdr res))) + (* sign (cdr res))))))) + x) + (calc-symbolic-mode (signal 'inexact-result nil)) + ((and (Math-numberp x) (Math-numberp b)) + (math-with-extra-prec 2 + (math-div (math-ln-raw (math-float x)) + (math-log-base-raw b)))) + ((and (eq (car-safe x) 'sdev) + (Math-numberp b)) + (math-make-sdev (calcFunc-log (nth 1 x) b) + (math-div (nth 2 x) + (math-mul (nth 1 x) + (math-log-base-raw b))))) + ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) + (not (math-intv-constp x))) + (math-realp b)) + (math-make-intv (nth 1 x) + (calcFunc-log (nth 2 x) b) + (calcFunc-log (nth 3 x) b))) + ((or (eq (car-safe x) 'intv) (eq (car-safe b) 'intv)) + (math-div (calcFunc-ln x) (calcFunc-ln b))) + ((or (math-infinitep x) + (math-infinitep b)) + (math-div (calcFunc-ln x) (calcFunc-ln b))) + (t (if (Math-numberp b) + (calc-record-why 'numberp x) + (calc-record-why 'numberp b)) + (list 'calcFunc-log x b))) +) + +(defun calcFunc-alog (x &optional b) + (cond ((or (null b) (equal b '(var e var-e))) + (math-normalize (list 'calcFunc-exp x))) + (t (math-pow b x))) +) + +(defun calcFunc-ilog (x b) + (if (and (math-natnump x) (not (eq x 0)) + (math-natnump b) (not (eq b 0))) + (if (eq b 1) + (math-reject-arg x "*Logarithm base one") + (if (Math-natnum-lessp x b) + 0 + (cdr (math-integer-log x b)))) + (math-floor (calcFunc-log x b))) +) + +(defun math-integer-log (x b) + (let ((pows (list b)) + (pow (math-sqr b)) + next + sum n) + (while (not (Math-lessp x pow)) + (setq pows (cons pow pows) + pow (math-sqr pow))) + (setq n (lsh 1 (1- (length pows))) + sum n + pow (car pows)) + (while (and (setq pows (cdr pows)) + (Math-lessp pow x)) + (setq n (/ n 2) + next (math-mul pow (car pows))) + (or (Math-lessp x next) + (setq pow next + sum (+ sum n)))) + (cons (equal pow x) sum)) +) + + +(defun math-log-base-raw (b) ; [N N] + (if (not (and (equal (car math-log-base-cache) b) + (eq (nth 1 math-log-base-cache) calc-internal-prec))) + (setq math-log-base-cache (list b calc-internal-prec + (math-ln-raw (math-float b))))) + (nth 2 math-log-base-cache) +) +(setq math-log-base-cache nil) + +(defun calcFunc-lnp1 (x) ; [N N] [Public] + (cond ((Math-equal-int x -1) + (if calc-infinite-mode + '(neg (var inf var-inf)) + (math-reject-arg x "*Logarithm of zero"))) + ((eq x 0) 0) + ((math-zerop x) '(float 0 0)) + (calc-symbolic-mode (signal 'inexact-result nil)) + ((Math-numberp x) + (math-with-extra-prec 2 + (let ((x (math-float x))) + (if (and (eq (car x) 'float) + (math-lessp-float x '(float 5 -1)) + (math-lessp-float '(float -5 -1) x)) + (math-ln-plus-1-raw x) + (math-ln-raw (math-add-float x '(float 1 0))))))) + ((eq (car-safe x) 'sdev) + (math-make-sdev (calcFunc-lnp1 (nth 1 x)) + (math-div (nth 2 x) (math-add (nth 1 x) 1)))) + ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x)) + (not (math-intv-constp x)))) + (math-make-intv (nth 1 x) + (calcFunc-lnp1 (nth 2 x)) + (calcFunc-lnp1 (nth 3 x)))) + ((math-infinitep x) + (if (equal x '(var nan var-nan)) + x + '(var inf var-inf))) + (t (calc-record-why 'numberp x) + (list 'calcFunc-lnp1 x))) +) + +(defun math-ln-raw (x) ; [N N] --- must be float format! + (cond ((eq (car-safe x) 'cplx) + (list 'cplx + (math-mul-float (math-ln-raw + (math-add-float (math-sqr-float (nth 1 x)) + (math-sqr-float (nth 2 x)))) + '(float 5 -1)) + (math-arctan2-raw (nth 2 x) (nth 1 x)))) + ((eq (car x) 'polar) + (math-polar (list 'cplx + (math-ln-raw (nth 1 x)) + (math-to-radians (nth 2 x))))) + ((Math-equal-int x 1) + '(float 0 0)) + (calc-symbolic-mode (signal 'inexact-result nil)) + ((math-posp (nth 1 x)) ; positive and real + (let ((xdigs (1- (math-numdigs (nth 1 x))))) + (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs))) + (math-mul-float (math-float (+ (nth 2 x) xdigs)) + (math-ln-10))))) + ((math-zerop x) + (math-reject-arg x "*Logarithm of zero")) + ((eq calc-complex-mode 'polar) ; negative and real + (math-polar + (list 'cplx ; negative and real + (math-ln-raw (math-neg-float x)) + (math-pi)))) + (t (list 'cplx ; negative and real + (math-ln-raw (math-neg-float x)) + (math-pi)))) +) + +(defun math-ln-raw-2 (x) ; [F F] + (cond ((math-lessp-float '(float 14 -1) x) + (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1))) + (math-ln-2))) + (t ; now .7 < x <= 1.4 + (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0)) + (math-add-float x '(float 1 0)))))) +) + +(defun math-ln-raw-3 (x) ; [F F] + (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x)) + '(float 2 0)) +) + +;;; Compute ln((1+x)/(1-x)) +(defun math-ln-raw-series (sum n x xsqr) + (math-working "log" sum) + (let* ((nextx (math-mul-float x xsqr)) + (nextsum (math-add-float sum (math-div-float nextx (math-float n))))) + (if (math-nearly-equal-float sum nextsum) + sum + (math-ln-raw-series nextsum (+ n 2) nextx xsqr))) +) + +(defun math-ln-plus-1-raw (x) + (math-lnp1-series x 2 x (math-neg x)) +) + +(defun math-lnp1-series (sum n xpow x) + (math-working "lnp1" sum) + (let* ((nextx (math-mul-float xpow x)) + (nextsum (math-add-float sum (math-div-float nextx (math-float n))))) + (if (math-nearly-equal-float sum nextsum) + sum + (math-lnp1-series nextsum (1+ n) nextx x))) +) + +(math-defcache math-ln-10 (float (bigpos 018 684 045 994 092 585 302 2) -21) + (math-ln-raw-2 '(float 1 1))) + +(math-defcache math-ln-2 (float (bigpos 417 309 945 559 180 147 693) -21) + (math-ln-raw-3 (math-float '(frac 1 3)))) + + + +;;; Hyperbolic functions. + +(defun calcFunc-sinh (x) ; [N N] [Public] + (cond ((eq x 0) 0) + (math-expand-formulas + (math-normalize + (list '/ (list '- (list 'calcFunc-exp x) + (list 'calcFunc-exp (list 'neg x))) 2))) + ((Math-numberp x) + (if calc-symbolic-mode (signal 'inexact-result nil)) + (math-with-extra-prec 2 + (let ((expx (math-exp-raw (math-float x)))) + (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1))))) + ((eq (car-safe x) 'sdev) + (math-make-sdev (calcFunc-sinh (nth 1 x)) + (math-mul (nth 2 x) (calcFunc-cosh (nth 1 x))))) + ((eq (car x) 'intv) + (math-sort-intv (nth 1 x) + (calcFunc-sinh (nth 2 x)) + (calcFunc-sinh (nth 3 x)))) + ((or (equal x '(var inf var-inf)) + (equal x '(neg (var inf var-inf))) + (equal x '(var nan var-nan))) + x) + (t (calc-record-why 'numberp x) + (list 'calcFunc-sinh x))) +) +(put 'calcFunc-sinh 'math-expandable t) + +(defun calcFunc-cosh (x) ; [N N] [Public] + (cond ((eq x 0) 1) + (math-expand-formulas + (math-normalize + (list '/ (list '+ (list 'calcFunc-exp x) + (list 'calcFunc-exp (list 'neg x))) 2))) + ((Math-numberp x) + (if calc-symbolic-mode (signal 'inexact-result nil)) + (math-with-extra-prec 2 + (let ((expx (math-exp-raw (math-float x)))) + (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1))))) + ((eq (car-safe x) 'sdev) + (math-make-sdev (calcFunc-cosh (nth 1 x)) + (math-mul (nth 2 x) + (calcFunc-sinh (nth 1 x))))) + ((and (eq (car x) 'intv) (math-intv-constp x)) + (setq x (math-abs x)) + (math-sort-intv (nth 1 x) + (calcFunc-cosh (nth 2 x)) + (calcFunc-cosh (nth 3 x)))) + ((or (equal x '(var inf var-inf)) + (equal x '(neg (var inf var-inf))) + (equal x '(var nan var-nan))) + (math-abs x)) + (t (calc-record-why 'numberp x) + (list 'calcFunc-cosh x))) +) +(put 'calcFunc-cosh 'math-expandable t) + +(defun calcFunc-tanh (x) ; [N N] [Public] + (cond ((eq x 0) 0) + (math-expand-formulas + (math-normalize + (let ((expx (list 'calcFunc-exp x)) + (expmx (list 'calcFunc-exp (list 'neg x)))) + (math-normalize + (list '/ (list '- expx expmx) (list '+ expx expmx)))))) + ((Math-numberp x) + (if calc-symbolic-mode (signal 'inexact-result nil)) + (math-with-extra-prec 2 + (let* ((expx (calcFunc-exp (math-float x))) + (expmx (math-div 1 expx))) + (math-div (math-sub expx expmx) + (math-add expx expmx))))) + ((eq (car-safe x) 'sdev) + (math-make-sdev (calcFunc-tanh (nth 1 x)) + (math-div (nth 2 x) + (math-sqr (calcFunc-cosh (nth 1 x)))))) + ((eq (car x) 'intv) + (math-sort-intv (nth 1 x) + (calcFunc-tanh (nth 2 x)) + (calcFunc-tanh (nth 3 x)))) + ((equal x '(var inf var-inf)) + 1) + ((equal x '(neg (var inf var-inf))) + -1) + ((equal x '(var nan var-nan)) + x) + (t (calc-record-why 'numberp x) + (list 'calcFunc-tanh x))) +) +(put 'calcFunc-tanh 'math-expandable t) + +(defun calcFunc-arcsinh (x) ; [N N] [Public] + (cond ((eq x 0) 0) + (math-expand-formulas + (math-normalize + (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt + (list '+ (list '^ x 2) 1)))))) + ((Math-numberp x) + (if calc-symbolic-mode (signal 'inexact-result nil)) + (math-with-extra-prec 2 + (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x) + '(float 1 0))))))) + ((eq (car-safe x) 'sdev) + (math-make-sdev (calcFunc-arcsinh (nth 1 x)) + (math-div (nth 2 x) + (math-sqrt + (math-add (math-sqr (nth 1 x)) 1))))) + ((eq (car x) 'intv) + (math-sort-intv (nth 1 x) + (calcFunc-arcsinh (nth 2 x)) + (calcFunc-arcsinh (nth 3 x)))) + ((or (equal x '(var inf var-inf)) + (equal x '(neg (var inf var-inf))) + (equal x '(var nan var-nan))) + x) + (t (calc-record-why 'numberp x) + (list 'calcFunc-arcsinh x))) +) +(put 'calcFunc-arcsinh 'math-expandable t) + +(defun calcFunc-arccosh (x) ; [N N] [Public] + (cond ((eq x 1) 0) + ((and (eq x -1) calc-symbolic-mode) + '(var pi var-pi)) + ((and (eq x 0) calc-symbolic-mode) + (math-div (math-mul '(var pi var-pi) '(var i var-i)) 2)) + (math-expand-formulas + (math-normalize + (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt + (list '- (list '^ x 2) 1)))))) + ((Math-numberp x) + (if calc-symbolic-mode (signal 'inexact-result nil)) + (if (Math-equal-int x -1) + (math-imaginary (math-pi)) + (math-with-extra-prec 2 + (if (or t ; need to do this even in the real case! + (memq (car-safe x) '(cplx polar))) + (let ((xp1 (math-add 1 x))) ; this gets the branch cuts right + (math-ln-raw + (math-add x (math-mul xp1 + (math-sqrt-raw + (math-div (math-sub + x + '(float 1 0)) + xp1)))))) + (math-ln-raw + (math-add x (math-sqrt-raw (math-add (math-sqr x) + '(float -1 0))))))))) + ((eq (car-safe x) 'sdev) + (math-make-sdev (calcFunc-arccosh (nth 1 x)) + (math-div (nth 2 x) + (math-sqrt + (math-add (math-sqr (nth 1 x)) -1))))) + ((eq (car x) 'intv) + (math-sort-intv (nth 1 x) + (calcFunc-arccosh (nth 2 x)) + (calcFunc-arccosh (nth 3 x)))) + ((or (equal x '(var inf var-inf)) + (equal x '(neg (var inf var-inf))) + (equal x '(var nan var-nan))) + x) + (t (calc-record-why 'numberp x) + (list 'calcFunc-arccosh x))) +) +(put 'calcFunc-arccosh 'math-expandable t) + +(defun calcFunc-arctanh (x) ; [N N] [Public] + (cond ((eq x 0) 0) + ((and (Math-equal-int x 1) calc-infinite-mode) + '(var inf var-inf)) + ((and (Math-equal-int x -1) calc-infinite-mode) + '(neg (var inf var-inf))) + (math-expand-formulas + (list '/ (list '- + (list 'calcFunc-ln (list '+ 1 x)) + (list 'calcFunc-ln (list '- 1 x))) 2)) + ((Math-numberp x) + (if calc-symbolic-mode (signal 'inexact-result nil)) + (math-with-extra-prec 2 + (if (or (memq (car-safe x) '(cplx polar)) + (Math-lessp 1 x)) + (math-mul (math-sub (math-ln-raw (math-add '(float 1 0) x)) + (math-ln-raw (math-sub '(float 1 0) x))) + '(float 5 -1)) + (if (and (math-equal-int x 1) calc-infinite-mode) + '(var inf var-inf) + (if (and (math-equal-int x -1) calc-infinite-mode) + '(neg (var inf var-inf)) + (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x) + (math-sub 1 x))) + '(float 5 -1))))))) + ((eq (car-safe x) 'sdev) + (math-make-sdev (calcFunc-arctanh (nth 1 x)) + (math-div (nth 2 x) + (math-sub 1 (math-sqr (nth 1 x)))))) + ((eq (car x) 'intv) + (math-sort-intv (nth 1 x) + (calcFunc-arctanh (nth 2 x)) + (calcFunc-arctanh (nth 3 x)))) + ((equal x '(var nan var-nan)) + x) + (t (calc-record-why 'numberp x) + (list 'calcFunc-arctanh x))) +) +(put 'calcFunc-arctanh 'math-expandable t) + + +;;; Convert A from HMS or degrees to radians. +(defun calcFunc-rad (a) ; [R R] [Public] + (cond ((or (Math-numberp a) + (eq (car a) 'intv)) + (math-with-extra-prec 2 + (math-mul a (math-pi-over-180)))) + ((eq (car a) 'hms) + (math-from-hms a 'rad)) + ((eq (car a) 'sdev) + (math-make-sdev (calcFunc-rad (nth 1 a)) + (calcFunc-rad (nth 2 a)))) + (math-expand-formulas + (math-div (math-mul a '(var pi var-pi)) 180)) + ((math-infinitep a) a) + (t (list 'calcFunc-rad a))) +) +(put 'calcFunc-rad 'math-expandable t) + +;;; Convert A from HMS or radians to degrees. +(defun calcFunc-deg (a) ; [R R] [Public] + (cond ((or (Math-numberp a) + (eq (car a) 'intv)) + (math-with-extra-prec 2 + (math-div a (math-pi-over-180)))) + ((eq (car a) 'hms) + (math-from-hms a 'deg)) + ((eq (car a) 'sdev) + (math-make-sdev (calcFunc-deg (nth 1 a)) + (calcFunc-deg (nth 2 a)))) + (math-expand-formulas + (math-div (math-mul 180 a) '(var pi var-pi))) + ((math-infinitep a) a) + (t (list 'calcFunc-deg a))) +) +(put 'calcFunc-deg 'math-expandable t) + + + + |