Usenet 1994 October/usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso/misc/volume15/calc-1.05/part08

Lisp/Scheme | 1990-10-14 | 56KB | 1,854 lines

Newsgroups: comp.sources.misc
X-UNIX-From: [e-mail]
subject: v15i035: Patch for GNU Emacs Calc, version 1.04 -> 1.05, part 08/20
from: [e-mail] (David Gillespie)
Sender: [e-mail] (Brandon S. Allbery - comp.sources.misc)

Posting-number: Volume 15, Issue 35
Submitted-by: [e-mail] (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