v15i035: Patch for GNU Emacs Calc, version 1.04 -> 1.05, part 08/20
David Gillespie
daveg at csvax.cs.caltech.edu
Mon Oct 15 11:17:32 AEST 1990
Posting-number: Volume 15, Issue 35
Submitted-by: daveg at csvax.cs.caltech.edu (David Gillespie)
Archive-name: calc-1.05/part08
#!/bin/sh
# this is part 8 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc.patch 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
sed 's/^X//' << 'SHAR_EOF' >> calc.patch
X! (while (<= (setq n (1+ n)) 0)
X! (setq vec (cons start vec)
X! start (math-mul start (or incr 2)))))
X! (setq vec (nreverse vec)))
X! (if (>= n 0)
X! (while (> n 0)
X! (setq vec (cons n vec)
X! n (1- n)))
X! (let ((i -1))
X! (while (>= i n)
X! (setq vec (cons i vec)
X! i (1- i))))))
X! (cons 'vec vec)))
X )
X (fset 'calcFunc-index (symbol-function 'math-vec-index))
X
X+ ;;; Find an element in a vector.
X+ (defun math-vec-find (vec x &optional start)
X+ (setq start (if start (math-check-fixnum start) 1))
X+ (if (< start 1) (math-reject-arg start 'posp))
X+ (setq vec (nthcdr start vec))
X+ (let ((n start))
X+ (while (and vec (not (math-equal x (car vec))))
X+ (setq n (1+ n)
X+ vec (cdr vec)))
X+ (if vec n 0))
X+ )
X+ (fset 'calcFunc-find (symbol-function 'math-vec-find))
X+
X+ ;;; Return a subvector of a vector.
X+ (defun math-subvector (vec start &optional end)
X+ (setq start (math-check-fixnum start)
X+ end (math-check-fixnum (or end 0)))
X+ (or (math-vectorp vec) (math-reject-arg vec 'vectorp))
X+ (let ((len (1- (length vec))))
X+ (if (<= start 0)
X+ (setq start (+ len start 1)))
X+ (if (<= end 0)
X+ (setq end (+ len end 1)))
X+ (if (or (> start len)
X+ (<= end start))
X+ '(vec)
X+ (setq vec (nthcdr start vec))
X+ (if (<= end len)
X+ (let ((chop (nthcdr (- end start 1) (setq vec (copy-sequence vec)))))
X+ (setcdr chop nil)))
X+ (cons 'vec vec)))
X+ )
X+ (fset 'calcFunc-subvec (symbol-function 'math-subvector))
X+
X+ ;;; Reverse the order of the elements of a vector.
X+ (defun math-reverse-vector (vec)
X+ (if (math-vectorp vec)
X+ (cons 'vec (reverse (cdr vec)))
X+ (math-reject-arg vec 'vectorp))
X+ )
X+ (fset 'calcFunc-rev (symbol-function 'math-reverse-vector))
X+
X+ ;;; Compress a vector according to a mask vector.
X+ (defun math-vec-compress (mask vec)
X+ (if (math-numberp mask)
X+ (if (math-zerop mask)
X+ '(vec)
X+ vec)
X+ (or (math-vectorp mask) (math-reject-arg mask 'vectorp))
X+ (or (math-constp mask) (math-reject-arg mask 'constp))
X+ (or (math-vectorp vec) (math-reject-arg vec 'vectorp))
X+ (or (= (length mask) (length vec)) (math-dimension-error))
X+ (let ((new nil))
X+ (while (setq mask (cdr mask) vec (cdr vec))
X+ (or (math-zerop (car mask))
X+ (setq new (cons (car vec) new))))
X+ (cons 'vec (nreverse new))))
X+ )
X+ (fset 'calcFunc-vmask (symbol-function 'math-vec-compress))
X+
X+ ;;; Expand a vector according to a mask vector.
X+ (defun math-vec-expand (mask vec &optional filler)
X+ (or (math-vectorp mask) (math-reject-arg mask 'vectorp))
X+ (or (math-constp mask) (math-reject-arg mask 'constp))
X+ (or (math-vectorp vec) (math-reject-arg vec 'vectorp))
X+ (setq vec (cdr vec))
X+ (let ((new nil))
X+ (while (setq mask (cdr mask))
X+ (if (math-zerop (car mask))
X+ (setq new (cons (or filler (car mask)) new))
X+ (setq new (cons (or (car vec) (car mask)) new)
X+ vec (cdr vec))))
X+ (cons 'vec (nreverse new)))
X+ )
X+ (fset 'calcFunc-vexp (symbol-function 'math-vec-expand))
X+
X
X ;;; Compute the row and column norms of a vector or matrix. [Public]
X (defun math-rnorm (a)
X***************
X*** 6824,6832 ****
X )
X (fset 'calcFunc-rsort (symbol-function 'math-rsort-vector))
X
X
X ;;; Compile a histogram of data from a vector.
X! (defun math-histogram (vec wts n)
X (or (Math-vectorp vec)
X (math-reject-arg vec 'vectorp))
X (if (Math-vectorp wts)
X--- 11121,11151 ----
X )
X (fset 'calcFunc-rsort (symbol-function 'math-rsort-vector))
X
X+ (defun math-grade-up-vector (grade-vec)
X+ (if (math-vectorp grade-vec)
X+ (let* ((len (1- (length grade-vec))))
X+ (cons 'vec (sort (cdr (math-vec-index len)) 'math-grade-beforep)))
X+ (math-reject-arg grade-vec 'vectorp))
X+ )
X+ (fset 'calcFunc-grade (symbol-function 'math-grade-up-vector))
X+
X+ (defun math-grade-down-vector (grade-vec)
X+ (if (math-vectorp grade-vec)
X+ (let* ((len (1- (length grade-vec))))
X+ (cons 'vec (nreverse (sort (cdr (math-vec-index len))
X+ 'math-grade-beforep))))
X+ (math-reject-arg grade-vec 'vectorp))
X+ )
X+ (fset 'calcFunc-rgrade (symbol-function 'math-grade-down-vector))
X+
X+ (defun math-grade-beforep (i j)
X+ (math-beforep (nth i grade-vec) (nth j grade-vec))
X+ )
X+
X
X ;;; Compile a histogram of data from a vector.
X! (defun math-histogram (vec wts &optional n)
X! (or n (setq n wts wts 1))
X (or (Math-vectorp vec)
X (math-reject-arg vec 'vectorp))
X (if (Math-vectorp wts)
X***************
X*** 7064,7069 ****
X--- 11383,11390 ----
X m
X )
X
X+ ;;;; [calc-arith.el]
X+
X (defun math-abs-approx (a)
X (cond ((Math-negp a)
X (math-neg a))
X***************
X*** 7078,7089 ****
X ((eq (car a) 'intv)
X (math-max (math-abs (nth 2 a)) (math-abs (nth 3 a))))
X ((eq (car a) 'vec)
X! (math-cnorm a))
X ((eq (car a) 'calcFunc-abs)
X (car a))
X (t a))
X )
X
X (defun math-lud-solve (lud b)
X (and lud
X (let* ((x (math-copy-matrix b))
X--- 11399,11416 ----
X ((eq (car a) 'intv)
X (math-max (math-abs (nth 2 a)) (math-abs (nth 3 a))))
X ((eq (car a) 'vec)
X! (math-reduce-vec 'math-add-abs-approx a))
X ((eq (car a) 'calcFunc-abs)
X (car a))
X (t a))
X )
X
X+ (defun math-add-abs-approx (a b)
X+ (math-add (math-abs-approx a) (math-abs-approx b))
X+ )
X+
X+ ;;;; [calc-mat.el]
X+
X (defun math-lud-solve (lud b)
X (and lud
X (let* ((x (math-copy-matrix b))
X***************
X*** 7386,7392 ****
X
X ;;;; Interval forms.
X
X! ;;; Build an interval form. [X I X X]
X (defun math-make-intv (mask lo hi)
X (if (memq (car-safe lo) '(cplx polar mod sdev intv vec))
X (math-reject-arg lo 'realp))
X--- 11713,11719 ----
X
X ;;;; Interval forms.
X
X! ;;; Build an interval form. [X S X X]
X (defun math-make-intv (mask lo hi)
X (if (memq (car-safe lo) '(cplx polar mod sdev intv vec))
X (math-reject-arg lo 'realp))
X***************
X*** 7405,7410 ****
X--- 11732,11743 ----
X (list 'intv mask lo hi))))
X (list 'intv mask lo hi))
X )
X+ (defun calcFunc-intv (mask lo hi)
X+ (if (math-messy-integerp mask) (setq mask (math-trunc mask)))
X+ (or (and (integerp mask) (>= mask 0) (<= mask 3))
X+ (math-reject-arg mask 'range))
X+ (math-make-intv mask lo hi)
X+ )
X
X (defun math-sort-intv (mask lo hi)
X (if (Math-lessp hi lo)
X***************
X*** 7768,7784 ****
X (defun math-want-polar (a b)
X (cond ((eq (car-safe a) 'polar)
X (if (eq (car-safe b) 'cplx)
X! (eq car-complex-mode 'polar)
X t))
X ((eq (car-safe a) 'cplx)
X (if (eq (car-safe b) 'polar)
X! (eq car-complex-mode 'polar)
X nil))
X ((eq (car-safe b) 'polar)
X t)
X ((eq (car-safe b) 'cplx)
X nil)
X! (t (eq (car-complex-mode 'polar))))
X )
X
X ;;; Force A to be in the (-pi,pi] or (-180,180] range.
X--- 12101,12117 ----
X (defun math-want-polar (a b)
X (cond ((eq (car-safe a) 'polar)
X (if (eq (car-safe b) 'cplx)
X! (eq calc-complex-mode 'polar)
X t))
X ((eq (car-safe a) 'cplx)
X (if (eq (car-safe b) 'polar)
X! (eq calc-complex-mode 'polar)
X nil))
X ((eq (car-safe b) 'polar)
X t)
X ((eq (car-safe b) 'cplx)
X nil)
X! (t (eq calc-complex-mode 'polar)))
X )
X
X ;;; Force A to be in the (-pi,pi] or (-180,180] range.
X***************
X*** 8019,8031 ****
X
X ;;;; [calc-arith.el]
X
X (defun math-pow-fancy (a b)
X (cond ((and (Math-numberp a) (Math-numberp b))
X! (cond ((and (eq (car-safe b) 'frac)
X! (equal (nth 2 b) 2))
X! (math-ipow (math-sqrt-raw (math-float a)) (nth 1 b)))
X! ((equal b '(float 5 -1))
X! (math-sqrt-raw (math-float a)))
X (t
X (math-with-extra-prec 2
X (math-exp-raw
X--- 12352,12379 ----
X
X ;;;; [calc-arith.el]
X
X+ (defun math-mod-fancy (a b)
X+ (cond ((and (eq (car-safe a) 'mod) (Math-realp b) (math-posp b))
X+ (math-make-mod (nth 1 a) b))
X+ ((and (eq (car-safe a) 'intv) (math-constp a) (math-posp b))
X+ (math-mod-intv a b))
X+ (t
X+ (if (Math-anglep a)
X+ (calc-record-why 'anglep b)
X+ (calc-record-why 'anglep a))
X+ (list '% a b)))
X+ )
X+
X (defun math-pow-fancy (a b)
X (cond ((and (Math-numberp a) (Math-numberp b))
X! (cond ((memq (math-quarter-integer b) '(1 2 3))
X! (math-pow (math-sqrt (if (math-floatp b) (math-float a) a))
X! (math-mul 2 b)))
X! ((and (eq (car b) 'frac)
X! (integerp (nth 2 b))
X! (<= (nth 2 b) 10))
X! (math-ipow (math-nth-root a (nth 2 b))
X! (nth 1 b)))
X (t
X (math-with-extra-prec 2
X (math-exp-raw
X***************
X*** 8141,8146 ****
X--- 12489,12523 ----
X (math-reject-arg b 'numberp)))
X )
X
X+ (defun math-quarter-integer (x)
X+ (if (Math-integerp x)
X+ 0
X+ (if (math-negp x)
X+ (progn
X+ (setq x (math-quarter-integer (math-neg x)))
X+ (and x (- 4 x)))
X+ (if (eq (car x) 'frac)
X+ (if (eq (nth 2 x) 2)
X+ 2
X+ (and (eq (nth 2 x) 4)
X+ (progn
X+ (setq x (nth 1 x))
X+ (% (if (consp x) (nth 1 x) x) 4))))
X+ (if (eq (car x) 'float)
X+ (if (>= (nth 2 x) 0)
X+ 0
X+ (if (= (nth 2 x) -1)
X+ (progn
X+ (setq x (nth 1 x))
X+ (and (= (% (if (consp x) (nth 1 x) x) 10) 5) 2))
X+ (if (= (nth 2 x) -2)
X+ (progn
X+ (setq x (nth 1 x)
X+ x (% (if (consp x) (nth 1 x) x) 100))
X+ (if (= x 25) 1
X+ (if (= x 75) 3))))))))))
X+ )
X+
X ;;; This assumes A < M and M > 0.
X (defun math-pow-mod (a b m) ; [R R R R]
X (if (and (Math-integerp a) (Math-integerp b) (Math-integerp m))
X***************
X*** 8223,8237 ****
X ;;; with an overestimate always works, even using truncating integer division!
X (defun math-isqrt (a)
X (cond ((Math-zerop a) a)
X! ((Math-negp a)
X! (math-imaginary (math-isqrt (math-neg a))))
X ((integerp a)
X (math-isqrt-small a))
X- ((eq (car a) 'bigpos)
X- (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a))))))
X (t
X! (math-floor (math-sqrt a))))
X )
X
X ;;; This returns (flag . result) where the flag is T if A is a perfect square.
X (defun math-isqrt-bignum (a) ; [P.l L]
X--- 12600,12613 ----
X ;;; with an overestimate always works, even using truncating integer division!
X (defun math-isqrt (a)
X (cond ((Math-zerop a) a)
X! ((not (math-num-natnump a))
X! (math-reject-arg a 'natnump))
X ((integerp a)
X (math-isqrt-small a))
X (t
X! (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a)))))))
X )
X+ (fset 'calcFunc-isqrt (symbol-function 'math-isqrt))
X
X ;;; This returns (flag . result) where the flag is T if A is a perfect square.
X (defun math-isqrt-bignum (a) ; [P.l L]
X***************
X*** 8267,8272 ****
X--- 12643,12655 ----
X guess)))
X )
X
X+ (defun math-zerop-bignum (a)
X+ (and (eq (car a) 0)
X+ (progn
X+ (while (eq (car (setq a (cdr a))) 0))
X+ (null a)))
X+ )
X+
X (defun math-scale-bignum-3 (a n) ; [L L S]
X (while (> n 0)
X (setq a (cons 0 a)
X***************
X*** 8285,8290 ****
X--- 12668,12674 ----
X )
X
X
X+
X ;;;; [calc-ext.el]
X
X (defun math-inexact-result ()
X***************
X*** 8395,8413 ****
X
X ;;; True if A and B differ only in the last digit of precision. [P F F]
X (defun math-nearly-equal-float (a b)
X! (let ((diff (nth 1 (math-sub-float a b))))
X! (or (eq diff 0)
X! (and (not (consp diff))
X! (< diff 10)
X! (> diff -10)
X! (= diff (if (< (nth 2 a) (nth 2 b))
X! (nth 2 a) (nth 2 b)))
X! (or (= (math-numdigs (nth 1 a)) calc-internal-prec)
X! (= (math-numdigs (nth 1 b)) calc-internal-prec)))))
X! )
X!
X! (defun math-nearly-equal (a b) ; [P R R] [Public]
X! (math-nearly-equal-float (math-float a) (math-float b))
X )
X
X ;;; True if A is nearly zero compared to B. [P F F]
X--- 12779,12823 ----
X
X ;;; True if A and B differ only in the last digit of precision. [P F F]
X (defun math-nearly-equal-float (a b)
X! (let ((ediff (- (nth 2 a) (nth 2 b))))
X! (cond ((= ediff 0) ;; Expanded out for speed
X! (setq ediff (math-add (Math-integer-neg (nth 1 a)) (nth 1 b)))
X! (or (eq ediff 0)
X! (and (not (consp ediff))
X! (< ediff 10)
X! (> ediff -10)
X! (= (math-numdigs (nth 1 a)) calc-internal-prec))))
X! ((= ediff 1)
X! (setq ediff (math-add (Math-integer-neg (nth 1 b))
X! (math-scale-int (nth 1 a) 1)))
X! (and (not (consp ediff))
X! (< ediff 10)
X! (> ediff -10)
X! (= (math-numdigs (nth 1 b)) calc-internal-prec)))
X! ((= ediff -1)
X! (setq ediff (math-add (Math-integer-neg (nth 1 a))
X! (math-scale-int (nth 1 b) 1)))
X! (and (not (consp ediff))
X! (< ediff 10)
X! (> ediff -10)
X! (= (math-numdigs (nth 1 a)) calc-internal-prec)))))
X! )
X!
X! (defun math-nearly-equal (a b) ; [P N N] [Public]
X! (setq a (math-float a))
X! (setq b (math-float b))
X! (if (eq (car a) 'polar) (setq a (math-complex a)))
X! (if (eq (car b) 'polar) (setq b (math-complex b)))
X! (if (eq (car a) 'cplx)
X! (if (eq (car b) 'cplx)
X! (and (math-nearly-equal-float (nth 1 a) (nth 1 b))
X! (math-nearly-equal-float (nth 2 a) (nth 2 b)))
X! (and (math-nearly-equal-float (nth 1 a) b)
X! (math-nearly-zerop-float (nth 2 a) b)))
X! (if (eq (car b) 'cplx)
X! (and (math-nearly-equal-float a (nth 1 b))
X! (math-nearly-zerop-float a (nth 2 b)))
X! (math-nearly-equal-float a b)))
X )
X
X ;;; True if A is nearly zero compared to B. [P F F]
X***************
X*** 8417,8424 ****
X (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec))))
X )
X
X! (defun math-nearly-zerop (a b)
X! (math-nearly-zerop-float (math-float a) (math-float b))
X )
X
X ;;; This implementation could be improved, accuracy-wise.
X--- 12827,12841 ----
X (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec))))
X )
X
X! (defun math-nearly-zerop (a b) ; [P N R] [Public]
X! (setq a (math-float a))
X! (setq b (math-float b))
X! (if (eq (car a) 'cplx)
X! (and (math-nearly-zerop-float (nth 1 a) b)
X! (math-nearly-zerop-float (nth 2 a) b))
X! (if (eq (car a) 'polar)
X! (math-nearly-zerop-float (nth 1 a) b)
X! (math-nearly-zerop-float a b)))
X )
X
X ;;; This implementation could be improved, accuracy-wise.
X***************
X*** 8451,8456 ****
X--- 12868,12953 ----
X
X
X
X+ (defun math-nth-root (a n)
X+ (or
X+ (and (= n 2) (math-sqrt a))
X+ (and (Math-zerop a) a)
X+ (and (Math-negp a)
X+ (math-pow a (math-div-float '(float 1 0) (math-float n))))
X+ (and (Math-integerp a)
X+ (let ((root (math-nth-root-integer a n)))
X+ (if (car root)
X+ (cdr root)
X+ (math-nth-root-float (math-float a) n (math-float (cdr root))))))
X+ (and (eq (car-safe a) 'frac)
X+ (let* ((num-root (math-nth-root-integer (nth 1 a) n))
X+ (den-root (math-nth-root-integer (nth 2 a) n)))
X+ (if (and (car num-root) (car den-root))
X+ (list 'frac (cdr num-root) (cdr den-root))
X+ (math-nth-root-float (math-float a) n
X+ (math-div (math-float (cdr num-root))
X+ (math-float (cdr den-root)))))))
X+ (and (eq (car-safe a) 'float)
X+ (math-nth-root-float a n))
X+ (and (eq (car-safe a) 'polar)
X+ (list 'polar
X+ (math-nth-root (nth 1 a) n)
X+ (math-div (nth 2 a) n)))
X+ (math-pow a (math-div-float '(float 1 0) (math-float n))))
X+ )
X+
X+ (defun math-nth-root-float (a n &optional guess)
X+ (math-inexact-result)
X+ (math-with-extra-prec 1
X+ (let ((nf (math-float n))
X+ (nfm1 (math-float (1- n))))
X+ (math-nth-root-float-iter a (or guess
X+ (math-make-float
X+ 1 (/ (+ (math-numdigs (nth 1 a))
X+ (nth 2 a)
X+ (/ n 2))
X+ n))))))
X+ )
X+
X+ (defun math-nth-root-float-iter (a guess) ; uses "n", "nf", "nfm1"
X+ (math-working "root" guess)
X+ (let ((g2 (math-div-float (math-add-float (math-mul nfm1 guess)
X+ (math-div-float
X+ a (math-ipow guess (1- n))))
X+ nf)))
X+ (if (math-nearly-equal-float g2 guess)
X+ g2
X+ (math-nth-root-float-iter a g2)))
X+ )
X+
X+ (defun math-nth-root-integer (a n &optional guess) ; [I I S]
X+ (math-nth-root-int-iter a (or guess
X+ (math-scale-int 1 (/ (+ (math-numdigs a)
X+ (1- n))
X+ n))))
X+ )
X+
X+ (defun math-nth-root-int-iter (a guess) ; uses "n"
X+ (math-working "root" guess)
X+ (let* ((q (math-idivmod a (math-ipow guess (1- n))))
X+ (s (math-add (car q) (math-mul (1- n) guess)))
X+ (g2 (math-idivmod s n)))
X+ (if (Math-natnum-lessp (car g2) guess)
X+ (math-nth-root-int-iter a (car g2))
X+ (cons (and (equal (car g2) guess)
X+ (eq (cdr q) 0)
X+ (eq (cdr g2) 0))
X+ guess)))
X+ )
X+
X+ (defun calcFunc-nroot (x n)
X+ (calcFunc-pow x (if (integerp n)
X+ (math-make-frac 1 n)
X+ (math-div 1 n)))
X+ )
X+
X+
X+
X ;;;; [calc-arith.el]
X
X ;;; Compute the minimum of two real numbers. [R R R] [Public]
X***************
X*** 8565,8571 ****
X
X
X (defun math-trunc-fancy (a)
X! (cond ((eq (car a) 'cplx) (math-trunc (nth 1 a)))
X ((eq (car a) 'polar) (math-trunc (math-complex a)))
X ((eq (car a) 'hms) (list 'hms (nth 1 a) 0 0))
X ((eq (car a) 'mod)
X--- 13062,13069 ----
X
X
X (defun math-trunc-fancy (a)
X! (cond ((eq (car a) 'frac) (math-quotient (nth 1 a) (nth 2 a)))
X! ((eq (car a) 'cplx) (math-trunc (nth 1 a)))
X ((eq (car a) 'polar) (math-trunc (math-complex a)))
X ((eq (car a) 'hms) (list 'hms (nth 1 a) 0 0))
X ((eq (car a) 'mod)
X***************
X*** 8736,8741 ****
X--- 13234,13265 ----
X (fset 'calcFunc-scf (symbol-function 'math-scale-float))
X
X
X+ (defun math-increment (x &optional step relative-to)
X+ (or step (setq step 1))
X+ (cond ((not (Math-integerp step))
X+ (math-reject-arg step 'integerp))
X+ ((Math-integerp x)
X+ (math-add x step))
X+ ((eq (car x) 'float)
X+ (if (and (math-zerop x)
X+ (eq (car-safe relative-to) 'float))
X+ (math-mul step
X+ (math-scale-float relative-to (- 1 calc-internal-prec)))
X+ (math-add-float x (math-make-float
X+ step
X+ (+ (nth 2 x)
X+ (- (math-numdigs (nth 1 x))
X+ calc-internal-prec))))))
X+ (t
X+ (math-reject-arg x 'realp)))
X+ )
X+ (fset 'calcFunc-incr (symbol-function 'math-increment))
X+
X+ (defun calcFunc-decr (x &optional step relative-to)
X+ (math-increment x (math-neg (or step 1)) relative-to)
X+ )
X+
X+
X ;;;; [calc-frac.el]
X
X ;;; Convert a real value to fractional form. [T R I; T R F] [Public]
X***************
X*** 8814,8820 ****
X )
X
X
X! ;;;; [calc-ext.el]
X
X (defun math-clean-number (a &optional prec) ; [X X S] [Public]
X (if prec
X--- 13338,13344 ----
X )
X
X
X! ;;;; [calc-stuff.el]
X
X (defun math-clean-number (a &optional prec) ; [X X S] [Public]
X (if prec
X***************
X*** 8856,8862 ****
X (if (= res 0)
X 1
X (if (= res 2)
X! (list 'calcFunc-eq a b)
X 0)))
X )
X
X--- 13380,13388 ----
X (if (= res 0)
X 1
X (if (= res 2)
X! (if (and (math-looks-negp a) (math-looks-negp b))
X! (list 'calcFunc-eq (math-neg a) (math-neg b))
X! (list 'calcFunc-eq a b))
X 0)))
X )
X
X***************
X*** 8865,8871 ****
X (if (= res 0)
X 0
X (if (= res 2)
X! (list 'calcFunc-neq a b)
X 1)))
X )
X
X--- 13391,13399 ----
X (if (= res 0)
X 0
X (if (= res 2)
X! (if (and (math-looks-negp a) (math-looks-negp b))
X! (list 'calcFunc-neq (math-neg a) (math-neg b))
X! (list 'calcFunc-neq a b))
X 1)))
X )
X
X***************
X*** 8874,8880 ****
X (if (= res -1)
X 1
X (if (= res 2)
X! (list 'calcFunc-lt a b)
X 0)))
X )
X
X--- 13402,13410 ----
X (if (= res -1)
X 1
X (if (= res 2)
X! (if (and (math-looks-negp a) (math-looks-negp b))
X! (list 'calcFunc-gt (math-neg a) (math-neg b))
X! (list 'calcFunc-lt a b))
X 0)))
X )
X
X***************
X*** 8883,8889 ****
X (if (= res 1)
X 1
X (if (= res 2)
X! (list 'calcFunc-gt a b)
X 0)))
X )
X
X--- 13413,13421 ----
X (if (= res 1)
X 1
X (if (= res 2)
X! (if (and (math-looks-negp a) (math-looks-negp b))
X! (list 'calcFunc-lt (math-neg a) (math-neg b))
X! (list 'calcFunc-gt a b))
X 0)))
X )
X
X***************
X*** 8892,8898 ****
X (if (= res 1)
X 0
X (if (= res 2)
X! (list 'calcFunc-leq a b)
X 1)))
X )
X
X--- 13424,13432 ----
X (if (= res 1)
X 0
X (if (= res 2)
X! (if (and (math-looks-negp a) (math-looks-negp b))
X! (list 'calcFunc-geq (math-neg a) (math-neg b))
X! (list 'calcFunc-leq a b))
X 1)))
X )
X
X***************
X*** 8901,8907 ****
X (if (= res -1)
X 0
X (if (= res 2)
X! (list 'calcFunc-geq a b)
X 1)))
X )
X
X--- 13435,13443 ----
X (if (= res -1)
X 0
X (if (= res 2)
X! (if (and (math-looks-negp a) (math-looks-negp b))
X! (list 'calcFunc-leq (math-neg a) (math-neg b))
X! (list 'calcFunc-geq a b))
X 1)))
X )
X
X***************
X*** 8934,8940 ****
X 1
X (if (Math-numberp a)
X 0
X! (list 'calcFunc-lnot a)))
X )
X
X (defun calcFunc-if (c e1 e2)
X--- 13470,13480 ----
X 1
X (if (Math-numberp a)
X 0
X! (let ((op (and (= (length a) 3)
X! (assq (car a) calc-tweak-eqn-table))))
X! (if op
X! (cons (nth 2 op) (cdr a))
X! (list 'calcFunc-lnot a)))))
X )
X
X (defun calcFunc-if (c e1 e2)
X***************
X*** 8942,8948 ****
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--- 13482,13510 ----
X e2
X (if (Math-numberp c)
X e1
X! (or (and (Math-vectorp c)
X! (math-constp c)
X! (let ((ee1 (if (Math-vectorp e1)
X! (if (= (length c) (length e1))
X! (cdr e1)
X! (calc-record-why "Dimension error" e1))
X! (list e1)))
X! (ee2 (if (Math-vectorp e2)
X! (if (= (length c) (length e2))
X! (cdr e2)
X! (calc-record-why "Dimension error" e2))
X! (list e2))))
X! (and ee1 ee2
X! (cons 'vec (math-if-vector (cdr c) ee1 ee2)))))
X! (list 'calcFunc-if c e1 e2))))
X! )
X!
X! (defun math-if-vector (c e1 e2)
X! (and c
X! (cons (if (Math-zerop (car c)) (car e2) (car e1))
X! (math-if-vector (cdr c)
X! (or (cdr e1) e1)
X! (or (cdr e2) e2))))
X )
X
X (defun math-normalize-logical-op (a)
X***************
X*** 8953,8959 ****
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--- 13515,13529 ----
X (math-normalize (nth 3 a))
X (if (Math-numberp a1)
X (math-normalize (nth 2 a))
X! (if (and (Math-vectorp (nth 1 a))
X! (math-constp (nth 1 a)))
X! (calcFunc-if (nth 1 a)
X! (math-normalize (nth 2 a))
X! (math-normalize (nth 3 a)))
X! (let ((calc-simplify-mode 'none))
X! (list 'calcFunc-if a1
X! (math-normalize (nth 2 a))
X! (math-normalize (nth 3 a)))))))))
X a)
X )
X
X***************
X*** 8999,9014 ****
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--- 13569,13575 ----
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 (math-calcFunc-to-var func)))
X )
X
X (defun calcFunc-integer (a)
X***************
X*** 9043,9048 ****
X--- 13604,13637 ----
X 0))
X )
X
X+ (defun calcFunc-negative (a)
X+ (if (math-looks-negp a)
X+ 1
X+ (if (or (math-zerop a)
X+ (math-posp a))
X+ 0
X+ (list 'calcFunc-negative a)))
X+ )
X+
X+ (defun calcFunc-variable (a)
X+ (if (eq (car-safe a) 'var)
X+ 1
X+ (if (Math-objvecp a)
X+ 0
X+ (list 'calcFunc-variable a)))
X+ )
X+
X+ (defun calcFunc-nonvar (a)
X+ (if (eq (car-safe a) 'var)
X+ (list 'calcFunc-nonvar a)
X+ 1)
X+ )
X+
X+ (defun calcFunc-istrue (a)
X+ (if (math-is-true a)
X+ 1
X+ 0)
X+ )
X
X
X
X***************
X*** 9320,9337 ****
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--- 13909,13926 ----
X (fset 'calcFunc-tan (symbol-function 'math-tan))
X
X (defun math-sin-raw (x) ; [N N]
X! (cond ((eq (car 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-add-float expx expmx)
X '(float 5 -1)))
X (math-mul-float (cdr sc)
X! (math-mul-float (math-sub-float expx expmx)
X '(float 5 -1))))))
X! ((eq (car 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***************
X*** 9343,9349 ****
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--- 13932,13938 ----
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 (math-pi-over-2) x)))
X )
X
X ;;; This could use a smarter method: Reduce x as in math-sin-raw, then
X***************
X*** 9354,9361 ****
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--- 13943,13950 ----
X )
X
X (defun math-tan-raw (x) ; [N N]
X! (cond ((eq (car x) 'cplx)
X! (let* ((x (math-mul 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***************
X*** 9365,9373 ****
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--- 13954,13963 ----
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-sub-float expx
X! expmx)
X '(float 5 -1)) d)))))
X! ((eq (car x) 'polar)
X (math-polar (math-tan-raw (math-complex x))))
X (t
X (let ((sc (math-sin-cos-raw x)))
X***************
X*** 9476,9483 ****
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--- 14066,14073 ----
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 x) '(cplx polar))
X! (memq (car 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***************
X*** 9489,9495 ****
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--- 14079,14085 ----
X )
X
X (defun math-arctan-raw (x) ; [N N]
X! (cond ((memq (car x) '(cplx polar))
X (math-with-extra-prec 2 ; extra-extra
X (math-mul '(cplx 0 -1)
X (math-ln-raw (math-mul
X***************
X*** 9643,9653 ****
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--- 14233,14243 ----
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*** 9677,9683 ****
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--- 14267,14288 ----
X (fset 'calcFunc-ln (symbol-function 'math-ln))
X
X (defun math-log10 (x) ; [N N] [Public]
X! (cond ((math-equal-int x 1)
X! (if (math-floatp x) '(float 0 0) 0))
X! ((and (Math-integerp x)
X! (math-posp x)
X! (let ((res (math-integer-log x 10)))
X! (and (car res)
X! (setq x (cdr res)))))
X! x)
X! ((and (eq (car-safe x) 'frac)
X! (eq (nth 1 x) 1)
X! (let ((res (math-integer-log (nth 2 x) 10)))
X! (and (car res)
X! (setq x (- (cdr res))))))
X! x)
X! (calc-symbolic-mode (signal 'inexact-result nil))
X! ((Math-numberp x)
X (math-with-extra-prec 2
X (let ((xf (math-float x)))
X (if (eq (nth 1 xf) 0)
X***************
X*** 9718,9723 ****
X--- 14323,14352 ----
X (math-reject-arg x "Logarithm of zero"))
X ((math-zerop b)
X (math-reject-arg b "Logarithm of zero"))
X+ ((math-equal-int b 1)
X+ (math-reject-arg b "Logarithm base one"))
X+ ((math-equal-int x 1)
X+ (if (or (math-floatp a) (math-floatp b)) '(float 0 0) 0))
X+ ((and (Math-ratp x) (Math-ratp b)
X+ (math-posp x) (math-posp b)
X+ (let* ((sign 1) (inv nil)
X+ (xx (if (math-lessp 1 x)
X+ x
X+ (setq sign -1)
X+ (math-div 1 x)))
X+ (bb (if (math-lessp 1 b)
X+ b
X+ (setq sign (- sign))
X+ (math-div 1 b)))
X+ (res (if (math-lessp xx bb)
X+ (setq inv (math-integer-log bb xx))
X+ (math-integer-log xx bb))))
X+ (and (car res)
X+ (setq x (if inv
X+ (math-div 1 (* sign (cdr res)))
X+ (* sign (cdr res)))))))
X+ x)
X+ (calc-symbolic-mode (signal 'inexact-result nil))
X ((and (Math-numberp x) (Math-numberp b))
X (math-with-extra-prec 2
X (math-div (math-ln-raw (math-float x))
X***************
X*** 9749,9755 ****
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--- 14378,14384 ----
X (math-log x b)))
X (math-normalize (list 'calcFunc-ln x)))
X )
X! (defun calcFunc-alog (x &optional b)
X (if b
X (if (equal b '(var e var-e))
X (math-normalize (list 'calcFunc-exp x))
X***************
X*** 9757,9762 ****
X--- 14386,14425 ----
X (math-normalize (list 'calcFunc-exp x)))
X )
X
X+ (defun calcFunc-ilog (x b)
X+ (or (math-num-integerp x) (math-reject-arg x 'integerp))
X+ (or (math-num-integerp b) (math-reject-arg b 'integerp))
X+ (setq x (math-trunc x))
X+ (setq b (math-trunc b))
X+ (or (math-posp x) (math-reject-arg x 'posp))
X+ (or (math-posp b) (math-reject-arg b 'posp))
X+ (if (eq b 1) (math-reject-arg x "Logarithm of zero"))
X+ (if (Math-natnum-lessp x b)
X+ 0
X+ (cdr (math-integer-log x b)))
X+ )
X+
X+ (defun math-integer-log (x b)
X+ (let ((pows (list b))
X+ (pow (math-sqr b))
X+ next
X+ sum n)
X+ (while (not (math-lessp x pow))
X+ (setq pows (cons pow pows)
X+ pow (math-sqr pow)))
X+ (setq n (lsh 1 (1- (length pows)))
X+ sum n
X+ pow (car pows))
X+ (while (and (setq pows (cdr pows))
X+ (math-lessp pow x))
X+ (setq n (/ n 2)
X+ next (math-mul pow (car pows)))
X+ (or (math-lessp x next)
X+ (setq pow next
X+ sum (+ sum n))))
X+ (cons (equal pow x) sum))
X+ )
X+
X
X (defun math-log-base-raw (b) ; [N N]
X (if (not (and (equal (car math-log-base-cache) b)
X***************
X*** 10036,10041 ****
X--- 14699,15504 ----
X
X
X
X+ ;;;; [calc-funcs.el]
X+
X+ ;;; Sources: Numerical Recipes, Press et al;
X+ ;;; Handbook of Mathematical Functions, Abramowitz & Stegun.
X+
X+
X+ ;;; Gamma function.
X+
X+ (defun calcFunc-gamma (x)
X+ (or (math-numberp x) (math-reject-arg x 'numberp))
X+ (calcFunc-fact (math-add x -1))
X+ )
X+
X+ (defun math-gammap1-raw (x fprec nfprec) ; compute gamma(1 + x)
X+ (cond ((math-lessp-float (math-real-part x) fprec)
X+ (if (math-lessp-float (math-real-part x) nfprec)
X+ (math-neg (math-div
X+ (math-pi)
X+ (math-mul (math-gammap1-raw
X+ (math-add (math-neg x)
X+ '(float -1 0))
X+ fprec nfprec)
X+ (math-sin-raw
X+ (math-mul (math-pi) x)))))
X+ (let ((xplus1 (math-add x '(float 1 0))))
X+ (math-div (math-gammap1-raw xplus1 fprec nfprec) xplus1))))
X+ (t ; re(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+ '(float 0 0)
X+ 2))))))
X+ )
X+
X+ (defun math-gamma-series (sum x xinvsqr oterm n)
X+ (math-working "gamma" sum)
X+ (let* ((bn (math-bernoulli-number n))
X+ (term (math-mul (math-div-float (math-float (nth 1 bn))
X+ (math-float (* (nth 2 bn)
X+ (* n (1- n)))))
X+ x))
X+ (next (math-add sum term)))
X+ (if (math-nearly-equal sum next)
X+ next
X+ (if (> n (* 2 calc-internal-prec))
X+ (progn
X+ ;; Need this because series eventually diverges for large enough n.
X+ (calc-record-why "Gamma computation stopped early, not all digits may be valid")
X+ next)
X+ (math-gamma-series next (math-mul x xinvsqr) xinvsqr term (+ n 2)))))
X+ )
X+
X+
X+ ;;; Incomplete gamma function.
X+
X+ (defun calcFunc-gammaP (a x)
X+ (math-inexact-result)
X+ (or (Math-numberp a) (math-reject-arg a 'numberp))
X+ (or (math-numberp x) (math-reject-arg x 'numberp))
X+ (if (and (math-num-integerp a)
X+ (integerp (setq a (math-trunc a)))
X+ (> a 0) (< a 20))
X+ (math-sub 1 (calcFunc-gammaQ a x))
X+ (let ((math-current-gamma-value (calcFunc-gamma a)))
X+ (math-div (calcFunc-gammag a x) math-current-gamma-value)))
X+ )
X+
X+ (defun calcFunc-gammaQ (a x)
X+ (math-inexact-result)
X+ (or (Math-numberp a) (math-reject-arg a 'numberp))
X+ (or (math-numberp x) (math-reject-arg x 'numberp))
X+ (if (and (math-num-integerp a)
X+ (integerp (setq a (math-trunc a)))
X+ (> a 0) (< a 20))
X+ (let ((n 0)
X+ (sum '(float 1 0))
X+ (term '(float 1 0)))
X+ (math-with-extra-prec 1
X+ (while (< (setq n (1+ n)) a)
X+ (setq term (math-div (math-mul term x) n)
X+ sum (math-add sum term))
X+ (math-working "gamma" sum))
X+ (math-mul sum (math-exp (math-neg x)))))
X+ (let ((math-current-gamma-value (calcFunc-gamma a)))
X+ (math-div (calcFunc-gammaG a x) math-current-gamma-value)))
X+ )
X+
X+ (defun calcFunc-gammag (a x)
X+ (math-inexact-result)
X+ (or (Math-numberp a) (math-reject-arg a 'numberp))
X+ (or (Math-numberp x) (math-reject-arg x 'numberp))
X+ (math-with-extra-prec 2
X+ (setq a (math-float a))
X+ (setq x (math-float x))
X+ (if (or (math-negp (math-real-part a))
X+ (math-lessp-float (math-real-part x)
X+ (math-add-float (math-real-part a)
X+ '(float 1 0))))
X+ (math-inc-gamma-series a x)
X+ (math-sub (or math-current-gamma-value (calcFunc-gamma a))
X+ (math-inc-gamma-cfrac a x))))
X+ )
X+ (setq math-current-gamma-value nil)
X+
X+ (defun calcFunc-gammaG (a x)
X+ (math-inexact-result)
X+ (or (Math-numberp a) (math-reject-arg a 'numberp))
X+ (or (Math-numberp x) (math-reject-arg x 'numberp))
X+ (math-with-extra-prec 2
X+ (setq a (math-float a))
X+ (setq x (math-float x))
X+ (if (or (math-negp (math-real-part a))
X+ (math-lessp-float (math-real-part x)
X+ (math-add-float (math-abs-approx a)
X+ '(float 1 0))))
X+ (math-sub (or math-current-gamma-value (calcFunc-gamma a))
X+ (math-inc-gamma-series a x))
X+ (math-inc-gamma-cfrac a x)))
X+ )
X+
X+ (defun math-inc-gamma-series (a x)
X+ (math-mul (math-exp-raw (math-sub (math-mul a (math-ln-raw x)) x))
X+ (math-with-extra-prec 2
X+ (let ((start (math-div '(float 1 0) a)))
X+ (math-inc-gamma-series-step start start a x))))
X+ )
X+
X+ (defun math-inc-gamma-series-step (sum term a x)
X+ (math-working "gamma" sum)
X+ (setq a (math-add a '(float 1 0))
X+ term (math-div (math-mul term x) a))
X+ (let ((next (math-add sum term)))
X+ (if (math-nearly-equal sum next)
X+ next
X+ (math-inc-gamma-series-step next term a x)))
X+ )
X+
X+ (defun math-inc-gamma-cfrac (a x)
X+ (math-mul (math-exp-raw (math-sub (math-mul a (math-ln-raw x)) x))
X+ (math-inc-gamma-cfrac-step '(float 1 0) x '(float 0 0) '(float 1 0)
X+ '(float 1 0) '(float 1 0) '(float 0 0)
X+ a x))
X+ )
X+
X+ (defun math-inc-gamma-cfrac-step (a0 a1 b0 b1 n fac g a x)
X+ (let ((ana (math-sub n a))
X+ (anf (math-mul n fac)))
X+ (setq n (math-add n '(float 1 0))
X+ a0 (math-mul (math-add a1 (math-mul a0 ana)) fac)
X+ b0 (math-mul (math-add b1 (math-mul b0 ana)) fac)
X+ a1 (math-add (math-mul x a0) (math-mul anf a1))
X+ b1 (math-add (math-mul x b0) (math-mul anf b1)))
X+ (if (math-zerop a1)
X+ (math-inc-gamma-cfrac-step a0 a1 b0 b1 n fac g a x)
X+ (setq fac (math-div '(float 1 0) a1))
X+ (let ((next (math-mul b1 fac)))
X+ (math-working "gamma" next)
X+ (if (math-nearly-equal next g)
X+ next
X+ (math-inc-gamma-cfrac-step a0 a1 b0 b1 n fac next a x)))))
X+ )
X+
X+
X+ ;;; Error function.
X+
X+ (defun calcFunc-erf (x)
X+ (let ((math-current-gamma-value (math-sqrt-pi)))
X+ (math-to-same-complex-quad
X+ (math-div (calcFunc-gammag '(float 5 -1)
X+ (math-sqr (math-to-complex-quad-one x)))
X+ math-current-gamma-value)
X+ x))
X+ )
X+
X+ (defun calcFunc-erfc (x)
X+ (if (math-posp x)
X+ (let ((math-current-gamma-value (math-sqrt-pi)))
X+ (math-div (calcFunc-gammaG '(float 5 -1) (math-sqr x))
X+ math-current-gamma-value))
X+ (math-add '(float 1 0) (calcFunc-erf (math-neg x))))
X+ )
X+
X+ (defun math-to-complex-quad-one (x)
X+ (if (eq (car-safe x) 'polar) (setq x (math-complex x)))
X+ (if (eq (car-safe x) 'cplx)
X+ (list 'cplx (math-abs (nth 1 x)) (math-abs (nth 2 x)))
X+ x)
X+ )
X+
X+ (defun math-to-same-complex-quad (x y)
X+ (if (eq (car-safe y) 'cplx)
X+ (if (eq (car-safe x) 'cplx)
X+ (list 'cplx
X+ (if (math-negp (nth 1 y)) (math-neg (nth 1 x)) (nth 1 x))
X+ (if (math-negp (nth 2 y)) (math-neg (nth 2 x)) (nth 2 x)))
X+ (if (math-negp (nth 1 y)) (math-neg x) x))
X+ (if (math-negp y)
X+ (if (eq (car-safe x) 'cplx)
X+ (list 'cplx (math-neg (nth 1 x)) (nth 2 x))
X+ (math-neg x))
X+ x))
X+ )
X+
X+
X+ ;;; Beta function.
X+
X+ (defun calcFunc-beta (a b)
X+ (if (math-num-integerp a)
X+ (let ((am (math-add a -1)))
X+ (or (math-numberp b) (math-reject-arg b 'numberp))
X+ (math-div 1 (math-mul b (math-choose (math-add b am) am))))
X+ (if (math-num-integerp b)
X+ (calcFunc-beta b a)
X+ (math-div (math-mul (calcFunc-gamma a) (calcFunc-gamma b))
X+ (calcFunc-gamma (math-add a b)))))
X+ )
X+
X+
X+ ;;; Incomplete beta function.
X+
X+ (defun calcFunc-betaI (x a b)
X+ (cond ((math-zerop x)
X+ '(float 0 0))
X+ ((math-equal-int x 1)
X+ '(float 1 0))
X+ ((or (math-zerop a)
X+ (and (math-num-integerp a)
X+ (math-negp a)))
X+ (if (or (math-zerop b)
X+ (and (math-num-integerp b)
X+ (math-negp b)))
X+ (math-reject-arg b 'nonzerop)
X+ '(float 1 0)))
X+ ((or (math-zerop b)
X+ (and (math-num-integerp b)
X+ (math-negp b)))
X+ '(float 0 0))
X+ ((not (math-numberp a)) (math-reject-arg a 'numberp))
X+ ((not (math-numberp b)) (math-reject-arg b 'numberp))
X+ ((math-inexact-result))
X+ (t (let ((math-current-beta-value (calcFunc-beta a b)))
X+ (math-div (calcFunc-betaB x a b) math-current-beta-value))))
X+ )
X+
X+ (defun calcFunc-betaB (x a b)
X+ (cond
X+ ((math-zerop x)
X+ '(float 0 0))
X+ ((math-equal-int x 1)
X+ (calcFunc-beta a b))
X+ ((not (math-numberp x)) (math-reject-arg x 'numberp))
X+ ((not (math-numberp a)) (math-reject-arg a 'numberp))
X+ ((not (math-numberp b)) (math-reject-arg b 'numberp))
X+ ((math-zerop a) (math-reject-arg a 'nonzerop))
X+ ((math-zerop b) (math-reject-arg b 'nonzerop))
X+ ((and (math-num-integerp b)
X+ (if (math-negp b)
X+ (math-reject-arg b 'range)
X+ (Math-natnum-lessp (setq b (math-trunc b)) 20)))
X+ (and calc-symbolic-mode (or (math-floatp a) (math-floatp b))
X+ (math-inexact-result))
X+ (math-mul
X+ (math-with-extra-prec 2
X+ (let* ((i 0)
X+ (term 1)
X+ (sum (math-div term a)))
X+ (while (< (setq i (1+ i)) b)
X+ (setq term (math-mul (math-div (math-mul term (- i b)) i) x)
X+ sum (math-add sum (math-div term (math-add a i))))
X+ (math-working "beta" sum))
X+ sum))
X+ (math-pow x a)))
X+ ((and (math-num-integerp a)
X+ (if (math-negp a)
X+ (math-reject-arg a 'range)
X+ (Math-natnum-lessp (setq a (math-trunc a)) 20)))
X+ (math-sub (or math-current-beta-value (calcFunc-beta a b))
X+ (calcFunc-betaB (math-sub 1 x) b a)))
X+ (t
X+ (math-inexact-result)
X+ (math-with-extra-prec 2
X+ (setq x (math-float x))
X+ (setq a (math-float a))
X+ (setq b (math-float b))
X+ (let ((bt (math-exp-raw (math-add (math-mul a (math-ln-raw x))
X+ (math-mul b (math-ln-raw
X+ (math-sub '(float 1 0)
X+ x)))))))
X+ (if (math-lessp x (math-div (math-add a '(float 1 0))
X+ (math-add (math-add a b) '(float 2 0))))
X+ (math-div (math-mul bt (math-beta-cfrac a b x)) a)
X+ (math-sub (or math-current-beta-value (calcFunc-beta a b))
X+ (math-div (math-mul bt
X+ (math-beta-cfrac b a (math-sub 1 x)))
X+ b)))))))
X+ )
X+ (setq math-current-beta-value nil)
X+
X+ (defun math-beta-cfrac (a b x)
X+ (let ((qab (math-add a b))
X+ (qap (math-add a '(float 1 0)))
X+ (qam (math-add a '(float -1 0))))
X+ (math-beta-cfrac-step '(float 1 0)
X+ (math-sub '(float 1 0)
X+ (math-div (math-mul qab x) qap))
X+ '(float 1 0) '(float 1 0)
X+ '(float 1 0)
X+ qab qap qam a b x))
X+ )
X+
X+ (defun math-beta-cfrac-step (az bz am bm m qab qap qam a b x)
X+ (let* ((two-m (math-mul m '(float 2 0)))
X+ (d (math-div (math-mul (math-mul (math-sub b m) m) x)
X+ (math-mul (math-add qam two-m) (math-add a two-m))))
X+ (ap (math-add az (math-mul d am)))
X+ (bp (math-add bz (math-mul d bm)))
X+ (d2 (math-neg
X+ (math-div (math-mul (math-mul (math-add a m) (math-add qab m)) x)
X+ (math-mul (math-add qap two-m) (math-add a two-m)))))
X+ (app (math-add ap (math-mul d2 az)))
X+ (bpp (math-add bp (math-mul d2 bz)))
X+ (next (math-div app bpp)))
X+ (math-working "beta" next)
X+ (if (math-nearly-equal next az)
X+ next
X+ (math-beta-cfrac-step next '(float 1 0)
X+ (math-div ap bpp) (math-div bp bpp)
X+ (math-add m '(float 1 0))
X+ qab qap qam a b x)))
X+ )
X+
X+
X+ ;;; Bessel functions.
X+
X+ ;;; Should generalize this to handle arbitrary precision!
X+
X+ (defun calcFunc-besJ (v x)
X+ (or (math-numberp v) (math-reject-arg v 'numberp))
X+ (or (math-numberp x) (math-reject-arg x 'numberp))
X+ (let ((calc-internal-prec (min 8 calc-internal-prec)))
X+ (math-with-extra-prec 3
X+ (setq x (math-float (math-normalize x)))
X+ (setq v (math-float (math-normalize v)))
X+ (cond ((math-zerop x)
X+ (if (math-zerop v)
X+ '(float 1 0)
X+ '(float 0 0)))
X+ ((math-inexact-result))
X+ ((not (math-num-integerp v))
X+ (let ((start (math-div 1 (calcFunc-fact v))))
X+ (math-mul (math-besJ-series start start
X+ 0
X+ (math-mul '(float -25 -2)
X+ (math-sqr x))
X+ v)
X+ (math-pow (math-div x 2) v))))
X+ ((math-negp (setq v (math-trunc v)))
X+ (if (math-oddp v)
X+ (math-neg (calcFunc-besJ (math-neg v) x))
X+ (calcFunc-besJ (math-neg v) x)))
X+ ((eq v 0)
X+ (math-besJ0 x))
X+ ((eq v 1)
X+ (math-besJ1 x))
X+ ((math-lessp v (math-abs-approx x))
X+ (let ((j 0)
X+ (bjm (math-besJ0 x))
X+ (bj (math-besJ1 x))
X+ (two-over-x (math-div 2 x))
X+ bjp)
X+ (while (< (setq j (1+ j)) v)
X+ (setq bjp (math-sub (math-mul (math-mul j two-over-x) bj)
X+ bjm)
X+ bjm bj
X+ bj bjp))
X+ bj))
X+ (t
X+ (if (math-lessp 100 v) (math-reject-arg v 'range))
X+ (let* ((j (logior (+ v (math-isqrt-small (* 40 v))) 1))
X+ (two-over-x (math-div 2 x))
X+ (jsum nil)
X+ (bjp '(float 0 0))
X+ (sum '(float 0 0))
X+ (bj '(float 1 0))
X+ bjm ans)
X+ (while (> (setq j (1- j)) 0)
X+ (setq bjm (math-sub (math-mul (math-mul j two-over-x) bj)
X+ bjp)
X+ bjp bj
X+ bj bjm)
X+ (if (> (nth 2 (math-abs-approx bj)) 10)
X+ (setq bj (math-mul bj '(float 1 -10))
X+ bjp (math-mul bjp '(float 1 -10))
X+ ans (and ans (math-mul ans '(float 1 -10)))
X+ sum (math-mul sum '(float 1 -10))))
X+ (or (setq jsum (not jsum))
X+ (setq sum (math-add sum bj)))
X+ (if (= j v)
X+ (setq ans bjp)))
X+ (math-div ans (math-sub (math-mul 2 sum) bj)))))))
X+ )
X+
X+ (defun math-besJ-series (sum term k zz vk)
X+ (math-working "besJ" sum)
X+ (setq k (1+ k)
X+ vk (math-add 1 vk)
X+ term (math-div (math-mul term zz) (math-mul k vk)))
X+ (let ((next (math-add sum term)))
X+ (if (math-nearly-equal next sum)
X+ next
X+ (math-besJ-series next term k zz vk)))
X+ )
X+
X+ (defun math-besJ0 (x &optional yflag)
X+ (cond ((and (not yflag) (math-negp (math-real-part x)))
X+ (math-besJ0 (math-neg x)))
X+ ((math-lessp '(float 8 0) (math-abs-approx x))
X+ (let* ((z (math-div '(float 8 0) x))
X+ (y (math-sqr z))
X+ (xx (math-add x '(float (bigneg 164 398 785) -9)))
X+ (a1 (math-poly-eval y
X+ '((float (bigpos 211 887 093 2) -16)
X+ (float (bigneg 639 370 073 2) -15)
X+ (float (bigpos 407 510 734 2) -14)
X+ (float (bigneg 627 628 098 1) -12)
X+ (float 1 0))))
X+ (a2 (math-poly-eval y
X+ '((float (bigneg 152 935 934) -16)
X+ (float (bigpos 161 095 621 7) -16)
X+ (float (bigneg 651 147 911 6) -15)
X+ (float (bigpos 765 488 430 1) -13)
X+ (float (bigneg 995 499 562 1) -11))))
X+ (sc (math-sin-cos-raw xx)))
X+ (if yflag
X+ (setq sc (cons (math-neg (cdr sc)) (car sc))))
X+ (math-mul (math-sqrt
X+ (math-div '(float (bigpos 722 619 636) -9) x))
X+ (math-sub (math-mul (cdr sc) a1)
X+ (math-mul (car sc) (math-mul z a2))))))
X+ (t
X+ (let ((y (math-sqr x)))
X+ (math-div (math-poly-eval y
X+ '((float (bigneg 456 052 849 1) -7)
X+ (float (bigpos 017 233 739 7) -5)
X+ (float (bigneg 418 442 121 1) -2)
X+ (float (bigpos 407 196 516 6) -1)
X+ (float (bigneg 354 590 362 13) 0)
X+ (float (bigpos 574 490 568 57) 0)))
X+ (math-poly-eval y
X+ '((float 1 0)
X+ (float (bigpos 712 532 678 2) -7)
X+ (float (bigpos 853 264 927 5) -5)
X+ (float (bigpos 718 680 494 9) -3)
X+ (float (bigpos 985 532 029 1) 0)
X+ (float (bigpos 411 490 568 57) 0)))))))
X+ )
X+
X+ (defun math-besJ1 (x &optional yflag)
X+ (cond ((and (math-negp (math-real-part x)) (not yflag))
X+ (math-neg (math-besJ1 (math-neg x))))
X+ ((math-lessp '(float 8 0) (math-abs-approx x))
X+ (let* ((z (math-div '(float 8 0) x))
X+ (y (math-sqr z))
X+ (xx (math-add x '(float (bigneg 491 194 356 2) -9)))
X+ (a1 (math-poly-eval y
X+ '((float (bigneg 019 337 240) -15)
X+ (float (bigpos 174 520 457 2) -15)
X+ (float (bigneg 496 396 516 3) -14)
X+ (float 183105 -8)
X+ (float 1 0))))
X+ (a2 (math-poly-eval y
X+ '((float (bigpos 412 787 105) -15)
X+ (float (bigneg 987 228 88) -14)
X+ (float (bigpos 096 199 449 8) -15)
X+ (float (bigneg 873 690 002 2) -13)
X+ (float (bigpos 995 499 687 4) -11))))
X+ (sc (math-sin-cos-raw xx)))
X+ (if yflag
X+ (setq sc (cons (math-neg (cdr sc)) (car sc)))
X+ (if (math-negp x)
X+ (setq sc (cons (math-neg (car sc)) (math-neg (cdr sc))))))
X+ (math-mul (math-sqrt (math-div '(float (bigpos 722 619 636) -9) x))
X+ (math-sub (math-mul (cdr sc) a1)
X+ (math-mul (car sc) (math-mul z a2))))))
X+ (t
X+ (let ((y (math-sqr x)))
X+ (math-mul
X+ x
X+ (math-div (math-poly-eval y
X+ '((float (bigneg 606 036 016 3) -8)
X+ (float (bigpos 826 044 157) -4)
X+ (float (bigneg 439 611 972 2) -3)
X+ (float (bigpos 531 968 423 2) -1)
X+ (float (bigneg 235 059 895 7) 0)
X+ (float (bigpos 232 614 362 72) 0)))
X+ (math-poly-eval y
X+ '((float 1 0)
X+ (float (bigpos 397 991 769 3) -7)
X+ (float (bigpos 394 743 944 9) -5)
X+ (float (bigpos 474 330 858 1) -2)
X+ (float (bigpos 178 535 300 2) 0)
X+ (float (bigpos 442 228 725 144)
X+ 0))))))))
X+ )
X+
X+ (defun calcFunc-besY (v x)
X+ (math-inexact-result)
X+ (or (math-numberp v) (math-reject-arg v 'numberp))
X+ (or (math-numberp x) (math-reject-arg x 'numberp))
X+ (let ((calc-internal-prec (min 8 calc-internal-prec)))
X+ (math-with-extra-prec 3
X+ (setq x (math-float (math-normalize x)))
X+ (setq v (math-float (math-normalize v)))
X+ (cond ((not (math-num-integerp v))
X+ (let ((sc (math-sin-cos-raw (math-mul v (math-pi)))))
X+ (math-div (math-sub (math-mul (calcFunc-besJ v x) (cdr sc))
X+ (calcFunc-besJ (math-neg v) x))
X+ (car sc))))
X+ ((math-negp (setq v (math-trunc v)))
X+ (if (math-oddp v)
X+ (math-neg (calcFunc-besY (math-neg v) x))
X+ (calcFunc-besY (math-neg v) x)))
X+ ((eq v 0)
X+ (math-besY0 x))
X+ ((eq v 1)
X+ (math-besY1 x))
X+ (t
X+ (let ((j 0)
X+ (bym (math-besY0 x))
X+ (by (math-besY1 x))
X+ (two-over-x (math-div 2 x))
X+ byp)
X+ (while (< (setq j (1+ j)) v)
X+ (setq byp (math-sub (math-mul (math-mul j two-over-x) by)
X+ bym)
X+ bym by
X+ by byp))
X+ by)))))
X+ )
X+
X+ (defun math-besY0 (x)
X+ (cond ((math-lessp (math-abs-approx x) '(float 8 0))
X+ (let ((y (math-sqr x)))
X+ (math-add
X+ (math-div (math-poly-eval y
X+ '((float (bigpos 733 622 284 2) -7)
X+ (float (bigneg 757 792 632 8) -5)
X+ (float (bigpos 129 988 087 1) -2)
X+ (float (bigneg 036 598 123 5) -1)
X+ (float (bigpos 065 834 062 7) 0)
X+ (float (bigneg 389 821 957 2) 0)))
X+ (math-poly-eval y
X+ '((float 1 0)
X+ (float (bigpos 244 030 261 2) -7)
X+ (float (bigpos 647 472 474) -4)
X+ (float (bigpos 438 466 189 7) -3)
X+ (float (bigpos 648 499 452 7) -1)
X+ (float (bigpos 269 544 076 40) 0))))
X+ (math-mul '(float (bigpos 772 619 636) -9)
X+ (math-mul (math-besJ0 x) (math-ln-raw x))))))
X+ ((math-negp (math-real-part x))
X+ (math-add (math-besJ0 (math-neg x) t)
X+ (math-mul '(cplx 0 2)
X+ (math-besJ0 (math-neg x)))))
X+ (t
X+ (math-besJ0 x t)))
X+ )
X+
X+ (defun math-besY1 (x)
X+ (cond ((math-lessp (math-abs-approx x) '(float 8 0))
X+ (let ((y (math-sqr x)))
X+ (math-add
X+ (math-mul
X+ x
X+ (math-div (math-poly-eval y
X+ '((float (bigpos 935 937 511 8) -6)
X+ (float (bigneg 726 922 237 4) -3)
X+ (float (bigpos 551 264 349 7) -1)
X+ (float (bigneg 139 438 153 5) 1)
X+ (float (bigpos 439 527 127) 4)
X+ (float (bigneg 943 604 900 4) 3)))
X+ (math-poly-eval y
X+ '((float 1 0)
X+ (float (bigpos 885 632 549 3) -7)
X+ (float (bigpos 605 042 102) -3)
X+ (float (bigpos 002 904 245 2) -2)
X+ (float (bigpos 367 650 733 3) 0)
X+ (float (bigpos 664 419 244 4) 2)
X+ (float (bigpos 057 958 249) 5)))))
X+ (math-mul '(float (bigpos 772 619 636) -9)
X+ (math-sub (math-mul (math-besJ1 x) (math-ln-raw x))
X+ (math-div 1 x))))))
X+ ((math-negp (math-real-part x))
X+ (math-neg
X+ (math-add (math-besJ1 (math-neg x) t)
X+ (math-mul '(cplx 0 2)
X+ (math-besJ1 (math-neg x))))))
X+ (t
X+ (math-besJ1 x t)))
X+ )
X+
X+ (defun math-poly-eval (x coefs)
X+ (let ((accum (car coefs)))
X+ (while (setq coefs (cdr coefs))
X+ (setq accum (math-add (car coefs) (math-mul accum x))))
X+ accum)
X+ )
X+
X+
X+ ;;;; Bernoulli and Euler polynomials and numbers.
X+
X+ (defun calcFunc-bern (n &optional x)
X+ (if (and x (not (math-zerop x)))
X+ (if (and calc-symbolic-mode (math-floatp x))
X+ (math-inexact-result)
X+ (math-build-polynomial-expr (math-bernoulli-coefs n) x))
X+ (or (math-num-natnump n) (math-reject-arg n 'natnump))
X+ (if (consp n)
X+ (progn
X+ (math-inexact-result)
X+ (math-float (math-bernoulli-number (math-trunc n))))
SHAR_EOF
echo "End of part 8, continue with part 9"
echo "9" > s2_seq_.tmp
exit 0
More information about the Comp.sources.misc
mailing list