daveg@csvax.caltech.edu (David Gillespie) (06/06/90)
Posting-number: Volume 13, Issue 34
Submitted-by: daveg@csvax.caltech.edu (David Gillespie)
Archive-name: gmcalc/part08
---- Cut Here and unpack ----
#!/bin/sh
# this is part 8 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc-ext.el continued
#
CurArch=8
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
exit 1; fi
( read Scheck
if test "$Scheck" != $CurArch
then echo "Please unpack part $Scheck next!"
exit 1;
else exit 0; fi
) < s2_seq_.tmp || exit 1
echo "x - Continuing file calc-ext.el"
sed 's/^X//' << 'SHAR_EOF' >> calc-ext.el
X (- (+ (math-numdigs (nth 1 a))
X (nth 2 a))
X (1+ tol)))))
X ((not (eq (car tol) 'float))
X (if (Math-realp tol)
X (math-to-fraction a (math-float tol))
X (math-reject-arg tol 'realp)))
X ((Math-negp tol)
X (math-to-fraction a (math-neg tol)))
X ((Math-zerop tol)
X (math-to-fraction a 0))
X ((not (math-lessp-float tol '(float 1 0)))
X (math-trunc a))
X ((Math-zerop a)
X 0)
X (t
X (let ((cfrac (math-continued-fraction a tol))
X (calc-prefer-frac t))
X (math-eval-continued-fraction cfrac))))
X)
X(fset 'calcFunc-frac (symbol-function 'math-to-fraction))
X
X(defun math-continued-fraction (a tol)
X (let ((calc-internal-prec (+ calc-internal-prec 2)))
X (let ((cfrac nil)
X (aa a)
X (calc-prefer-frac nil)
X int)
X (while (or (null cfrac)
X (and (not (Math-zerop aa))
X (not (math-lessp-float
X (math-abs
X (math-sub a
X (let ((f (math-eval-continued-fraction
X cfrac)))
X (math-working "Fractionalize" f)
X f)))
X tol))))
X (setq int (math-trunc aa)
X aa (math-sub aa int)
X cfrac (cons int cfrac))
X (or (Math-zerop aa)
X (setq aa (math-div 1 aa))))
X cfrac))
X)
X
X(defun math-eval-continued-fraction (cf)
X (let ((n (car cf))
X (d 1)
X temp)
X (while (setq cf (cdr cf))
X (setq temp (math-add (math-mul (car cf) n) d)
X d n
X n temp))
X (math-div n d))
X)
X
X
X(defun math-clean-number (a &optional prec) ; [X X S] [Public]
X (if prec
X (cond ((Math-messy-integerp prec)
X (math-clean-number a (math-trunc prec)))
X ((or (not (integerp prec))
X (< prec 3))
X (calc-record-why "Precision must be an integer 3 or above")
X (list 'calcFunc-clean a prec))
X ((not (Math-objvecp a))
X (list 'calcFunc-clean a prec))
X (t (let ((calc-internal-prec prec))
X (math-clean-number (math-normalize a)))))
X (cond ((eq (car-safe a) 'polar)
X (let ((theta (math-mod (nth 2 a)
X (if (eq calc-angle-mode 'rad)
X (math-two-pi)
X 360))))
X (math-neg
X (math-neg
X (math-normalize
X (list 'polar (nth 1 a) theta))))))
X ((Math-vectorp a)
X (math-map-vec 'math-clean-number a))
X ((Math-objectp a) a)
X (t (list 'calcFunc-clean a))))
X)
X(fset 'calcFunc-clean (symbol-function 'math-clean-number))
X
X
X
X
X;;;; Logical operations.
X
X(defun calcFunc-eq (a b)
X (let ((res (math-compare a b)))
X (if (= res 0)
X 1
X (if (= res 2)
X (list 'calcFunc-eq a b)
X 0)))
X)
X
X(defun calcFunc-neq (a b)
X (let ((res (math-compare a b)))
X (if (= res 0)
X 0
X (if (= res 2)
X (list 'calcFunc-neq a b)
X 1)))
X)
X
X(defun calcFunc-lt (a b)
X (let ((res (math-compare a b)))
X (if (= res -1)
X 1
X (if (= res 2)
X (list 'calcFunc-lt a b)
X 0)))
X)
X
X(defun calcFunc-gt (a b)
X (let ((res (math-compare a b)))
X (if (= res 1)
X 1
X (if (= res 2)
X (list 'calcFunc-gt a b)
X 0)))
X)
X
X(defun calcFunc-leq (a b)
X (let ((res (math-compare a b)))
X (if (= res 1)
X 0
X (if (= res 2)
X (list 'calcFunc-leq a b)
X 1)))
X)
X
X(defun calcFunc-geq (a b)
X (let ((res (math-compare a b)))
X (if (= res -1)
X 0
X (if (= res 2)
X (list 'calcFunc-geq a b)
X 1)))
X)
X
X(defun calcFunc-land (a b)
X (cond ((Math-zerop a)
X a)
X ((Math-zerop b)
X b)
X ((Math-numberp a)
X b)
X ((Math-numberp b)
X a)
X (t (list 'calcFunc-land a b)))
X)
X
X(defun calcFunc-lor (a b)
X (cond ((Math-zerop a)
X b)
X ((Math-zerop b)
X a)
X ((Math-numberp a)
X a)
X ((Math-numberp b)
X b)
X (t (list 'calcFunc-lor a b)))
X)
X
X(defun calcFunc-lnot (a)
X (if (Math-zerop a)
X 1
X (if (Math-numberp a)
X 0
X (list 'calcFunc-lnot a)))
X)
X
X(defun calcFunc-if (c e1 e2)
X (if (Math-zerop c)
X e2
X (if (Math-numberp c)
X e1
X (list 'calcFunc-if c e1 e2)))
X)
X
X(defun math-normalize-logical-op (a)
X (or (and (eq (car a) 'calcFunc-if)
X (= (length a) 4)
X (let ((a1 (math-normalize (nth 1 a))))
X (if (Math-zerop a1)
X (math-normalize (nth 3 a))
X (if (Math-numberp a1)
X (math-normalize (nth 2 a))
X (list 'calcFunc-if a1 (nth 2 a) (nth 3 a))))))
X a)
X)
X
X(defun calcFunc-in (a b)
X (or (and (eq (car-safe b) 'vec)
X (let ((bb b))
X (while (and (setq bb (cdr bb))
X (not (if (memq (car-safe (car bb)) '(vec intv))
X (eq (calcFunc-in a (car bb)) 1)
X (math-equal a (car bb))))))
X (if bb 1 (and (math-constp a) (math-constp bb) 0))))
X (and (eq (car-safe b) 'intv)
X (let ((res (math-compare a (nth 2 b))))
X (cond ((= res -1)
X 0)
X ((= res 0)
X (if (memq (nth 1 b) '(2 3)) 1 0))
X ((/= res 1)
X nil)
X ((= (setq res (math-compare a (nth 3 b))) 1)
X 0)
X ((= res 0)
X (if (memq (nth 1 b) '(1 3)) 1 0))
X ((/= res -1)
X nil)
X (t 1))))
X (and (math-equal a b)
X 1)
X (and (math-constp a) (math-constp b)
X 0)
X (list 'calcFunc-in a b))
X)
X
X(defun calcFunc-typeof (a)
X (cond ((Math-integerp a) 1)
X ((eq (car a) 'frac) 2)
X ((eq (car a) 'float) 3)
X ((eq (car a) 'hms) 4)
X ((eq (car a) 'cplx) 5)
X ((eq (car a) 'polar) 6)
X ((eq (car a) 'sdev) 7)
X ((eq (car a) 'intv) 8)
X ((eq (car a) 'mod) 9)
X ((eq (car a) 'var) 100)
X ((eq (car a) 'vec) (if (math-matrixp a) 102 101))
X (t (let ((func (assq (car a) '( ( + . calcFunc-add )
X ( - . calcFunc-sub )
X ( * . calcFunc-mul )
X ( / . calcFunc-div )
X ( ^ . calcFunc-pow )
X ( % . calcFunc-mod )
X ( neg . calcFunc-neg )
X ( | . calcFunc-vconcat ) ))))
X (setq func (if func (cdr func) (car a)))
X (math-calcFunc-to-var func))))
X)
X
X(defun calcFunc-integer (a)
X (if (Math-integerp a)
X 1
X (if (Math-objvecp a)
X 0
X (list 'calcFunc-integer a)))
X)
X
X(defun calcFunc-real (a)
X (if (Math-realp a)
X 1
X (if (Math-objvecp a)
X 0
X (list 'calcFunc-real a)))
X)
X
X(defun calcFunc-constant (a)
X (if (math-constp a)
X 1
X (if (Math-objvecp a)
X 0
X (list 'calcFunc-constant a)))
X)
X
X(defun calcFunc-refers (a b)
X (if (math-expr-contains a b)
X 1
X (if (eq (car-safe a) 'var)
X (list 'calcFunc-refers a b)
X 0))
X)
X
X
X
X
X;;;; Complex numbers.
X
X(defun math-to-polar (a) ; [C N] [Public]
X (cond ((Math-vectorp a)
X (math-map-vec 'math-to-polar a))
X ((Math-realp a) a)
X ((Math-numberp a)
X (math-normalize (math-polar a)))
X (t (list 'calcFunc-polar)))
X)
X(fset 'calcFunc-polar (symbol-function 'math-to-polar))
X
X(defun math-to-rectangular (a) ; [N N] [Public]
X (cond ((Math-vectorp a)
X (math-map-vec 'math-to-rectangular a))
X ((Math-realp a) a)
X ((Math-numberp a)
X (math-normalize (math-complex a)))
X (t (list 'calcFunc-rect)))
X)
X(fset 'calcFunc-rect (symbol-function 'math-to-rectangular))
X
X;;; Compute the complex conjugate of A. [O O] [Public]
X(defun math-conj (a)
X (cond ((math-real-objectp a)
X a)
X ((eq (car a) 'cplx)
X (list 'cplx (nth 1 a) (math-neg (nth 2 a))))
X ((eq (car a) 'polar)
X (list 'polar (nth 1 a) (math-neg (nth 2 a))))
X ((eq (car a) 'vec)
X (math-map-vec 'math-conj a))
X ((eq (car a) 'calcFunc-conj)
X (nth 1 a))
X (t (calc-record-why 'numberp a)
X (list 'calcFunc-conj a)))
X)
X(fset 'calcFunc-conj (symbol-function 'math-conj))
X
X;;; Compute the absolute value squared of A. [F N] [Public]
X(defun math-abssqr (a)
X (cond ((Math-realp a)
X (math-sqr a))
X ((eq (car a) 'cplx)
X (math-add (math-sqr (nth 1 a)) (math-sqr (nth 2 a))))
X ((eq (car a) 'polar)
X (math-sqr (nth 1 a)))
X ((and (memq (car a) '(sdev intv)) (math-constp a))
X (math-sqr (math-abs a)))
X ((eq (car a) 'vec)
X (math-reduce-vec 'math-add (math-map-vec 'math-abssqr a)))
X (t (calc-record-why 'numvecp a)
X (list 'calcFunc-abssqr a)))
X)
X(fset 'calcFunc-abssqr (symbol-function 'math-abssqr))
X
X;;; Compute the complex argument of A. [F N] [Public]
X(defun math-cplx-arg (a)
X (cond ((Math-anglep a)
X (if (math-negp a) (math-half-circle nil) 0))
X ((eq (car-safe a) 'cplx)
X (math-arctan2 (nth 2 a) (nth 1 a)))
X ((eq (car-safe a) 'polar)
X (nth 2 a))
X ((eq (car a) 'vec)
X (math-map-vec 'math-cplx-arg a))
X (t (calc-record-why 'numvecp a)
X (list 'calcFunc-arg a)))
X)
X(fset 'calcFunc-arg (symbol-function 'math-cplx-arg))
X
X;;; Extract the real or complex part of a complex number. [R N] [Public]
X;;; Also extracts the real part of a modulo form.
X(defun math-real-part (a)
X (cond ((memq (car-safe a) '(mod sdev))
X (nth 1 a))
X ((math-real-objectp a) a)
X ((eq (car a) 'cplx)
X (nth 1 a))
X ((eq (car a) 'polar)
X (math-mul (nth 1 a) (math-cos (nth 2 a))))
X ((eq (car a) 'vec)
X (math-map-vec 'math-real-part a))
X (t (calc-record-why 'numberp a)
X (list 'calcFunc-re a)))
X)
X(fset 'calcFunc-re (symbol-function 'math-real-part))
X
X(defun math-imag-part (a)
X (cond ((math-real-objectp a)
X (if (math-floatp a) '(float 0 0) 0))
X ((eq (car a) 'cplx)
X (nth 2 a))
X ((eq (car a) 'polar)
X (math-mul (nth 1 a) (math-sin (nth 2 a))))
X ((eq (car a) 'vec)
X (math-map-vec 'math-imag-part a))
X (t (calc-record-why 'numberp a)
X (list 'calcFunc-im a)))
X)
X(fset 'calcFunc-im (symbol-function 'math-imag-part))
X
X
X
X;;;; Transcendental functions.
X
X;;; All of these functions are defined on the complex plane.
X;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
X
X;;; Most functions increase calc-internal-prec by 2 digits, then round
X;;; down afterward. "-raw" functions use the current precision, require
X;;; their arguments to be in float (or complex float) format, and always
X;;; work in radians (where applicable).
X
X(defun math-to-radians (a) ; [N N]
X (cond ((eq (car-safe a) 'hms)
X (math-from-hms a 'rad))
X ((memq calc-angle-mode '(deg hms))
X (math-mul a (math-pi-over-180)))
X (t a))
X)
X
X(defun math-from-radians (a) ; [N N]
X (cond ((eq calc-angle-mode 'deg)
X (if (math-constp a)
X (math-div a (math-pi-over-180))
X (list 'calcFunc-deg a)))
X ((eq calc-angle-mode 'hms)
X (math-to-hms a 'rad))
X (t a))
X)
X
X(defun math-to-radians-2 (a) ; [N N]
X (cond ((eq (car-safe a) 'hms)
X (math-from-hms a 'rad))
X ((memq calc-angle-mode '(deg hms))
X (if calc-symbolic-mode
X (math-div (math-mul a '(var pi var-pi)) 180)
X (math-mul a (math-pi-over-180))))
X (t a))
X)
X
X(defun math-from-radians-2 (a) ; [N N]
X (cond ((memq calc-angle-mode '(deg hms))
X (if calc-symbolic-mode
X (math-div (math-mul 180 a) '(var pi var-pi))
X (math-div a (math-pi-over-180))))
X (t a))
X)
X
X
X
X;;; Sine, cosine, and tangent.
X
X(defun math-sin (x) ; [N N] [Public]
X (cond ((Math-scalarp x)
X (math-with-extra-prec 2
X (math-sin-raw (math-to-radians (math-float x)))))
X ((eq (car x) 'sdev)
X (if (math-constp x)
X (math-with-extra-prec 2
X (let* ((xx (math-to-radians (math-float (nth 1 x))))
X (xs (math-to-radians (math-float (nth 2 x))))
X (sc (math-sin-cos-raw xx)))
X (math-make-sdev (car sc)
X (math-mul xs (math-abs (cdr sc))))))
X (math-make-sdev (math-sin (nth 1 x))
X (math-mul (nth 2 x)
X (math-abs (math-cos (nth 1 x)))))))
X ((and (eq (car x) 'intv) (math-constp x))
X (math-cos (math-sub x (math-quarter-circle nil))))
X (t (calc-record-why 'scalarp x)
X (list 'calcFunc-sin x)))
X)
X(fset 'calcFunc-sin (symbol-function 'math-sin))
X
X(defun math-cos (x) ; [N N] [Public]
X (cond ((Math-scalarp x)
X (math-with-extra-prec 2
X (math-cos-raw (math-to-radians (math-float x)))))
X ((eq (car x) 'sdev)
X (if (math-constp x)
X (math-with-extra-prec 2
X (let* ((xx (math-to-radians (math-float (nth 1 x))))
X (xs (math-to-radians (math-float (nth 2 x))))
X (sc (math-sin-cos-raw xx)))
X (math-make-sdev (cdr sc)
X (math-mul xs (math-abs (car sc))))))
X (math-make-sdev (math-cos (nth 1 x))
X (math-mul (nth 2 x)
X (math-abs (math-sin (nth 1 x)))))))
X ((and (eq (car x) 'intv) (math-constp x))
X (math-with-extra-prec 2
X (let* ((xx (math-to-radians (math-float x)))
X (na (math-floor (math-div (nth 2 xx) (math-pi))))
X (nb (math-floor (math-div (nth 3 xx) (math-pi))))
X (span (math-sub nb na)))
X (if (memq span '(0 1))
X (let ((int (math-sort-intv (nth 1 x)
X (math-cos-raw (nth 2 xx))
X (math-cos-raw (nth 3 xx)))))
X (if (eq span 1)
X (if (math-evenp na)
X (math-make-intv (logior (nth 1 x) 2)
X -1
X (nth 3 int))
X (math-make-intv (logior (nth 1 x) 1)
X (nth 2 int)
X 1))
X int))
X (list 'intv 3 -1 1)))))
X (t (calc-record-why 'scalarp x)
X (list 'calcFunc-cos x)))
X)
X(fset 'calcFunc-cos (symbol-function 'math-cos))
X
X(defun math-sin-cos (x) ; [V N] [Public]
X (if (Math-scalarp x)
X (math-with-extra-prec 2
X (let ((sc (math-sin-cos-raw (math-to-radians (math-float x)))))
X (list 'vec (cdr sc) (car sc)))) ; the vector [cos, sin]
X (list 'vec (math-sin x) (math-cos x)))
X)
X(fset 'calcFunc-sincos (symbol-function 'math-sin-cos))
X
X(defun math-tan (x) ; [N N] [Public]
X (cond ((Math-scalarp x)
X (math-with-extra-prec 2
X (math-tan-raw (math-to-radians (math-float x)))))
X ((eq (car x) 'sdev)
X (if (math-constp x)
X (math-with-extra-prec 2
X (let* ((xx (math-to-radians (math-float (nth 1 x))))
X (xs (math-to-radians (math-float (nth 2 x))))
X (sc (math-sin-cos-raw xx)))
X (if (math-zerop (cdr sc))
X (progn
X (calc-record-why "Division by zero")
X (list 'calcFunc-tan x))
X (math-make-sdev (math-div-float (car sc) (cdr sc))
X (math-div-float xs (math-sqr (cdr sc)))))))
X (math-make-sdev (math-tan (nth 1 x))
X (math-div (nth 2 x)
X (math-sqr (math-cos (nth 1 x)))))))
X ((and (eq (car x) 'intv) (math-constp x))
X (or (math-with-extra-prec 2
X (let* ((xx (math-to-radians (math-float x)))
X (na (math-floor (math-div (math-sub (nth 2 xx)
X (math-pi-over-2))
X (math-pi))))
X (nb (math-floor (math-div (math-sub (nth 3 xx)
X (math-pi-over-2))
X (math-pi)))))
X (and (equal na nb)
X (math-sort-intv (nth 1 x)
X (math-tan-raw (nth 2 xx))
X (math-tan-raw (nth 3 xx))))))
X (progn
X (calc-record-why "Infinite interval" x)
X (list 'calcFunc-tan x))))
X (t (calc-record-why 'scalarp x)
X (list 'calcFunc-tan x)))
X)
X(fset 'calcFunc-tan (symbol-function 'math-tan))
X
X(defun math-sin-raw (x) ; [N N]
X (cond ((eq (car-safe x) 'cplx)
X (let* ((expx (math-exp-raw (nth 2 x)))
X (expmx (math-div-float '(float 1 0) expx))
X (sc (math-sin-cos-raw (nth 1 x))))
X (list 'cplx
X (math-mul-float (car sc)
X (math-mul-float (math-sub expx expmx)
X '(float 5 -1)))
X (math-mul-float (cdr sc)
X (math-mul-float (math-add-float expx expmx)
X '(float 5 -1))))))
X ((eq (car-safe x) 'polar)
X (math-polar (math-sin-raw (math-complex x))))
X ((Math-integer-negp (nth 1 x))
X (math-neg-float (math-sin-raw (math-neg-float x))))
X ((math-lessp-float '(float 7 0) x) ; avoid inf loops due to roundoff
X (math-sin-raw (math-mod x (math-two-pi))))
X (t (math-sin-raw-2 x x)))
X)
X
X(defun math-cos-raw (x) ; [N N]
X (if (eq (car-safe x) 'polar)
X (math-polar (math-cos-raw (math-complex x)))
X (math-sin-raw (math-sub-float (math-pi-over-2) x)))
X)
X
X;;; This could use a smarter method: Reduce x as in math-sin-raw, then
X;;; compute either sin(x) or cos(x), whichever is smaller, and compute
X;;; the other using the identity sin(x)^2 + cos(x)^2 = 1.
X(defun math-sin-cos-raw (x) ; [F.F F] (result is (sin x . cos x))
X (cons (math-sin-raw x) (math-cos-raw x))
X)
X
X(defun math-tan-raw (x) ; [N N]
X (cond ((eq (car-safe x) 'cplx)
X (let* ((x (math-mul-float x '(float 2 0)))
X (expx (math-exp-raw (nth 2 x)))
X (expmx (math-div-float '(float 1 0) expx))
X (sc (math-sin-cos-raw (nth 1 x)))
X (d (math-add-float (cdr sc)
X (math-mul-float (math-add-float expx expmx)
X '(float 5 -1)))))
X (and (not (eq (nth 1 d) 0))
X (list 'cplx
X (math-div-float (car sc) d)
X (math-div-float (math-mul-float (math-add expx expmx)
X '(float 5 -1)) d)))))
X ((eq (car-safe x) 'polar)
X (math-polar (math-tan-raw (math-complex x))))
X (t
X (let ((sc (math-sin-cos-raw x)))
X (if (eq (nth 1 (cdr sc)) 0)
X (math-reject-arg x "Division by zero")
X (math-div-float (car sc) (cdr sc))))))
X)
X
X(defun math-sin-raw-2 (x orgx) ; This avoids poss of inf recursion. [F F]
X (let ((xmpo2 (math-sub-float (math-pi-over-2) x)))
X (cond ((Math-integer-negp (nth 1 xmpo2))
X (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi))
X orgx)))
X ((math-lessp-float (math-pi-over-4) x)
X (math-cos-raw-2 xmpo2 orgx))
X ((math-lessp-float x (math-neg (math-pi-over-4)))
X (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx)))
X ((math-nearly-zerop-float x orgx) '(float 0 0))
X (calc-symbolic-mode (signal 'inexact-result nil))
X (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x))))))
X)
X
X(defun math-cos-raw-2 (x orgx) ; [F F]
X (cond ((math-nearly-zerop-float x orgx) '(float 1 0))
X (calc-symbolic-mode (signal 'inexact-result nil))
X (t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
X (math-sin-series
X (math-add-float '(float 1 0)
X (math-mul-float xnegsqr '(float 5 -1)))
X 24 5 xnegsqr xnegsqr))))
X)
X
X(defun math-sin-series (sum nfac n x xnegsqr)
X (math-working "sin" sum)
X (let* ((nextx (math-mul-float x xnegsqr))
X (nextsum (math-add-float sum (math-div-float nextx
X (math-float nfac)))))
X (if (math-nearly-equal-float sum nextsum)
X sum
X (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
X (+ n 2) nextx xnegsqr)))
X)
X
X
X;;; Inverse sine, cosine, tangent.
X
X(defun math-arcsin (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2
X (math-from-radians (math-arcsin-raw (math-float x)))))
X ((and (eq (car x) 'sdev)
X (or (not (math-constp (nth 1 x)))
X (not (Math-lessp 1 (math-abs (nth 1 x))))))
X (math-make-sdev (math-arcsin (nth 1 x))
X (math-from-radians
X (math-div (nth 2 x)
X (math-sqrt
X (math-sub 1 (math-sqr (nth 1 x))))))))
X ((and (eq (car x) 'intv) (math-constp x))
X (math-sort-intv (nth 1 x)
X (math-arcsin (nth 2 x))
X (math-arcsin (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-arcsin x)))
X)
X(fset 'calcFunc-arcsin (symbol-function 'math-arcsin))
X
X(defun math-arccos (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2
X (math-from-radians (math-arccos-raw (math-float x)))))
X ((and (eq (car x) 'sdev)
X (or (not (math-constp (nth 1 x)))
X (not (Math-lessp 1 (math-abs (nth 1 x))))))
X (math-make-sdev (math-arccos (nth 1 x))
X (math-from-radians
X (math-div (nth 2 x)
X (math-sqrt
X (math-sub 1 (math-sqr (nth 1 x))))))))
X ((and (eq (car x) 'intv) (math-constp x))
X (math-sort-intv (nth 1 x)
X (math-arccos (nth 2 x))
X (math-arccos (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-arccos x)))
X)
X(fset 'calcFunc-arccos (symbol-function 'math-arccos))
X
X(defun math-arctan (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2
X (math-from-radians (math-arctan-raw (math-float x)))))
X ((eq (car x) 'sdev)
X (math-make-sdev (math-arctan (nth 1 x))
X (math-from-radians
X (math-div (nth 2 x)
X (math-add 1 (math-sqr (nth 1 x)))))))
X ((and (eq (car x) 'intv) (math-constp x))
X (math-sort-intv (nth 1 x)
X (math-arctan (nth 2 x))
X (math-arctan (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-arctan x)))
X)
X(fset 'calcFunc-arctan (symbol-function 'math-arctan))
X
X(defun math-arcsin-raw (x) ; [N N]
X (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x)))))
X (if (or (memq (car-safe x) '(cplx polar))
X (memq (car-safe a) '(cplx polar)))
X (math-with-extra-prec 2 ; use extra precision for difficult case
X (math-mul '(cplx 0 -1)
X (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a))))
X (math-arctan2-raw x a)))
X)
X
X(defun math-arccos-raw (x) ; [N N]
X (math-sub (math-pi-over-2) (math-arcsin-raw x))
X)
X
X(defun math-arctan-raw (x) ; [N N]
X (cond ((memq (car-safe x) '(cplx polar))
X (math-with-extra-prec 2 ; extra-extra
X (math-mul '(cplx 0 -1)
X (math-ln-raw (math-mul
X (math-add 1 (math-mul '(cplx 0 1) x))
X (math-sqrt-raw
X (math-div 1 (math-add
X 1 (math-sqr x)))))))))
X ((Math-integer-negp (nth 1 x))
X (math-neg-float (math-arctan-raw (math-neg-float x))))
X ((math-zerop x) x)
X ((math-equal-int x 1) (math-pi-over-4))
X ((math-equal-int x -1) (math-neg (math-pi-over-4)))
X (calc-symbolic-mode (signal 'inexact-result nil))
X ((math-lessp-float '(float 414214 -6) x) ; if x > sqrt(2) - 1, reduce
X (if (math-lessp-float '(float 1 0) x)
X (math-sub-float (math-mul-float (math-pi) '(float 5 -1))
X (math-arctan-raw (math-div-float '(float 1 0) x)))
X (math-sub-float (math-mul-float (math-pi) '(float 25 -2))
X (math-arctan-raw (math-div-float
X (math-sub-float '(float 1 0) x)
X (math-add-float '(float 1 0)
X x))))))
X (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x)))))
X)
X
X(defun math-arctan-series (sum n x xnegsqr)
X (math-working "arctan" sum)
X (let* ((nextx (math-mul-float x xnegsqr))
X (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
X (if (math-nearly-equal-float sum nextsum)
X sum
X (math-arctan-series nextsum (+ n 2) nextx xnegsqr)))
X)
X
X(defun math-arctan2 (y x) ; [F R R] [Public]
X (if (Math-anglep y)
X (if (Math-anglep x)
X (math-with-extra-prec 2
X (math-from-radians (math-arctan2-raw (math-float y)
X (math-float x))))
X (calc-record-why 'anglep x)
X (list 'calcFunc-arctan2 y x))
X (calc-record-why 'anglep y)
X (list 'calcFunc-arctan2 y x))
X)
X(fset 'calcFunc-arctan2 (symbol-function 'math-arctan2))
X
X(defun math-arctan2-raw (y x) ; [F R R]
X (cond ((math-zerop y)
X (if (math-negp x) (math-pi) 0))
X ((math-zerop x)
X (if (math-posp y)
X (math-pi-over-2)
X (math-neg (math-pi-over-2))))
X ((math-posp x)
X (math-arctan-raw (math-div-float y x)))
X ((math-posp y)
X (math-add-float (math-arctan-raw (math-div-float y x))
X (math-pi)))
X (t
X (math-sub-float (math-arctan-raw (math-div-float y x))
X (math-pi))))
X)
X
X(defun math-arc-sin-cos (x) ; [V N] [Public]
X (if (and (Math-vectorp x)
X (= (length x) 3))
X (math-arctan2 (nth 2 x) (nth 1 x))
X (math-reject-arg x "Two-element vector expected"))
X)
X(fset 'calcFunc-arcsincos (symbol-function 'math-arc-sin-cos))
X
X
X
X;;; Exponential function.
X
X(defun math-exp (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2 (math-exp-raw (math-float x))))
X ((eq (car-safe x) 'sdev)
X (let ((ex (math-exp (nth 1 x))))
X (math-make-sdev ex (math-mul (nth 2 x) ex))))
X ((eq (car-safe x) 'intv)
X (math-make-intv (nth 1 x) (math-exp (nth 2 x)) (math-exp (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-exp x)))
X)
X(fset 'calcFunc-exp (symbol-function 'math-exp))
X
X(defun math-exp-minus-1 (x) ; [N N] [Public]
X (cond ((math-zerop x) '(float 0 0))
X (calc-symbolic-mode (signal 'inexact-result nil))
X ((Math-numberp x)
X (math-with-extra-prec 2
X (let ((x (math-float x)))
X (if (and (eq (car x) 'float)
X (math-lessp-float x '(float 1 0))
X (math-lessp-float '(float -1 0) x))
X (math-exp-minus-1-raw x)
X (math-add (math-exp-raw x) -1)))))
X ((eq (car-safe x) 'sdev)
X (if (math-constp x)
X (let ((ex (math-exp-minus-1 (nth 1 x))))
X (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1))))
X (math-make-sdev (math-exp-minus-1 (nth 1 x))
X (math-mul (nth 2 x) (math-exp (nth 1 x))))))
X ((eq (car-safe x) 'intv)
X (math-make-intv (nth 1 x)
X (math-exp-minus-1 (nth 2 x))
X (math-exp-minus-1 (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-expm1 x)))
X)
X(fset 'calcFunc-expm1 (symbol-function 'math-exp-minus-1))
X
X(defun math-exp10 (x) ; [N N] [Public]
X (math-pow '(float 1 1) x)
X)
X(fset 'calcFunc-exp10 (symbol-function 'math-exp10))
X
X(defun math-exp-raw (x) ; [N N]
X (cond ((math-zerop x) '(float 1 0))
X (calc-symbolic-mode (signal 'inexact-result nil))
X ((eq (car x) 'cplx)
X (let ((expx (math-exp-raw (nth 1 x)))
X (sc (math-sin-cos-raw (nth 2 x))))
X (list 'cplx
X (math-mul-float expx (cdr sc))
X (math-mul-float expx (car sc)))))
X ((eq (car x) 'polar)
X (let ((xc (math-complex x)))
X (list 'polar
X (math-exp-raw (nth 1 x))
X (nth 2 x))))
X ((or (math-lessp-float '(float 5 -1) x)
X (math-lessp-float x '(float -5 -1)))
X (let* ((two-x (math-mul-float x '(float 2 0)))
X (hint (math-scale-int (nth 1 two-x) (nth 2 two-x)))
X (hfrac (math-sub-float x (math-mul-float (math-float hint)
X '(float 5 -1)))))
X (math-mul-float (math-ipow (math-sqrt-e) hint)
X (math-add-float '(float 1 0)
X (math-exp-minus-1-raw hfrac)))))
X (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x))))
X)
X
X(defun math-exp-minus-1-raw (x) ; [F F]
X (math-exp-series x 2 3 x x)
X)
X
X(defun math-exp-series (sum nfac n xpow x)
X (math-working "exp" sum)
X (let* ((nextx (math-mul-float xpow x))
X (nextsum (math-add-float sum (math-div-float nextx
X (math-float nfac)))))
X (if (math-nearly-equal-float sum nextsum)
X sum
X (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x)))
X)
X
X
X
X;;; Logarithms.
X
X(defun math-ln (x) ; [N N] [Public]
X (cond ((math-zerop x)
X (math-reject-arg x "Logarithm of zero"))
X ((Math-numberp x)
X (math-with-extra-prec 2 (math-ln-raw (math-float x))))
X ((and (eq (car-safe x) 'sdev)
X (or (not (math-constp (nth 1 x)))
X (math-posp (nth 1 x))))
X (math-make-sdev (math-ln (nth 1 x))
X (math-div (nth 2 x) (nth 1 x))))
X ((and (eq (car-safe x) 'intv) (Math-posp (nth 2 x)))
X (math-make-intv (nth 1 x) (math-ln (nth 2 x)) (math-ln (nth 3 x))))
X ((equal x '(var e var-e))
X 1)
X ((and (eq (car-safe x) '^)
X (equal (nth 1 x) '(var e var-e)))
X (nth 2 x))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-ln x)))
X)
X(fset 'calcFunc-ln (symbol-function 'math-ln))
X
X(defun math-log10 (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2
X (let ((xf (math-float x)))
X (if (eq (nth 1 xf) 0)
X (math-reject-arg x "Logarithm of zero"))
X (if (Math-integer-posp (nth 1 xf))
X (if (eq (nth 1 xf) 1) ; log10(1*10^n) = n
X (math-float (nth 2 xf))
X (let ((xdigs (1- (math-numdigs (nth 1 xf)))))
X (math-add-float
X (math-div-float (math-ln-raw-2
X (list 'float (nth 1 xf) (- xdigs)))
X (math-ln-10))
X (math-float (+ (nth 2 xf) xdigs)))))
X (math-div (math-ln xf) (math-ln-10))))))
X ((and (eq (car-safe x) 'sdev)
X (or (not (math-constp (nth 1 x)))
X (math-posp (nth 1 x))))
X (math-make-sdev (math-log10 (nth 1 x))
X (math-div (nth 2 x)
X (math-mul (nth 1 x) (math-ln 10)))))
X ((and (eq (car-safe x) 'intv) (Math-posp (nth 2 x)))
X (math-make-intv (nth 1 x)
X (math-log10 (nth 2 x))
X (math-log10 (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-log10 x)))
X)
X(fset 'calcFunc-log10 (symbol-function 'math-log10))
X
X(defun calcFunc-pow10 (x)
X (math-pow '(float 1 1) x)
X)
X
X(defun math-log (x b) ; [N N N] [Public]
X (cond ((or (eq b 10) (equal b '(float 1 1)))
X (math-log10 x))
X ((math-zerop x)
X (math-reject-arg x "Logarithm of zero"))
X ((math-zerop b)
X (math-reject-arg b "Logarithm of zero"))
X ((and (Math-numberp x) (Math-numberp b))
X (math-with-extra-prec 2
X (math-div (math-ln-raw (math-float x))
X (math-log-base-raw b))))
X ((and (eq (car-safe x) 'sdev)
X (or (not (math-constp (nth 1 x)))
X (math-posp (nth 1 x)))
X (Math-numberp b))
X (math-make-sdev (math-log (nth 1 x) b)
X (math-div (nth 2 x)
X (math-mul (nth 1 x)
X (math-log-base-raw b)))))
X ((and (eq (car-safe x) 'intv) (Math-posp (nth 2 x)))
X (math-make-intv (nth 1 x)
X (math-log (nth 2 x) b)
X (math-log (nth 3 x) b)))
X (t (if (Math-numberp b)
X (calc-record-why 'numberp x)
X (calc-record-why 'numberp b))
X (list 'calcFunc-log x b)))
X)
X
X(defun calcFunc-log (x &optional b)
X (if b
X (if (or (eq b 10) (equal b '(float 1 1)))
X (math-normalize (list 'calcFunc-log10 x))
X (if (equal b '(var e var-e))
X (math-normalize (list 'calcFunc-ln x))
X (math-log x b)))
X (math-normalize (list 'calcFunc-ln x)))
X)
X(defun calcFunc-ilog (x &optional b)
X (if b
X (if (equal b '(var e var-e))
X (math-normalize (list 'calcFunc-exp x))
X (math-pow b x))
X (math-normalize (list 'calcFunc-exp x)))
X)
X
X
X(defun math-log-base-raw (b) ; [N N]
X (if (not (and (equal (car math-log-base-cache) b)
X (eq (nth 1 math-log-base-cache) calc-internal-prec)))
X (setq math-log-base-cache (list b calc-internal-prec
X (math-ln-raw (math-float b)))))
X (nth 2 math-log-base-cache)
X)
X(setq math-log-base-cache nil)
X
X(defun math-ln-plus-1 (x) ; [N N] [Public]
X (cond ((Math-equal-int x -1) (math-reject-arg x "Logarithm of zero"))
X ((math-zerop x) '(float 0 0))
X (calc-symbolic-mode (signal 'inexact-result nil))
X ((Math-numberp x)
X (math-with-extra-prec 2
X (let ((x (math-float x)))
X (if (and (eq (car x) 'float)
X (math-lessp-float x '(float 5 -1))
X (math-lessp-float '(float -5 -1) x))
X (math-ln-plus-1-raw x)
X (math-ln-raw (math-add-float x '(float 1 0)))))))
X ((and (eq (car-safe x) 'sdev)
X (or (not (math-constp (nth 1 x)))
X (Math-lessp -1 (nth 1 x))))
X (math-make-sdev (math-ln-plus-1 (nth 1 x))
X (math-div (nth 2 x) (math-add (nth 1 x) 1))))
X ((and (eq (car-safe x) 'intv) (Math-posp (nth 2 x)))
X (math-make-intv (nth 1 x)
X (math-ln-plus-1 (nth 2 x))
X (math-ln-plus-1 (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-lnp1 x)))
X)
X(fset 'calcFunc-lnp1 (symbol-function 'math-ln-plus-1))
X
X(defun math-ln-raw (x) ; [N N] --- must be float format!
X (cond ((eq (car-safe x) 'cplx)
X (list 'cplx
X (math-mul-float (math-ln-raw
X (math-add-float (math-sqr-float (nth 1 x))
X (math-sqr-float (nth 2 x))))
X '(float 5 -1))
X (math-arctan2-raw (nth 2 x) (nth 1 x))))
X ((eq (car x) 'polar)
X (math-polar (list 'cplx
X (math-ln-raw (nth 1 x))
X (nth 2 x))))
X ((Math-equal-int x 1)
X '(float 0 0))
X (calc-symbolic-mode (signal 'inexact-result nil))
X ((math-posp (nth 1 x)) ; positive and real
X (let ((xdigs (1- (math-numdigs (nth 1 x)))))
X (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs)))
X (math-mul-float (math-float (+ (nth 2 x) xdigs))
X (math-ln-10)))))
X ((math-zerop x)
X (error "Logarithm of zero"))
X ((eq calc-complex-mode 'polar) ; negative and real
X (math-polar
X (list 'cplx ; negative and real
X (math-ln-raw (math-neg-float x))
X (math-pi))))
X (t (list 'cplx ; negative and real
X (math-ln-raw (math-neg-float x))
X (math-pi))))
X)
X
X(defun math-ln-raw-2 (x) ; [F F]
X (cond ((math-lessp-float '(float 14 -1) x)
X (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1)))
X (math-ln-2)))
X (t ; now .7 < x <= 1.4
X (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
X (math-add-float x '(float 1 0))))))
X)
X
X(defun math-ln-raw-3 (x) ; [F F]
X (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
X '(float 2 0))
X)
X
X;;; Compute ln((1+x)/(1-x))
X(defun math-ln-raw-series (sum n x xsqr)
X (math-working "log" sum)
X (let* ((nextx (math-mul-float x xsqr))
X (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
X (if (math-nearly-equal-float sum nextsum)
X sum
X (math-ln-raw-series nextsum (+ n 2) nextx xsqr)))
X)
X
X(defun math-ln-plus-1-raw (x)
X (math-lnp1-series x 2 x (math-neg x))
X)
X
X(defun math-lnp1-series (sum n xpow x)
X (math-working "lnp1" sum)
X (let* ((nextx (math-mul-float xpow x))
X (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
X (if (math-nearly-equal-float sum nextsum)
X sum
X (math-lnp1-series nextsum (1+ n) nextx x)))
X)
X
X(math-defcache math-ln-10 (float (bigpos 018 684 045 994 092 585 302 2) -21)
X (math-ln-raw-2 '(float 1 1)))
X
X(math-defcache math-ln-2 (float (bigpos 417 309 945 559 180 147 693) -21)
X (math-ln-raw-3 (math-float '(frac 1 3))))
X
X
X
X;;; Hyperbolic functions.
X
X(defun math-sinh (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2
X (let ((expx (math-exp-raw (math-float x))))
X (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1)))))
X ((eq (car-safe x) 'sdev)
X (math-make-sdev (math-sinh (nth 1 x))
X (math-mul (nth 2 x) (math-cosh (nth 1 x)))))
X ((and (eq (car x) 'intv) (math-constp x))
X (math-sort-intv (nth 1 x)
X (math-sinh (nth 2 x))
X (math-sinh (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-sinh x)))
X)
X(fset 'calcFunc-sinh (symbol-function 'math-sinh))
X
X(defun math-cosh (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2
X (let ((expx (math-exp-raw (math-float x))))
X (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1)))))
X ((eq (car-safe x) 'sdev)
X (math-make-sdev (math-cosh (nth 1 x))
X (math-mul (nth 2 x)
X (math-abs (math-sinh (nth 1 x))))))
X ((and (eq (car x) 'intv) (math-constp x))
X (setq x (math-abs x))
X (math-sort-intv (nth 1 x)
X (math-cosh (nth 2 x))
X (math-cosh (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-cosh x)))
X)
X(fset 'calcFunc-cosh (symbol-function 'math-cosh))
X
X(defun math-tanh (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2
X (let* ((expx (math-exp (math-float x)))
X (expmx (math-div 1 expx)))
X (math-div (math-sub expx expmx)
X (math-add expx expmx)))))
X ((eq (car-safe x) 'sdev)
X (math-make-sdev (math-tanh (nth 1 x))
X (math-div (nth 2 x)
X (math-sqr (math-cosh (nth 1 x))))))
X ((and (eq (car x) 'intv) (math-constp x))
X (math-sort-intv (nth 1 x)
X (math-tanh (nth 2 x))
X (math-tanh (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-tanh x)))
X)
X(fset 'calcFunc-tanh (symbol-function 'math-tanh))
X
X(defun math-arcsinh (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2
X (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x)
X '(float 1 0)))))))
X ((eq (car-safe x) 'sdev)
X (math-make-sdev (math-arcsinh (nth 1 x))
X (math-div (nth 2 x)
X (math-sqrt
X (math-add (math-sqr (nth 1 x)) 1)))))
X ((and (eq (car x) 'intv) (math-constp x))
X (math-sort-intv (nth 1 x)
X (math-arcsinh (nth 2 x))
X (math-arcsinh (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-arcsinh x)))
X)
X(fset 'calcFunc-arcsinh (symbol-function 'math-arcsinh))
X
X(defun math-arccosh (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2
X (if (or t ; need to do this even in the real case!
X (memq (car-safe x) '(cplx polar)))
X (let ((xp1 (math-add 1 x))) ; this gets the branch cuts right
X (math-ln-raw
X (math-add x (math-mul xp1
X (math-sqrt-raw (math-div (math-sub
X x
X '(float 1 0))
X xp1))))))
X (math-ln-raw
X (math-add x (math-sqrt-raw (math-add (math-sqr x)
X '(float -1 0))))))))
X ((and (eq (car-safe x) 'sdev)
X (or (not (math-constp (nth 1 x)))
X (not (Math-lessp (nth 1 x) 1))))
X (math-make-sdev (math-arccosh (nth 1 x))
X (math-div (nth 2 x)
X (math-sqrt
X (math-add (math-sqr (nth 1 x)) -1)))))
X ((and (eq (car x) 'intv) (math-constp x))
X (math-sort-intv (nth 1 x)
X (math-arccosh (nth 2 x))
X (math-arccosh (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-arccosh x)))
X)
X(fset 'calcFunc-arccosh (symbol-function 'math-arccosh))
X
X(defun math-arctanh (x) ; [N N] [Public]
X (cond ((Math-numberp x)
X (math-with-extra-prec 2
X (if (memq (car-safe x) '(cplx polar))
X (math-ln-raw
X (math-mul (math-add 1 x)
X (math-sqrt-raw
X (math-div '(float 1 0) (math-sub 1 (math-sqr x))))))
X (math-mul (math-ln-raw (math-div (math-add '(float 1 0) x)
X (math-sub 1 x)))
X '(float 5 -1)))))
X ((and (eq (car-safe x) 'sdev)
X (or (not (math-constp (nth 1 x)))
X (Math-lessp (math-abs (nth 1 x)) 1)))
X (math-make-sdev (math-arctanh (nth 1 x))
X (math-div (nth 2 x)
X (math-sub 1 (math-sqr (nth 1 x))))))
X ((and (eq (car x) 'intv) (math-constp x))
X (math-sort-intv (nth 1 x)
X (math-arctanh (nth 2 x))
X (math-arctanh (nth 3 x))))
X (t (calc-record-why 'numberp x)
X (list 'calcFunc-arctanh x)))
X)
X(fset 'calcFunc-arctanh (symbol-function 'math-arctanh))
X
X
X;;; Convert A from HMS or degrees to radians.
X(defun math-deg-to-rad (a) ; [R R] [Public]
X (cond ((or (Math-numberp a)
X (eq (car a) 'intv))
X (math-mul a (math-pi-over-180)))
X ((eq (car a) 'hms)
X (math-from-hms a 'rad))
X ((eq (car a) 'sdev)
X (math-make-sdev (math-deg-to-rad (nth 1 a))
X (math-deg-to-rad (nth 2 a))))
X (t (list 'calcFunc-rad a)))
X)
X(fset 'calcFunc-rad (symbol-function 'math-deg-to-rad))
X
X;;; Convert A from HMS or radians to degrees.
X(defun math-rad-to-deg (a) ; [R R] [Public]
X (cond ((or (Math-numberp a)
X (eq (car a) 'intv))
X (math-div a (math-pi-over-180)))
X ((eq (car a) 'hms)
X (math-from-hms a 'deg))
X ((eq (car a) 'sdev)
X (math-make-sdev (math-rad-to-deg (nth 1 a))
X (math-rad-to-deg (nth 2 a))))
X (t (list 'calcFunc-deg a)))
X)
X(fset 'calcFunc-deg (symbol-function 'math-rad-to-deg))
X
X
X
X
X;;;; Number theory.
X
X(defun calcFunc-idiv (a b) ; [I I I] [Public]
X (cond ((and (Math-natnump a) (Math-natnump b) (not (eq b 0)))
X (math-quotient a b))
X ((Math-realp a)
X (if (Math-realp b)
X (let ((calc-prefer-frac t))
X (math-floor (math-div a b)))
X (math-reject-arg b 'realp)))
X ((eq (car-safe a) 'hms)
X (if (eq (car-safe b) 'hms)
X (let ((calc-prefer-frac t))
X (math-floor (math-div a b)))
X (math-reject-arg b 'hmsp)))
X ((and (or (eq (car-safe a) 'intv) (Math-realp a))
X (or (eq (car-safe b) 'intv) (Math-realp b)))
X (math-floor (math-div a b)))
X (t (math-reject-arg a 'anglep)))
X)
X
X(defun calcFunc-fdiv (a b) ; [R I I] [Public]
X (if (Math-num-integerp a)
X (if (Math-num-integerp b)
X (if (Math-zerop b)
X (math-reject-arg a "Division by zero")
X (math-make-frac (math-trunc a) (math-trunc b)))
X (math-reject-arg b 'integerp))
X (math-reject-arg a 'integerp))
X)
X
X(defun math-lcm (a b)
X (let ((g (math-gcd a b)))
X (if (Math-numberp g)
X (math-div (math-mul a b) g)
X (list 'calcFunc-lcm a b)))
X)
X(fset 'calcFunc-lcm (symbol-function 'math-lcm))
X
X(defun math-extended-gcd (a b) ; Knuth section 4.5.2
X (cond
X ((not (Math-integerp a))
X (if (Math-messy-integerp a)
X (math-extended-gcd (math-trunc a) b)
X (calc-record-why 'integerp a)
X (list 'calcFunc-egcd a b)))
X ((not (Math-integerp b))
X (if (Math-messy-integerp b)
X (math-extended-gcd a (math-trunc b))
X (calc-record-why 'integerp b)
X (list 'calcFunc-egcd a b)))
X (t
X (let ((u1 1) (u2 0) (u3 a)
X (v1 0) (v2 1) (v3 b)
X t1 t2 q)
X (while (not (eq v3 0))
X (setq q (math-idivmod u3 v3)
X t1 (math-sub u1 (math-mul v1 (car q)))
X t2 (math-sub u2 (math-mul v2 (car q)))
X u1 v1 u2 v2 u3 v3
X v1 t1 v2 t2 v3 (cdr q)))
X (list 'vec u3 u1 u2))))
X)
X(fset 'calcFunc-egcd (symbol-function 'math-extended-gcd))
X
X
X;;; Factorial and related functions.
X
X(defun math-factorial (n) ; [I I] [F F] [Public]
X (let (temp)
X (cond ((Math-integer-negp n) (list 'calcFunc-fact n))
X ((Math-zerop n) 1)
X ((integerp n) (math-factorial-iter (1- n) 2 1))
X ((and (math-messy-integerp n)
X (Math-lessp (setq temp (math-trunc n)) 100))
X (if (natnump temp)
X (math-with-extra-prec 1
X (math-factorial-iter (1- temp) 2 '(float 1 0)))
X (list 'calcFunc-fact max)))
X ((math-realp n)
X (math-with-extra-prec 3
X (math-gammap1-raw (math-float n))))
X (t (calc-record-why 'realp n)
X (list 'calcFunc-fact n))))
X)
X(fset 'calcFunc-fact (symbol-function 'math-factorial))
X
X(defun math-factorial-iter (count n f)
X (if (= (% n 5) 1)
X (math-working (format "factorial(%d)" (1- n)) f))
X (if (> count 0)
X (math-factorial-iter (1- count) (1+ n) (math-mul n f))
X f)
X)
X
X(math-defcache math-sqrt-two-pi nil
X (math-sqrt (math-two-pi)))
X
X(defun math-gammap1-raw (x) ; compute gamma(1 + x)
X (cond ((math-lessp-float x '(float 1 1))
X (if (math-lessp-float x '(float -10 0))
X (setq x (math-neg-float
X (math-div-float
X (math-pi)
X (math-mul-float (math-gammap1-raw
X (math-add-float (math-neg-float x)
X '(float -1 0)))
X (math-sin-raw
X (math-mul (math-pi) x)))))))
X (let ((xplus1 (math-add-float x '(float 1 0))))
X (math-div-float (math-gammap1-raw xplus1) xplus1)))
X (t ; x now >= 10.0
X (let ((xinv (math-div 1 x))
X (lnx (math-ln-raw x)))
X (math-mul (math-sqrt-two-pi)
X (math-exp-raw
X (math-gamma-series
X (math-sub (math-mul (math-add x '(float 5 -1))
X lnx)
X x)
X xinv
X (math-sqr xinv)
X 2))))))
X)
X
X(defun calcFunc-gamma (x)
X (calcFunc-fact (math-add x -1))
X)
X
X(defun math-gamma-series (sum x xinvsqr n)
X (math-working "gamma" sum)
X (let* ((bn (math-bernoulli-number n)) ; this will always be a "frac" form.
X (next (math-add-float
X sum
X (math-mul-float (math-div-float (math-float (nth 1 bn))
X (math-float (* (nth 2 bn)
X (* n (1- n)))))
X x))))
X (if (math-nearly-equal-float sum next)
X next
X (if (= n 24)
X (progn
X (calc-record-why "Gamma computation stopped early, not all digits may be valid")
X next)
X (math-gamma-series next (math-mul-float x xinvsqr) xinvsqr (+ n 2)))))
X)
X
X(defun math-bernoulli-number (n)
X (if (= n 1)
X '(frac -1 2)
X (if (= (% n 2) 1)
X 0
X (aref '[ 1 (frac 1 6) (frac -1 30) (frac 1 42) (frac -1 30)
X (frac 5 66) (frac -691 2730) (frac 7 6) (frac -3617 510)
X (frac 43867 798) (frac -174611 330) (frac 854513 138)
X (frac (bigneg 91 364 236) 2730) ]
X (/ n 2))))
X)
X;;; To come up with more, we could use this rule:
X;;; Bn = n! bn
X;;; bn = - sum_k=0^n-1 bk / (n-k+1)!
X
X(defun math-double-factorial (n) ; [I I] [F F] [Public]
X (cond ((Math-integer-negp n) (list 'calcFunc-dfact n))
X ((Math-zerop n) 1)
X ((integerp n) (math-double-factorial-iter n (+ 2 (% n 2)) 1))
X ((math-messy-integerp n)
X (let ((temp (math-trunc n)))
X (if (natnump temp)
X (math-with-extra-prec 1
X (math-double-factorial-iter temp (+ 2 (% temp 2))
X '(float 1 0)))
X (list 'calcFunc-dfact max))))
X (t (calc-record-why 'natnump n)
X (list 'calcFunc-dfact n)))
X)
X(fset 'calcFunc-dfact (symbol-function 'math-double-factorial))
X
X(defun math-double-factorial-iter (max n f)
X (if (< (% n 10) 2)
X (math-working (format "dfact(%d)" (- n 2)) f))
X (if (<= n max)
X (math-double-factorial-iter max (+ n 2) (math-mul n f))
X f)
X)
X
X(defun math-permutations (n m) ; [I I I] [F F F] [Public]
X (cond ((and (integerp n) (integerp m) (<= m n) (>= m 0))
X (math-factorial-iter n (1+ (- n m)) 1))
X ((or (not (math-num-integerp n))
X (not (math-num-integerp m)))
X (or (math-realp n) (math-reject-arg n 'realp))
X (or (math-realp m) (math-reject-arg m 'realp))
X (and (math-num-integerp n) (math-negp n) (math-reject-arg n 'range))
X (and (math-num-integerp m) (math-negp m) (math-reject-arg m 'range))
X (math-div (math-factorial n) (math-factorial m)))
X (t
X (let ((tn (math-trunc n))
X (tm (math-trunc m)))
X (or (integerp tn) (math-reject-arg tn 'fixnump))
X (or (integerp tm) (math-reject-arg tm 'fixnump))
X (or (and (<= tm tn) (>= tm 0)) (math-reject-arg tm 'range))
X (math-with-extra-prec 1
X (math-factorial-iter tm (1+ (- tn tm)) '(float 1 0))))))
X)
X(fset 'calcFunc-perm (symbol-function 'math-permutations))
X
X(defun math-choose (n m) ; [I I I] [F F F] [Public]
X (cond ((and (integerp n) (integerp m) (<= m n) (>= m 0))
X (math-choose-iter m n 1 1))
X ((not (math-realp n))
X (math-reject-arg n 'realp))
X ((not (math-realp m))
X (math-reject-arg m 'realp))
X ((not (math-num-integerp m))
X (if (and (math-num-integerp n) (math-negp n))
X (list 'calcFunc-choose n m)
X (math-div (math-factorial (math-float n))
X (math-mul (math-factorial m)
X (math-factorial (math-sub n m))))))
X ((math-negp m) 0)
X ((math-negp n)
X (let ((val (math-choose (math-add (math-add n m) -1) m)))
X (if (math-evenp (math-trunc m))
X val
X (math-neg val))))
X ((and (math-num-integerp n)
X (Math-lessp n m))
X 0)
X (t
X (let ((tm (math-trunc m)))
X (or (integerp tm) (math-reject-arg tm 'fixnump))
X (if (> tm 100)
X (math-div (math-factorial (math-float n))
X (math-mul (math-factorial (math-float m))
X (math-factorial (math-float
X (math-sub n m)))))
X (math-with-extra-prec 1
X (math-choose-float-iter tm n 1 '(float 1 0)))))))
X)
X(fset 'calcFunc-choose (symbol-function 'math-choose))
X
X(defun math-choose-iter (m n i c)
X (if (= (% i 5) 1)
X (math-working (format "choose(%d)" (1- i)) c))
X (if (<= i m)
X (math-choose-iter m (1- n) (1+ i)
X (math-quotient (math-mul c n) i))
X c)
X)
X
X(defun math-choose-float-iter (count n i c)
X (if (= (% i 5) 1)
X (math-working (format "choose(%d)" (1- i)) c))
X (if (> count 0)
X (math-choose-float-iter (1- count) (math-sub n 1) (1+ i)
X (math-div (math-mul c n) i))
X c)
X)
X
X
X;;; Produce a random digit in the range 0..999.
X;;; Avoid various pitfalls that may lurk in the built-in (random) function!
X(defun math-random-digit ()
X (prog1
X (% (lsh (random math-first-random-flag) -4) 1000)
X (setq math-first-random-flag nil))
X)
X(setq math-first-random-flag t)
X
X;;; Produce an N-digit random integer.
X(defun math-random-digits (n)
X (cond ((<= n 6)
X (math-scale-right (+ (* (math-random-digit) 1000) (math-random-digit))
X (- 6 n)))
X (t (let* ((slop (% (- 900003 n) 3))
X (i (/ (+ n slop) 3))
X (digs nil))
X (while (> i 0)
X (setq digs (cons (math-random-digit) digs)
X i (1- i)))
X (math-normalize (math-scale-right (cons 'bigpos digs)
X slop)))))
X)
X
X;;; Produce a uniformly-distributed random float 0 <= N < 1.
X(defun math-random-float ()
X (math-make-float (math-random-digits calc-internal-prec)
X (- calc-internal-prec))
X)
X
X;;; Produce a Gaussian-distributed random float with mean=0, sigma=1.
X(defun math-gaussian-float ()
X (math-with-extra-prec 2
X (if (and math-gaussian-cache
X (= (car math-gaussian-cache) calc-internal-prec))
X (prog1
X (cdr math-gaussian-cache)
X (setq math-gaussian-cache nil))
X (let* ((v1 (math-add (math-mul (math-random-float) 2) -1))
X (v2 (math-add (math-mul (math-random-float) 2) -1))
X (r (math-add (math-sqr v1) (math-sqr v2))))
X (while (or (not (Math-lessp r 1)) (math-zerop r))
X (setq v1 (math-add (math-mul (math-random-float) 2) -1)
X v2 (math-add (math-mul (math-random-float) 2) -1)
X r (math-add (math-sqr v1) (math-sqr v2))))
X (let ((fac (math-sqrt (math-mul (math-div (math-ln r) r) -2))))
X (setq math-gaussian-cache (cons calc-internal-prec
X (math-mul v1 fac)))
X (math-mul v2 fac)))))
X)
X(setq math-gaussian-cache nil)
X
X;;; Produce a random integer or real 0 <= N < MAX.
X(defun math-random (max)
X (cond ((Math-zerop max)
X (math-gaussian-float))
X ((Math-integerp max)
X (let* ((digs (math-numdigs max))
X (r (math-random-digits (+ digs 3))))
X (math-mod r max)))
X ((Math-realp max)
X (math-mul (math-random-float) max))
X ((and (eq (car max) 'intv) (math-constp max)
X (Math-lessp (nth 2 max) (nth 3 max)))
X (if (math-floatp max)
X (let ((val (math-add (math-mul (math-random-float)
X (math-sub (nth 3 max) (nth 2 max)))
X (nth 2 max))))
X (if (or (and (memq (nth 1 max) '(0 1)) ; almost not worth
X (math-equal val (nth 2 max))) ; checking!
X (and (memq (nth 1 max) '(0 2))
X (math-equal val (nth 3 max))))
X (math-random max)
X val))
X (let ((lo (if (memq (nth 1 max) '(0 1))
X (math-add (nth 2 max) 1) (nth 2 max)))
X (hi (if (memq (nth 1 max) '(1 3))
X (math-add (nth 3 max) 1) (nth 3 max))))
X (if (Math-lessp lo hi)
X (math-add (math-random (math-sub hi lo)) lo)
X (math-reject-arg max "Empty interval")))))
X (t (math-reject-arg max 'realp)))
X)
X(fset 'calcFunc-random (symbol-function 'math-random))
X
X
X;;; Check if the integer N is prime. [X I]
X;;; Return (nil) if non-prime,
X;;; (nil N) if non-prime with known factor N,
X;;; (nil unknown) if non-prime with no known factors,
X;;; (t) if prime,
X;;; (maybe N P) if probably prime (after N iters with probability P%)
X(defun math-prime-test (n iters)
X (if (and (Math-vectorp n) (cdr n))
X (setq n (nth (1- (length n)) n)))
X (if (Math-messy-integerp n)
X (setq n (math-trunc n)))
X (let ((res))
X (while (> iters 0)
X (setq res
X (cond ((and (integerp n) (<= n 5003))
X (list (= (math-next-small-prime n) n)))
X ((not (Math-integerp n))
X (error "Argument must be an integer"))
X ((Math-integer-negp n)
X '(nil))
X ((Math-natnum-lessp n '(bigpos 0 0 8))
X (setq n (math-fixnum n))
X (let ((i -1) v)
X (while (and (> (% n (setq v (aref math-primes-table
X (setq i (1+ i)))))
X 0)
X (< (* v v) n)))
X (if (= (% n v) 0)
X (list nil v)
X '(t))))
X ((not (equal n (car math-prime-test-cache)))
X (cond ((= (% (nth 1 n) 2) 0) '(nil 2))
X ((= (% (nth 1 n) 5) 0) '(nil 5))
X (t (let ((dig (cdr n)) (sum 0))
X (while dig
X (if (cdr dig)
X (setq sum (% (+ (+ sum (car dig))
X (* (nth 1 dig) 1000))
X 111111)
X dig (cdr (cdr dig)))
X (setq sum (% (+ sum (car dig)) 111111)
X dig nil)))
X (cond ((= (% sum 3) 0) '(nil 3))
X ((= (% sum 7) 0) '(nil 7))
X ((= (% sum 11) 0) '(nil 11))
X ((= (% sum 13) 0) '(nil 13))
X ((= (% sum 37) 0) '(nil 37))
X (t
X (setq math-prime-test-cache-k 1
X math-prime-test-cache-q
X (math-div2 n)
X math-prime-test-cache-nm1
X (math-add n -1))
X (while (math-evenp
X math-prime-test-cache-q)
X (setq math-prime-test-cache-k
X (1+ math-prime-test-cache-k)
X math-prime-test-cache-q
X (math-div2
X math-prime-test-cache-q)))
X (setq iters (1+ iters))
X (list 'maybe
X 0
X (math-sub
X 100
X (math-div
X '(float 232 0)
X (math-numdigs n))))))))))
X ((not (eq (car (nth 1 math-prime-test-cache)) 'maybe))
X (nth 1 math-prime-test-cache))
X (t ; Fermat step
X (let* ((x (math-add (math-random (math-add n -2)) 2))
X (y (math-pow-mod x math-prime-test-cache-q n))
X (j 0))
X (while (and (not (eq y 1))
X (not (equal y math-prime-test-cache-nm1))
X (< (setq j (1+ j)) math-prime-test-cache-k))
X (setq y (math-mod (math-mul y y) n)))
X (if (or (equal y math-prime-test-cache-nm1)
X (and (eq y 1) (eq j 0)))
X (list 'maybe
X (1+ (nth 1 (nth 1 math-prime-test-cache)))
X (math-mul (nth 2 (nth 1 math-prime-test-cache))
X '(float 25 -2)))
X '(nil unknown))))))
X (setq math-prime-test-cache (list n res)
X iters (if (eq (car res) 'maybe)
X (1- iters)
X 0)))
X res)
X)
X(defvar math-prime-test-cache '(-1))
X
X;;; Theory: summing base-10^6 digits modulo 111111 is "casting out 999999s".
X;;; Initial probability that N is prime is 1/ln(N) = log10(e)/log10(N).
X;;; After culling [2,3,5,7,11,13,37], probability of primality is 5.36 x more.
X;;; Initial reported probability of non-primality is thus 100% - this.
X;;; Each Fermat step multiplies this probability by 25%.
X;;; The Fermat step is algorithm P from Knuth section 4.5.4.
X
X
X(defun math-prime-factors (n)
X (setq math-prime-factors-finished t)
X (if (Math-messy-integerp n)
X (setq n (math-trunc n)))
X (if (and (Math-natnump n) (Math-natnum-lessp 2 n))
X (let (factors res p (i 0))
X (while (and (not (eq n 1))
X (< i (length math-primes-table)))
X (setq p (aref math-primes-table i))
X (while (eq (cdr (setq res (cond ((eq n p) (cons 1 0))
X ((eq n 1) (cons 0 1))
X ((consp n) (math-idivmod n p))
X (t (cons (/ n p) (% n p))))))
X 0)
X (math-working "factor" p)
X (setq factors (nconc factors (list p))
X n (car res)))
X (or (eq n 1)
X (Math-natnum-lessp p (car res))
X (setq factors (nconc factors (list n))
X n 1))
X (setq i (1+ i)))
X (or (setq math-prime-factors-finished (eq n 1))
X (setq factors (nconc factors (list n))))
X (cons 'vec factors))
X (calc-record-why 'integerp n)
X (list 'calcFunc-prfac n))
X)
X(fset 'calcFunc-prfac (symbol-function 'math-prime-factors))
X
X(defun math-totient (n)
X (if (Math-messy-integerp n)
X (setq n (math-trunc n)))
X (if (Math-natnump n)
X (if (Math-natnum-lessp n 2)
X (if (Math-negp n)
X (math-totient (math-abs n))
X n)
X (let ((factors (cdr (math-prime-factors n)))
X p)
X (if math-prime-factors-finished
X (progn
X (while factors
X (setq p (car factors)
X n (math-mul (math-div n p) (math-add p -1)))
X (while (equal p (car factors))
X (setq factors (cdr factors))))
X n)
X (calc-record-why "Number too big to factor" n)
X (list 'calcFunc-totient n))))
X (calc-record-why 'natnump n)
X (list 'calcFunc-totient n))
X)
X(fset 'calcFunc-totient (symbol-function 'math-totient))
X
X(defun math-moebius (n)
X (if (Math-messy-integerp n)
X (setq n (math-trunc n)))
X (if (and (Math-natnump n) (not (eq n 0)))
X (if (Math-natnum-lessp n 2)
X (if (Math-negp n)
X (math-moebius (math-abs n))
X 1)
X (let ((factors (cdr (math-prime-factors n)))
X (mu 1))
X (if math-prime-factors-finished
X (progn
X (while factors
X (setq mu (if (equal (car factors) (nth 1 factors))
X 0 (math-neg mu))
X factors (cdr factors)))
X mu)
X (calc-record-why "Number too big to factor" n)
X (list 'calcFunc-moebius n))))
X (calc-record-why 'natnump n)
X (list 'calcFunc-moebius n))
X)
X(fset 'calcFunc-moebius (symbol-function 'math-moebius))
X
X
X(defun math-next-prime (n iters)
X (if (Math-integerp n)
X (if (Math-integer-negp n)
X 2
X (if (and (integerp n) (< n 5003))
X (math-next-small-prime (1+ n))
X (if (math-evenp n)
X (setq n (math-add n -1)))
X (let (res)
X (while (not (car (setq res (math-prime-test
X (setq n (math-add n 2)) iters)))))
X (if (and calc-verbose-nextprime
X (eq (car res) 'maybe))
X (calc-report-prime-test res)))
X n))
X (if (Math-realp n)
X (math-next-prime (math-trunc n) iters)
X (math-reject-arg n 'integerp)))
X)
X(fset 'calcFunc-nextprime (symbol-function 'math-next-prime))
X(setq calc-verbose-nextprime nil)
X
X(defun math-prev-prime (n iters)
X (if (Math-integerp n)
X (if (Math-lessp n 4)
X 2
X (if (math-evenp n)
X (setq n (math-add n 1)))
X (let (res)
X (while (not (car (setq res (math-prime-test
X (setq n (math-add n -2)) iters)))))
X (if (and calc-verbose-nextprime
X (eq (car res) 'maybe))
X (calc-report-prime-test res)))
X n)
X (if (Math-realp n)
X (math-prev-prime (math-ceiling n) iters)
X (math-reject-arg n 'integerp)))
X)
X(fset 'calcFunc-prevprime (symbol-function 'math-prev-prime))
X
X(defun math-next-small-prime (n)
X (if (and (integerp n) (> n 2))
X (let ((lo -1)
X (hi (length math-primes-table))
X mid)
X (while (> (- hi lo) 1)
X (if (> n (aref math-primes-table
X (setq mid (ash (+ lo hi) -1))))
X (setq lo mid)
X (setq hi mid)))
X (aref math-primes-table hi))
X 2)
X)
X
X(defconst math-primes-table
X [2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89
X 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181
X 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277
X 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383
X 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487
X 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601
X 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709
X 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827
X 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947
X 953 967 971 977 983 991 997 1009 1013 1019 1021 1031 1033 1039 1049
X 1051 1061 1063 1069 1087 1091 1093 1097 1103 1109 1117 1123 1129 1151
SHAR_EOF
echo "End of part 8"
echo "File calc-ext.el is continued in part 9"
echo "9" > s2_seq_.tmp
exit 0