daveg@csvax.cs.caltech.edu (David Gillespie) (10/15/90)
Posting-number: Volume 15, Issue 35 Submitted-by: daveg@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