Usenet 1994 January

home *** CD-ROM | disk | FTP | other *** search

/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume13 / gmcalc / part07 < prev next >

Wrap

Text File | 1990-06-05 | 57.1 KB | 1,916 lines

Newsgroups: comp.sources.misc From: daveg@csvax.caltech.edu (David Gillespie) Subject: v13i033: Emacs Calculator 1.01, part 07/19 Sender: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc) Posting-number: Volume 13, Issue 33 Submitted-by: daveg@csvax.caltech.edu (David Gillespie) Archive-name: gmcalc/part07 ---- Cut Here and unpack ---- #!/bin/sh # this is part 7 of a multipart archive # do not concatenate these parts, unpack them in order with /bin/sh # file calc-ext.el continued # CurArch=7 if test ! -r s2_seq_.tmp then echo "Please unpack part 1 first!" exit 1; fi ( read Scheck if test "$Scheck" != $CurArch then echo "Please unpack part $Scheck next!" exit 1; else exit 0; fi ) < s2_seq_.tmp || exit 1 echo "x - Continuing file calc-ext.el" sed 's/^X//' << 'SHAR_EOF' >> calc-ext.el X X;;; Create a vector of consecutive integers. [Public] X(defun math-vec-index (n) X (and (not (integerp n)) X (setq n (math-check-fixnum n))) X (or (natnump n) (math-reject-arg n 'natnump)) X (let ((vec nil)) X (while (> n 0) X (setq vec (cons n vec) X n (1- n))) X (cons 'vec vec)) X) X(fset 'calcFunc-index (symbol-function 'math-vec-index)) X X X;;; Compute the row and column norms of a vector or matrix. [Public] X(defun math-rnorm (a) X (if (and (Math-vectorp a) X (math-constp a)) X (if (math-matrixp a) X (math-reduce-vec 'math-max (math-map-vec 'math-cnorm a)) X (math-reduce-vec 'math-max (math-map-vec 'math-abs a))) X (calc-record-why 'vectorp a) X (list 'calcFunc-rnorm a)) X) X(fset 'calcFunc-rnorm (symbol-function 'math-rnorm)) X X(defun math-cnorm (a) X (if (and (Math-vectorp a) X (math-constp a)) X (if (math-matrixp a) X (math-reduce-vec 'math-max X (math-reduce-cols 'math-add-abs a)) X (math-reduce-vec 'math-add-abs a)) X (calc-record-why 'vectorp a) X (list 'calcFunc-cnorm a)) X) X(fset 'calcFunc-cnorm (symbol-function 'math-cnorm)) X X(defun math-add-abs (a b) X (math-add (math-abs a) (math-abs b)) X) X X X;;; Sort the elements of a vector into increasing order. X(defun math-sort-vector (vec) ; [Public] X (if (math-vectorp vec) X (cons 'vec (sort (copy-sequence (cdr vec)) 'math-beforep)) X (math-reject-arg vec 'vectorp)) X) X(fset 'calcFunc-sort (symbol-function 'math-sort-vector)) X X(defun math-rsort-vector (vec) ; [Public] X (if (math-vectorp vec) X (cons 'vec (nreverse (sort (copy-sequence (cdr vec)) 'math-beforep))) X (math-reject-arg vec 'vectorp)) 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 (or (= (length vec) (length wts)) X (math-dimension-error))) X (or (natnump n) X (math-reject-arg n 'natnump)) X (let ((res (make-vector n 0)) X (vp vec) X (wvec (Math-vectorp wts)) X (wp wts) X bin) X (while (setq vp (cdr vp)) X (setq bin (car vp)) X (or (natnump bin) X (setq bin (math-floor bin))) X (and (natnump bin) X (< bin n) X (aset res bin (math-add (aref res bin) X (if wvec (car (setq wp (cdr wp))) wts))))) X (cons 'vec (append res nil))) X) X(fset 'calcFunc-histogram (symbol-function 'math-histogram)) X X X(defun math-matrix-trace (mat) ; [Public] X (if (math-square-matrixp mat) X (math-matrix-trace-step 2 (1- (length mat)) mat (nth 1 (nth 1 mat))) X (math-reject-arg mat 'square-matrixp)) X) X(fset 'calcFunc-tr (symbol-function 'math-matrix-trace)) X X(defun math-matrix-trace-step (n size mat sum) X (if (<= n size) X (math-matrix-trace-step (1+ n) size mat X (math-add sum (nth n (nth n mat)))) X sum) X) X X X;;; Matrix inverse and determinant. X(defun math-matrix-inv-raw (m) X (let ((n (1- (length m)))) X (if (<= n 3) X (let ((det (math-det-raw m))) X (and (not (math-zerop det)) X (math-div X (cond ((= n 1) 1) X ((= n 2) X (list 'vec X (list 'vec X (nth 2 (nth 2 m)) X (math-neg (nth 2 (nth 1 m)))) X (list 'vec X (math-neg (nth 1 (nth 2 m))) X (nth 1 (nth 1 m))))) X ((= n 3) X (list 'vec X (list 'vec X (math-sub (math-mul (nth 3 (nth 3 m)) X (nth 2 (nth 2 m))) X (math-mul (nth 3 (nth 2 m)) X (nth 2 (nth 3 m)))) X (math-sub (math-mul (nth 3 (nth 1 m)) X (nth 2 (nth 3 m))) X (math-mul (nth 3 (nth 3 m)) X (nth 2 (nth 1 m)))) X (math-sub (math-mul (nth 3 (nth 2 m)) X (nth 2 (nth 1 m))) X (math-mul (nth 3 (nth 1 m)) X (nth 2 (nth 2 m))))) X (list 'vec X (math-sub (math-mul (nth 3 (nth 2 m)) X (nth 1 (nth 3 m))) X (math-mul (nth 3 (nth 3 m)) X (nth 1 (nth 2 m)))) X (math-sub (math-mul (nth 3 (nth 3 m)) X (nth 1 (nth 1 m))) X (math-mul (nth 3 (nth 1 m)) X (nth 1 (nth 3 m)))) X (math-sub (math-mul (nth 3 (nth 1 m)) X (nth 1 (nth 2 m))) X (math-mul (nth 3 (nth 2 m)) X (nth 1 (nth 1 m))))) X (list 'vec X (math-sub (math-mul (nth 2 (nth 3 m)) X (nth 1 (nth 2 m))) X (math-mul (nth 2 (nth 2 m)) X (nth 1 (nth 3 m)))) X (math-sub (math-mul (nth 2 (nth 1 m)) X (nth 1 (nth 3 m))) X (math-mul (nth 2 (nth 3 m)) X (nth 1 (nth 1 m)))) X (math-sub (math-mul (nth 2 (nth 2 m)) X (nth 1 (nth 1 m))) X (math-mul (nth 2 (nth 1 m)) X (nth 1 (nth 2 m)))))))) X det))) X (let ((lud (math-matrix-lud m))) X (and lud X (math-lud-solve lud (math-diag-matrix 1 n)))))) X) X X(defun math-matrix-det (m) X (if (math-square-matrixp m) X (math-with-extra-prec 2 (math-det-raw m)) X (math-reject-arg m 'square-matrixp)) X) X(fset 'calcFunc-det (symbol-function 'math-matrix-det)) X X(defun math-det-raw (m) X (let ((n (1- (length m)))) X (cond ((= n 1) X (nth 1 (nth 1 m))) X ((= n 2) X (math-sub (math-mul (nth 1 (nth 1 m)) X (nth 2 (nth 2 m))) X (math-mul (nth 2 (nth 1 m)) X (nth 1 (nth 2 m))))) X ((= n 3) X (math-sub X (math-sub X (math-sub X (math-add X (math-add X (math-mul (nth 1 (nth 1 m)) X (math-mul (nth 2 (nth 2 m)) X (nth 3 (nth 3 m)))) X (math-mul (nth 2 (nth 1 m)) X (math-mul (nth 3 (nth 2 m)) X (nth 1 (nth 3 m))))) X (math-mul (nth 3 (nth 1 m)) X (math-mul (nth 1 (nth 2 m)) X (nth 2 (nth 3 m))))) X (math-mul (nth 3 (nth 1 m)) X (math-mul (nth 2 (nth 2 m)) X (nth 1 (nth 3 m))))) X (math-mul (nth 1 (nth 1 m)) X (math-mul (nth 3 (nth 2 m)) X (nth 2 (nth 3 m))))) X (math-mul (nth 2 (nth 1 m)) X (math-mul (nth 1 (nth 2 m)) X (nth 3 (nth 3 m)))))) X (t (let ((lud (math-matrix-lud m))) X (if lud X (let ((lu (car lud))) X (math-det-step n (nth 2 lud))) X 0))))) X) X X(defun math-det-step (n prod) X (if (> n 0) X (math-det-step (1- n) (math-mul prod (nth n (nth n lu)))) X prod) X) X X;;; This returns a list (LU index d), or NIL if not possible. X;;; Argument M must be a square matrix. X(defun math-matrix-lud (m) X (let ((old (assoc m math-lud-cache)) X (context (list calc-internal-prec calc-prefer-frac))) X (if (and old (equal (nth 1 old) context)) X (cdr (cdr old)) X (let* ((lud (catch 'singular (math-do-matrix-lud m))) X (entry (cons context lud))) X (if old X (setcdr old entry) X (setq math-lud-cache (cons (cons m entry) math-lud-cache))) X lud))) X) X(defvar math-lud-cache nil) X X;;; Numerical Recipes section 2.3; implicit pivoting omitted. X(defun math-do-matrix-lud (m) X (let* ((lu (math-copy-matrix m)) X (n (1- (length lu))) X i (j 1) k imax sum big X (d 1) (index nil)) X (while (<= j n) X (setq i 1 X big 0 X imax j) X (while (< i j) X (math-working "LUD step" (format "%d/%d" j i)) X (setq sum (nth j (nth i lu)) X k 1) X (while (< k i) X (setq sum (math-sub sum (math-mul (nth k (nth i lu)) X (nth j (nth k lu)))) X k (1+ k))) X (setcar (nthcdr j (nth i lu)) sum) X (setq i (1+ i))) X (while (<= i n) X (math-working "LUD step" (format "%d/%d" j i)) X (setq sum (nth j (nth i lu)) X k 1) X (while (< k j) X (setq sum (math-sub sum (math-mul (nth k (nth i lu)) X (nth j (nth k lu)))) X k (1+ k))) X (setcar (nthcdr j (nth i lu)) sum) X (let ((dum (math-abs-approx sum))) X (if (Math-lessp big dum) X (setq big dum X imax i))) X (setq i (1+ i))) X (if (> imax j) X (setq lu (math-swap-rows lu j imax) X d (- d))) X (setq index (cons imax index)) X (let ((pivot (nth j (nth j lu)))) X (if (math-zerop pivot) X (throw 'singular nil) X (setq i j) X (while (<= (setq i (1+ i)) n) X (setcar (nthcdr j (nth i lu)) X (math-div (nth j (nth i lu)) pivot))))) X (setq j (1+ j))) X (list lu (nreverse index) d)) X) X X(defun math-swap-rows (m r1 r2) X (or (= r1 r2) X (let* ((r1prev (nthcdr (1- r1) m)) X (row1 (cdr r1prev)) X (r2prev (nthcdr (1- r2) m)) X (row2 (cdr r2prev)) X (r2next (cdr row2))) X (setcdr r2prev row1) X (setcdr r1prev row2) X (setcdr row2 (cdr row1)) X (setcdr row1 r2next))) X m X) X X(defun math-abs-approx (a) X (cond ((Math-negp a) X (math-neg a)) X ((Math-anglep a) X a) X ((eq (car a) 'cplx) X (math-add (math-abs (nth 1 a)) (math-abs (nth 2 a)))) X ((eq (car a) 'polar) X (nth 1 a)) X ((eq (car a) 'sdev) X (math-abs (nth 1 a))) 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 (n (1- (length x))) X (m (1- (length (nth 1 x)))) X (lu (car lud)) X (col 1) X i j ip ii index sum) X (while (<= col m) X (math-working "LUD solver step" col) X (setq i 1 X ii nil X index (nth 1 lud)) X (while (<= i n) X (setq ip (car index) X index (cdr index) X sum (nth col (nth ip x))) X (setcar (nthcdr col (nth ip x)) (nth col (nth i x))) X (if (null ii) X (or (math-zerop sum) X (setq ii i)) X (setq j ii) X (while (< j i) X (setq sum (math-sub sum (math-mul (nth j (nth i lu)) X (nth col (nth j x)))) X j (1+ j)))) X (setcar (nthcdr col (nth i x)) sum) X (setq i (1+ i))) X (while (>= (setq i (1- i)) 1) X (setq sum (nth col (nth i x)) X j i) X (while (<= (setq j (1+ j)) n) X (setq sum (math-sub sum (math-mul (nth j (nth i lu)) X (nth col (nth j x)))))) X (setcar (nthcdr col (nth i x)) X (math-div sum (nth i (nth i lu))))) X (setq col (1+ col))) X x)) X) X X(defun calcFunc-lud (m) X (if (math-square-matrixp m) X (or (math-with-extra-prec 2 X (let ((lud (math-matrix-lud m))) X (and lud X (let* ((lmat (math-copy-matrix (car lud))) X (umat (math-copy-matrix (car lud))) X (n (1- (length (car lud)))) X (perm (math-diag-matrix 1 n)) X i (j 1)) X (while (<= j n) X (setq i 1) X (while (< i j) X (setcar (nthcdr j (nth i lmat)) 0) X (setq i (1+ i))) X (setcar (nthcdr j (nth j lmat)) 1) X (while (<= (setq i (1+ i)) n) X (setcar (nthcdr j (nth i umat)) 0)) X (setq j (1+ j))) X (while (>= (setq j (1- j)) 1) X (let ((pos (nth (1- j) (nth 1 lud)))) X (or (= pos j) X (setq perm (math-swap-rows perm j pos))))) X (list 'vec perm lmat umat))))) X (math-reject-arg m "Singular matrix")) X (math-reject-arg m 'square-matrixp)) X) X X;;; Compute a right-handed vector cross product. [O O O] [Public] X(defun math-cross (a b) X (if (and (eq (car-safe a) 'vec) X (= (length a) 4)) X (if (and (eq (car-safe b) 'vec) X (= (length b) 4)) X (list 'vec X (math-sub (math-mul (nth 2 a) (nth 3 b)) X (math-mul (nth 3 a) (nth 2 b))) X (math-sub (math-mul (nth 3 a) (nth 1 b)) X (math-mul (nth 1 a) (nth 3 b))) X (math-sub (math-mul (nth 1 a) (nth 2 b)) X (math-mul (nth 2 a) (nth 1 b)))) X (math-reject-arg b "Three-vector expected")) X (math-reject-arg a "Three-vector expected")) X) X(fset 'calcFunc-cross (symbol-function 'math-cross)) X X X X X;;;; Hours-minutes-seconds forms. X X(defun math-normalize-hms (a) X (let ((h (math-normalize (nth 1 a))) X (m (math-normalize (nth 2 a))) X (s (let ((calc-internal-prec (max (- calc-internal-prec 4) 3))) X (math-normalize (nth 3 a))))) X (if (math-negp h) X (progn X (if (math-posp s) X (setq s (math-add s -60) X m (math-add m 1))) X (if (math-posp m) X (setq m (math-add m -60) X h (math-add h 1))) X (if (not (math-lessp -60 s)) X (setq s (math-add s 60) X m (math-add m -1))) X (if (not (math-lessp -60 m)) X (setq m (math-add m 60) X h (math-add h -1)))) X (if (math-negp s) X (setq s (math-add s 60) X m (math-add m -1))) X (if (math-negp m) X (setq m (math-add m 60) X h (math-add h -1))) X (if (not (math-lessp s 60)) X (setq s (math-add s -60) X m (math-add m 1))) X (if (not (math-lessp m 60)) X (setq m (math-add m -60) X h (math-add h 1)))) X (if (and (eq (car-safe s) 'float) X (<= (+ (math-numdigs (nth 1 s)) (nth 2 s)) X (- 2 calc-internal-prec))) X (setq s 0)) X (list 'hms h m s)) X) X X;;; Convert A from ANG or current angular mode to HMS format. X(defun math-to-hms (a &optional ang) ; [X R] [Public] X (cond ((eq (car-safe a) 'hms) a) X ((eq (car-safe a) 'sdev) X (math-make-sdev (math-to-hms (nth 1 a)) X (math-to-hms (nth 2 a)))) X ((not (Math-numberp a)) X (list 'calcFunc-hms a)) X ((math-negp a) X (math-neg (math-to-hms (math-neg a) ang))) X ((eq (or ang calc-angle-mode) 'rad) X (math-to-hms (math-div a (math-pi-over-180)) 'deg)) X ((memq (car-safe a) '(cplx polar)) a) X (t X ;(setq a (let ((calc-internal-prec (max (1- calc-internal-prec) 3))) X ; (math-normalize a))) X (math-normalize X (let* ((b (math-mul a 3600)) X (hm (math-trunc (math-div b 60))) X (hmd (math-idivmod hm 60))) X (list 'hms X (car hmd) X (cdr hmd) X (math-sub b (math-mul hm 60))))))) X) X(fset 'calcFunc-hms (symbol-function 'math-to-hms)) X X;;; Convert A from HMS format to ANG or current angular mode. X(defun math-from-hms (a &optional ang) ; [R X] [Public] X (cond ((not (eq (car-safe a) 'hms)) X (if (Math-numberp a) X a X (if (eq (car-safe a) 'sdev) X (math-make-sdev (math-from-hms (nth 1 a)) X (math-from-hms (nth 2 a))) X (if (eq (or ang calc-angle-mode) 'rad) X (list '>rad a) X (list '>deg a))))) X ((math-negp a) X (math-neg (math-from-hms (math-neg a) ang))) X ((eq (or ang calc-angle-mode) 'rad) X (math-mul (math-from-hms a 'deg) (math-pi-over-180))) X ((memq (car-safe a) '(cplx polar)) a) X (t X (math-add (math-div (math-add (math-div (nth 3 a) X '(float 6 1)) X (nth 2 a)) X 60) X (nth 1 a)))) X) X X X X;;;; Complex numbers. X X(defun math-normalize-polar (a) X (let ((r (math-normalize (nth 1 a))) X (th (math-normalize (nth 2 a)))) X (cond ((math-zerop r) X '(polar 0 0)) X ((or (math-zerop th)) X r) X ((and (not (eq calc-angle-mode 'rad)) X (or (equal th '(float 18 1)) X (equal th 180))) X (math-neg r)) X ((math-negp r) X (math-neg (list 'polar (math-neg r) th))) X (t X (list 'polar r th)))) X) X X X;;; Coerce A to be complex (rectangular form). [c N] X(defun math-complex (a) X (cond ((eq (car-safe a) 'cplx) a) X ((eq (car-safe a) 'polar) X (if (math-zerop (nth 1 a)) X (nth 1 a) X (let ((sc (math-sin-cos (nth 2 a)))) X (list 'cplx X (math-mul (nth 1 a) (nth 1 sc)) X (math-mul (nth 1 a) (nth 2 sc)))))) X (t (list 'cplx a 0))) X) X X;;; Coerce A to be complex (polar form). [c N] X(defun math-polar (a) X (cond ((eq (car-safe a) 'polar) a) X ((math-zerop a) '(polar 0 0)) X (t X (list 'polar X (math-abs a) X (math-cplx-arg a)))) X) X X;;; Multiply A by the imaginary constant i. [N N] [Public] X(defun math-imaginary (a) X (if (and (Math-objvecp a) X (not calc-symbolic-mode)) X (math-mul a X (if (or (eq (car-safe a) 'polar) X (and (not (eq (car-safe a) 'cplx)) X (eq calc-complex-mode 'polar))) X (list 'polar 1 (math-quarter-circle nil)) X '(cplx 0 1))) X (math-mul a '(var i var-i))) X) X X X;;;; Error forms. X X;;; Build a standard deviation form. [X X X] X(defun math-make-sdev (x sigma) X (if (memq (car-safe x) '(cplx polar mod sdev intv vec)) X (math-reject-arg x 'realp)) X (if (memq (car-safe sigma) '(cplx polar mod sdev intv vec)) X (math-reject-arg sigma 'realp)) X (if (math-negp sigma) X (list 'sdev x (math-neg sigma)) X (if (and (math-zerop sigma) (Math-scalarp x)) X x X (list 'sdev x sigma))) X) X(defun calcFunc-sdev (x sigma) X (math-make-sdev x sigma) X) X X X X;;;; Modulo forms. X X(defun math-normalize-mod (a) X (let ((n (math-normalize (nth 1 a))) X (m (math-normalize (nth 2 a)))) X (if (and (math-anglep n) (math-anglep m) (math-posp m)) X (math-make-mod n m) X (if (math-anglep n) X (if (math-anglep m) X (calc-record-why "Modulus must be positive" m) X (calc-record-why "Modulus must be real" m)) X (calc-record-why "Value must be real" n)) X (list 'calcFunc-makemod n m))) X) X X;;; Build a modulo form. [N R R] X(defun math-make-mod (n m) X (setq calc-previous-modulo m) X (and n X (if (not (and (Math-anglep n) (Math-anglep m))) X (math-reject-arg n 'anglep) X (if (or (Math-negp n) X (not (Math-lessp n m))) X (list 'mod (math-mod n m) m) X (list 'mod n m)))) X) X(defun calcFunc-makemod (n m) X (math-make-mod n m) X) X X 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 (if (memq (car-safe hi) '(cplx polar mod sdev intv vec)) X (math-reject-arg hi 'realp)) X (if (and (Math-realp lo) (Math-realp hi)) X (let ((cmp (math-compare lo hi))) X (if (= cmp 0) X (if (= mask 3) X lo X (list 'intv mask lo hi)) X (if (> cmp 0) X (if (= mask 3) X (list 'intv 2 lo lo) X (list 'intv mask lo lo)) X (list 'intv mask lo hi)))) X (list 'intv mask lo hi)) X) X X(defun math-sort-intv (mask lo hi) X (if (Math-lessp hi lo) X (math-make-intv (aref [0 2 1 3] mask) hi lo) X (math-make-intv mask lo hi)) X) X X X X;;;; Arithmetic. X X(defun math-neg-fancy (a) X (cond ((eq (car a) 'polar) X (list 'polar X (nth 1 a) X (if (math-posp (nth 2 a)) X (math-sub (nth 2 a) (math-half-circle nil)) X (math-add (nth 2 a) (math-half-circle nil))))) X ((eq (car a) 'mod) X (if (math-zerop (nth 1 a)) X a X (list 'mod (math-sub (nth 2 a) (nth 1 a)) (nth 2 a)))) X ((eq (car a) 'sdev) X (list 'sdev (math-neg (nth 1 a)) (nth 2 a))) X ((eq (car a) 'intv) X (math-make-intv (aref [0 2 1 3] (nth 1 a)) X (math-neg (nth 3 a)) X (math-neg (nth 2 a)))) X ((eq (car a) '-) X (math-sub (nth 2 a) (nth 1 a))) X ((and (memq (car a) '(* /)) X (math-looks-negp (nth 1 a))) X (list (car a) (math-neg (nth 1 a)) (nth 2 a))) X ((and (memq (car a) '(* /)) X (math-looks-negp (nth 2 a))) X (list (car a) (nth 1 a) (math-neg (nth 2 a)))) X ((and (memq (car a) '(* /)) X (or (math-numberp (nth 1 a)) X (and (eq (car (nth 1 a)) '*) X (math-numberp (nth 1 (nth 1 a)))))) X (list (car a) (math-neg (nth 1 a)) (nth 2 a))) X ((and (eq (car a) '/) X (or (math-numberp (nth 2 a)) X (and (eq (car (nth 2 a)) '*) X (math-numberp (nth 1 (nth 2 a)))))) X (list (car a) (nth 1 a) (math-neg (nth 2 a)))) X ((eq (car a) 'neg) X (nth 1 a)) X (t (list 'neg a))) X) X X(defun math-neg-float (a) X (list 'float (Math-integer-neg (nth 1 a)) (nth 2 a)) X) X X(defun math-add-objects-fancy (a b) X (cond ((and (Math-numberp a) (Math-numberp b)) X (setq aa (math-complex a) X bb (math-complex b)) X (math-normalize X (let ((res (list 'cplx X (math-add (nth 1 aa) (nth 1 bb)) X (math-add (nth 2 aa) (nth 2 bb))))) X (if (math-want-polar a b) X (math-polar res) X res)))) X ((or (Math-vectorp a) (Math-vectorp b)) X (math-map-vec-2 'math-add a b)) X ((eq (car-safe a) 'sdev) X (if (eq (car-safe b) 'sdev) X (math-make-sdev (math-add (nth 1 a) (nth 1 b)) X (math-hypot (nth 2 a) (nth 2 b))) X (and (or (Math-anglep b) X (not (Math-objvecp b))) X (math-make-sdev (math-add (nth 1 a) b) X (nth 2 a))))) X ((and (eq (car-safe b) 'sdev) X (or (Math-anglep a) X (not (Math-objvecp a)))) X (math-make-sdev (math-add a (nth 1 b)) X (nth 2 b))) X ((eq (car-safe a) 'intv) X (if (eq (car-safe b) 'intv) X (math-make-intv (logand (nth 1 a) (nth 1 b)) X (math-add (nth 2 a) (nth 2 b)) X (math-add (nth 3 a) (nth 3 b))) X (and (or (Math-anglep b) X (not (Math-objvecp b))) X (math-make-intv (nth 1 a) X (math-add (nth 2 a) b) X (math-add (nth 3 a) b))))) X ((and (eq (car-safe b) 'intv) X (or (Math-anglep a) X (not (Math-objvecp a)))) X (math-make-intv (nth 1 b) X (math-add a (nth 2 b)) X (math-add a (nth 3 b)))) X ((and (eq (car-safe a) 'mod) X (eq (car-safe b) 'mod) X (equal (nth 2 a) (nth 2 b))) X (math-make-mod (math-add (nth 1 a) (nth 1 b)) (nth 2 a))) X ((and (eq (car-safe a) 'mod) X (Math-anglep b)) X (math-make-mod (math-add (nth 1 a) b) (nth 2 a))) X ((and (eq (car-safe b) 'mod) X (Math-anglep a)) X (math-make-mod (math-add a (nth 1 b)) (nth 2 b))) X ((and (or (eq (car-safe a) 'hms) (eq (car-safe b) 'hms)) X (and (Math-anglep a) (Math-anglep b))) X (or (eq (car-safe a) 'hms) (setq a (math-to-hms a))) X (or (eq (car-safe b) 'hms) (setq b (math-to-hms b))) X (math-normalize X (if (math-negp a) X (math-neg (math-add (math-neg a) (math-neg b))) X (if (math-negp b) X (let* ((s (math-add (nth 3 a) (nth 3 b))) X (m (math-add (nth 2 a) (nth 2 b))) X (h (math-add (nth 1 a) (nth 1 b)))) X (if (math-negp s) X (setq s (math-add s 60) X m (math-add m -1))) X (if (math-negp m) X (setq m (math-add m 60) X h (math-add h -1))) X (if (math-negp h) X (math-add b a) X (list 'hms h m s))) X (let* ((s (math-add (nth 3 a) (nth 3 b))) X (m (math-add (nth 2 a) (nth 2 b))) X (h (math-add (nth 1 a) (nth 1 b)))) X (list 'hms h m s)))))) X (t (calc-record-why "Incompatible arguments" a b))) X) X X(defun math-add-symb-fancy (a b) X (or (and (eq (car-safe b) '+) X (math-add (math-add a (nth 1 b)) X (nth 2 b))) X (and (eq (car-safe b) '-) X (math-sub (math-add a (nth 1 b)) X (nth 2 b))) X (and (eq (car-safe b) 'neg) X (eq (car-safe (nth 1 b)) '+) X (math-sub (math-sub a (nth 1 (nth 1 b))) X (nth 2 (nth 1 b)))) X (cond X ((eq (car-safe a) '+) X (let ((temp (math-combine-sum (nth 2 a) b nil nil t))) X (and temp X (math-add (nth 1 a) temp)))) X ((eq (car-safe a) '-) X (let ((temp (math-combine-sum (nth 2 a) b t nil t))) X (and temp X (math-add (nth 1 a) temp)))) X ((and (Math-objectp a) (Math-objectp b)) X nil) X (t X (math-combine-sum a b nil nil nil))) X (and (Math-looks-negp b) X (list '- a (math-neg b))) X (and (Math-looks-negp a) X (list '- b (math-neg a))) X (list '+ a b)) X) X X X(defun math-mul-objects-fancy (a b) X (cond ((and (Math-numberp a) (Math-numberp b)) X (math-normalize X (if (math-want-polar a b) X (let ((a (math-polar a)) X (b (math-polar b))) X (list 'polar X (math-mul (nth 1 a) (nth 1 b)) X (math-fix-circular (math-add (nth 2 a) (nth 2 b))))) X (setq a (math-complex a) X b (math-complex b)) X (list 'cplx X (math-sub (math-mul (nth 1 a) (nth 1 b)) X (math-mul (nth 2 a) (nth 2 b))) X (math-add (math-mul (nth 1 a) (nth 2 b)) X (math-mul (nth 2 a) (nth 1 b))))))) X ((Math-vectorp a) X (if (Math-vectorp b) X (if (math-matrixp a) X (if (math-matrixp b) X (cons 'vec (math-mul-mats (cdr a) X (mapcar 'cdr X (cdr b)))) X (math-mat-col X (cons 'vec X (if (= (length (nth 1 a)) 2) X (math-mul-mats (cdr a) X (mapcar 'cdr X (cdr (math-row-matrix X b)))) X (math-mul-mats (cdr a) X (mapcar 'cdr X (cdr (math-col-matrix X b)))))) X 1)) X (if (math-matrixp b) X (cons 'vec (math-mul-mat-row a (mapcar 'cdr (cdr b)))) X (car (math-mul-mat-row a X (mapcar 'cdr X (cdr (math-col-matrix X b))))))) X (math-map-vec-2 'math-mul a b))) X ((Math-vectorp b) X (math-map-vec-2 'math-mul a b)) X ((eq (car-safe a) 'sdev) X (if (eq (car-safe b) 'sdev) X (math-make-sdev (math-mul (nth 1 a) (nth 1 b)) X (math-hypot (math-mul (nth 2 a) (nth 1 b)) X (math-mul (nth 2 b) (nth 1 a)))) X (and (or (Math-anglep b) X (not (Math-objvecp b))) X (math-make-sdev (math-mul (nth 1 a) b) X (math-abs (math-mul (nth 2 a) b)))))) X ((and (eq (car-safe b) 'sdev) X (or (Math-anglep a) X (not (Math-objvecp a)))) X (math-make-sdev (math-mul a (nth 1 b)) X (math-abs (math-mul a (nth 2 b))))) X ((and (eq (car-safe a) 'intv) (Math-anglep b)) X (if (Math-negp b) X (math-neg (math-mul a (math-neg b))) X (math-make-intv (nth 1 a) X (math-mul (nth 2 a) b) X (math-mul (nth 3 a) b)))) X ((and (eq (car-safe b) 'intv) (Math-anglep a)) X (math-mul b a)) X ((and (eq (car-safe a) 'intv) (math-constp a) X (eq (car-safe b) 'intv) (math-constp b)) X (let ((lo (math-mul a (nth 2 b))) X (hi (math-mul a (nth 3 b)))) X (and (Math-anglep lo) X (setq lo (list 'intv (if (memq (nth 1 b) '(2 3)) 3 0) lo lo))) X (and (Math-anglep hi) X (setq hi (list 'intv (if (memq (nth 1 b) '(1 3)) 3 0) hi hi))) X (math-combine-intervals (nth 2 lo) (and (memq (nth 1 b) '(2 3)) X (memq (nth 1 lo) '(2 3))) X (nth 3 lo) (and (memq (nth 1 b) '(2 3)) X (memq (nth 1 lo) '(1 3))) X (nth 2 hi) (and (memq (nth 1 b) '(1 3)) X (memq (nth 1 hi) '(2 3))) X (nth 3 hi) (and (memq (nth 1 b) '(1 3)) X (memq (nth 1 hi) '(1 3)))))) X ((and (eq (car-safe a) 'mod) X (eq (car-safe b) 'mod) X (equal (nth 2 a) (nth 2 b))) X (math-make-mod (math-mul (nth 1 a) (nth 1 b)) (nth 2 a))) X ((and (eq (car-safe a) 'mod) X (Math-anglep b)) X (math-make-mod (math-mul (nth 1 a) b) (nth 2 a))) X ((and (eq (car-safe b) 'mod) X (Math-anglep a)) X (math-make-mod (math-mul a (nth 1 b)) (nth 2 b))) X ((and (eq (car-safe a) 'hms) (Math-realp b)) X (math-with-extra-prec 2 X (math-to-hms (math-mul (math-from-hms a 'deg) b) 'deg))) X ((and (eq (car-safe b) 'hms) (Math-realp a)) X (math-mul b a)) X (t (calc-record-why "Incompatible arguments" a b))) X) X X;;; Fast function to multiply floating-point numbers. X(defun math-mul-float (a b) ; [F F F] X (math-make-float (math-mul (nth 1 a) (nth 1 b)) X (+ (nth 2 a) (nth 2 b))) X) X X(defun math-sqr-float (a) ; [F F] X (math-make-float (math-mul (nth 1 a) (nth 1 a)) X (+ (nth 2 a) (nth 2 a))) X) X X(defun math-combine-intervals (a am b bm c cm d dm) X (let (res) X (if (= (setq res (math-compare a c)) 1) X (setq a c am cm) X (if (= res 0) X (setq am (or am cm)))) X (if (= (setq res (math-compare b d)) -1) X (setq b d bm dm) X (if (= res 0) X (setq bm (or bm dm)))) X (math-make-intv (+ (if am 2 0) (if bm 1 0)) a b)) X) X X(defun math-mul-symb-fancy (a b) X (or (and (Math-equal-int a 1) X b) X (and (Math-equal-int a -1) X (math-neg b)) X (and (Math-numberp b) X (math-mul b a)) X (and (eq (car-safe a) 'neg) X (math-neg (math-mul (nth 1 a) b))) X (and (eq (car-safe b) 'neg) X (math-neg (math-mul a (nth 1 b)))) X (and (eq (car-safe a) '*) X (math-mul (nth 1 a) X (math-mul (nth 2 a) b))) X (and (eq (car-safe a) '^) X (Math-looks-negp (nth 2 a)) X (not (and (eq (car-safe b) '^) (Math-looks-negp (nth 2 b)))) X (math-div b (math-normalize X (list '^ (nth 1 a) (math-neg (nth 2 a)))))) X (and (eq (car-safe b) '^) X (Math-looks-negp (nth 2 b)) X (not (and (eq (car-safe a) '^) (Math-looks-negp (nth 2 a)))) X (math-div a (math-normalize X (list '^ (nth 1 b) (math-neg (nth 2 b)))))) X (and (eq (car-safe a) '/) X (math-div (math-mul (nth 1 a) b) (nth 2 a))) X (and (eq (car-safe b) '/) X (math-div (math-mul a (nth 1 b)) (nth 2 b))) X (and (eq (car-safe b) '+) X (Math-numberp a) X (or (Math-numberp (nth 1 b)) X (Math-numberp (nth 2 b))) X (math-add (math-mul a (nth 1 b)) X (math-mul a (nth 2 b)))) X (and (eq (car-safe b) '-) X (Math-numberp a) X (or (Math-numberp (nth 1 b)) X (Math-numberp (nth 2 b))) X (math-sub (math-mul a (nth 1 b)) X (math-mul a (nth 2 b)))) X (and (eq (car-safe b) '*) X (Math-numberp (nth 1 b)) X (not (Math-numberp a)) X (math-mul (nth 1 b) (math-mul a (nth 2 b)))) X (and (or t ; this seems more reasonable... X (eq (car-safe a) '-) X (math-looks-negp a)) X (math-looks-negp b) X (math-mul (math-neg a) (math-neg b))) X (and (eq (car-safe b) '-) X (math-looks-negp a) X (math-mul (math-neg a) (math-neg b))) X (cond X ((eq (car-safe b) '*) X (let ((temp (math-combine-prod a (nth 1 b) nil nil t))) X (and temp X (math-mul temp (nth 2 b))))) X (t X (math-combine-prod a b nil nil nil))) X (list '* a b)) X) X 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(defun math-fix-circular (a &optional dir) ; [R R] X (cond ((eq calc-angle-mode 'deg) X (cond ((and (Math-lessp '(float 18 1) a) (not (eq dir 1))) X (math-fix-circular (math-add a '(float -36 1)) -1)) X ((or (Math-lessp '(float -18 1) a) (eq dir -1)) X a) X (t X (math-fix-circular (math-add a '(float 36 1)) 1)))) X ((eq calc-angle-mode 'hms) X (cond ((and (Math-lessp 180 (nth 1 a)) (not (eq dir 1))) X (math-fix-circular (math-add a '(float -36 1)) -1)) X ((or (Math-lessp -180 (nth 1 a)) (eq dir -1)) X a) X (t X (math-fix-circular (math-add a '(float 36 1)) 1)))) X (t X (cond ((and (Math-lessp (math-pi) a) (not (eq dir 1))) X (math-fix-circular (math-sub a (math-two-pi)) -1)) X ((or (Math-lessp (math-neg (math-pi)) a) (eq dir -1)) X a) X (t X (math-fix-circular (math-add a (math-two-pi)) 1))))) X) X X X(defun math-div-objects-fancy (a b) X (cond ((and (Math-numberp a) (Math-numberp b)) X (math-normalize X (cond ((math-want-polar a b) X (let ((a (math-polar a)) X (b (math-polar b))) X (list 'polar X (math-div (nth 1 a) (nth 1 b)) X (math-fix-circular (math-sub (nth 2 a) X (nth 2 b)))))) X ((Math-realp b) X (setq a (math-complex a)) X (list 'cplx (math-div (nth 1 a) b) X (math-div (nth 2 a) b))) X (t X (setq a (math-complex a) X b (math-complex b)) X (math-div X (list 'cplx X (math-add (math-mul (nth 1 a) (nth 1 b)) X (math-mul (nth 2 a) (nth 2 b))) X (math-sub (math-mul (nth 2 a) (nth 1 b)) X (math-mul (nth 1 a) (nth 2 b)))) X (math-add (math-sqr (nth 1 b)) X (math-sqr (nth 2 b)))))))) X ((math-matrixp b) X (if (math-square-matrixp b) X (let ((n1 (length b))) X (if (Math-vectorp a) X (if (math-matrixp a) X (if (= (length a) n1) X (math-lud-solve (math-matrix-lud b) a) X (if (= (length (nth 1 a)) n1) X (math-transpose X (math-lud-solve (math-matrix-lud X (math-transpose b)) X (math-transpose a))) X (math-dimension-error))) X (if (= (length a) n1) X (math-mat-col (math-lud-solve (math-matrix-lud b) X (math-col-matrix a)) X 1) X (math-dimension-error))) X (if (Math-equal-int a 1) X (math-inv b) X (math-mul a (math-inv b))))) X (math-reject-arg b 'square-matrixp))) X ((Math-vectorp a) X (math-map-vec-2 'math-div a b)) X ((eq (car-safe a) 'sdev) X (if (eq (car-safe b) 'sdev) X (let ((x (math-div (nth 1 a) (nth 1 b)))) X (math-make-sdev X x X (math-div X (math-hypot (nth 2 a) (math-mul (nth 2 b) x)) X (math-abs (nth 1 b))))) X (and (or (Math-anglep b) X (not (Math-objvecp b))) X (math-make-sdev (math-div (nth 1 a) b) X (math-abs (math-div (nth 2 a) b)))))) X ((and (eq (car-safe b) 'sdev) X (or (Math-anglep a) X (not (Math-objvecp a)))) X (let ((x (math-div a (nth 1 b)))) X (math-make-sdev x X (math-abs (math-div (math-mul (nth 2 b) x) X (nth 1 b)))))) X ((and (eq (car-safe a) 'intv) (Math-anglep b)) X (if (Math-negp b) X (math-neg (math-div a (math-neg b))) X (math-make-intv (nth 1 a) X (math-div (nth 2 a) b) X (math-div (nth 3 a) b)))) X ((and (eq (car-safe b) 'intv) (Math-anglep a)) X (if (Math-posp (nth 2 b)) X (if (Math-negp a) X (math-neg (math-div (math-neg a) b)) X (math-make-intv (aref [0 2 1 3] (nth 1 b)) X (math-div a (nth 3 b)) X (math-div a (nth 2 b)))) X (if (Math-negp (nth 3 b)) X (math-neg (math-div a (math-neg b))) X (calc-record-why "Division by zero" b) X nil))) X ((and (eq (car-safe a) 'intv) (math-constp a) X (eq (car-safe b) 'intv) (math-constp b)) X (if (or (Math-posp (nth 2 b)) (Math-negp (nth 3 b))) X (let ((lo (math-div a (nth 2 b))) X (hi (math-div a (nth 3 b)))) X (and (Math-anglep lo) X (setq lo (list 'intv (if (memq (nth 1 b) '(2 3)) 3 0) X lo lo))) X (and (Math-anglep hi) X (setq hi (list 'intv (if (memq (nth 1 b) '(1 3)) 3 0) X hi hi))) X (math-combine-intervals X (nth 2 lo) (and (memq (nth 1 b) '(2 3)) X (memq (nth 1 lo) '(2 3))) X (nth 3 lo) (and (memq (nth 1 b) '(2 3)) X (memq (nth 1 lo) '(1 3))) X (nth 2 hi) (and (memq (nth 1 b) '(1 3)) X (memq (nth 1 hi) '(2 3))) X (nth 3 hi) (and (memq (nth 1 b) '(1 3)) X (memq (nth 1 hi) '(1 3))))) X (calc-record-why "Division by zero" b) X nil)) X ((and (eq (car-safe a) 'mod) X (eq (car-safe b) 'mod) X (equal (nth 2 a) (nth 2 b))) X (math-make-mod (math-div-mod (nth 1 a) (nth 1 b) (nth 2 a)) X (nth 2 a))) X ((and (eq (car-safe a) 'mod) X (Math-anglep b)) X (math-make-mod (math-div-mod (nth 1 a) b (nth 2 a)) (nth 2 a))) X ((and (eq (car-safe b) 'mod) X (Math-anglep a)) X (math-make-mod (math-div-mod a (nth 1 b) (nth 2 b)) (nth 2 b))) X ((eq (car-safe a) 'hms) X (if (eq (car-safe b) 'hms) X (math-with-extra-prec 1 X (math-div (math-from-hms a 'deg) X (math-from-hms b 'deg))) X (math-with-extra-prec 2 X (math-to-hms (math-div (math-from-hms a 'deg) b) 'deg)))) X (t (calc-record-why "Incompatible arguments" a b))) X) X X(defun math-div-symb-fancy (a b) X (or (and (Math-equal-int b 1) a) X (and (Math-equal-int b -1) (math-neg a)) X (and (eq (car-safe b) '^) X (or (Math-looks-negp (nth 2 b)) (Math-equal-int a 1)) X (math-mul a (math-normalize X (list '^ (nth 1 b) (math-neg (nth 2 b)))))) X (and (eq (car-safe a) 'neg) X (math-neg (math-div (nth 1 a) b))) X (and (eq (car-safe b) 'neg) X (math-neg (math-div a (nth 1 b)))) X (and (eq (car-safe a) '/) X (math-div (nth 1 a) (math-mul (nth 2 a) b))) X (and (eq (car-safe b) '/) X (math-div (math-mul a (nth 2 b)) (nth 1 b))) X (and (eq (car-safe b) 'frac) X (math-mul a (math-make-frac (nth 2 b) (nth 1 b)))) X (and (eq (car-safe a) '+) X (or (Math-numberp (nth 1 a)) X (Math-numberp (nth 2 a))) X (Math-numberp b) X (math-add (math-div (nth 1 a) b) X (math-div (nth 2 a) b))) X (and (eq (car-safe a) '-) X (or (Math-numberp (nth 1 a)) X (Math-numberp (nth 2 a))) X (Math-numberp b) X (math-sub (math-div (nth 1 a) b) X (math-div (nth 2 a) b))) X (and (or (eq (car-safe a) '-) X (math-looks-negp a)) X (math-looks-negp b) X (math-div (math-neg a) (math-neg b))) X (and (eq (car-safe b) '-) X (math-looks-negp a) X (math-div (math-neg a) (math-neg b))) X (if (eq (car-safe a) '*) X (if (eq (car-safe b) '*) X (let ((c (math-combine-prod (nth 1 a) (nth 1 b) nil t t))) X (and c X (math-div (math-mul c (nth 2 a)) (nth 2 b)))) X (let ((c (math-combine-prod (nth 1 a) b nil t t))) X (and c X (math-mul c (nth 2 a))))) X (if (eq (car-safe b) '*) X (let ((c (math-combine-prod a (nth 1 b) nil t t))) X (and c X (math-div c (nth 2 b)))) X (math-combine-prod a b nil t nil))) X (list '/ a b)) X) X X(defun math-div-mod (a b m) ; [R R R R] (Returns nil if no solution) X (and (Math-integerp a) (Math-integerp b) (Math-integerp m) X (let ((u1 1) (u3 b) (v1 0) (v3 m)) X (while (not (eq v3 0)) ; See Knuth sec 4.5.2, exercise 15 X (let* ((q (math-idivmod u3 v3)) X (t1 (math-sub u1 (math-mul v1 (car q))))) X (setq u1 v1 u3 v3 v1 t1 v3 (cdr q)))) X (let ((q (math-idivmod a u3))) X (and (eq (cdr q) 0) X (math-mod (math-mul (car q) u1) m))))) X) X X(defun math-mod-intv (a b) X (let* ((q1 (math-floor (math-div (nth 2 a) b))) X (q2 (math-floor (math-div (nth 3 a) b))) X (m1 (math-sub (nth 2 a) (math-mul q1 b))) X (m2 (math-sub (nth 3 a) (math-mul q2 b)))) X (cond ((equal q1 q2) X (math-sort-intv (nth 1 a) m1 m2)) X ((and (math-equal-int (math-sub q2 q1) 1) X (math-zerop m2) X (memq (nth 1 a) '(0 2))) X (math-make-intv (nth 1 a) m1 b)) X (t X (math-make-intv 2 0 b)))) X) 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 (math-float (math-mul b (math-ln-raw (math-float a))))))))) X ((or (not (Math-objvecp a)) X (not (Math-objectp b))) X (cond ((and (eq (car-safe a) 'calcFunc-sqrt) X (math-evenp b)) X (math-pow (nth 1 a) (math-div2 b))) X ((eq (car-safe a) '*) X (math-mul (math-pow (nth 1 a) b) X (math-pow (nth 2 a) b))) X ((eq (car-safe a) '/) X (math-div (math-pow (nth 1 a) b) X (math-pow (nth 2 a) b))) X ((and (eq (car-safe a) '^) X (Math-integerp b)) X (math-pow (nth 1 a) (math-mul (nth 2 a) b))) X ((and (math-looks-negp a) X (Math-integerp b)) X (if (math-evenp b) X (math-pow (math-neg a) b) X (math-neg (math-pow (math-neg a) b)))) X (t (if (Math-objectp a) X (calc-record-why 'objectp b) X (calc-record-why 'objectp a)) X (list '^ a b)))) X ((and (eq (car-safe a) 'sdev) (eq (car-safe b) 'sdev)) X (if (and (math-constp a) (math-constp b)) X (math-with-extra-prec 2 X (let* ((ln (math-ln-raw (math-float (nth 1 a)))) X (pow (math-exp-raw X (math-float (math-mul (nth 1 b) ln))))) X (list 'sdev X pow X (math-mul X pow X (math-hypot (math-mul (nth 2 a) X (math-div (nth 1 b) X (nth 1 a))) X (math-mul (nth 2 b) ln)))))) X (let ((pow (math-pow (nth 1 a) (nth 1 b)))) X (list 'sdev X pow X (math-mul pow X (math-hypot (math-mul (nth 2 a) X (math-div (nth 1 b) X (nth 1 a))) X (math-mul (nth 2 b) X (math-ln X (nth 1 a))))))))) X ((and (eq (car-safe a) 'sdev) (Math-realp b)) X (if (math-constp a) X (math-with-extra-prec 2 X (let ((pow (math-pow (nth 1 a) (math-sub b 1)))) X (list 'sdev X (math-mul pow (nth 1 a)) X (math-mul pow (math-mul (nth 2 a) b))))) X (list 'sdev X (math-mul (math-pow (nth 1 a) b)) X (math-mul (math-pow (nth 1 a) (math-add b -1)) X (math-mul (nth 2 a) b))))) X ((and (eq (car-safe b) 'sdev) (Math-realp a)) X (math-with-extra-prec 2 X (let* ((ln (math-ln-raw (math-float a))) X (pow (math-exp (math-mul (nth 1 b) ln)))) X (list 'sdev X pow X (math-mul pow (math-mul (nth 2 b) ln)))))) X ((and (eq (car-safe a) 'intv) (math-constp a) X (Math-realp b) X (or (Math-posp (nth 2 a)) X (Math-natnump b) X (and (math-zerop (nth 2 a)) X (Math-posp b)))) X (if (math-evenp b) X (setq a (math-abs a))) X (math-sort-intv (nth 1 a) X (math-pow (nth 2 a) b) X (math-pow (nth 3 a) b))) X ((and (eq (car-safe b) 'intv) (math-constp b) X (Math-posp a)) X (math-sort-intv (nth 1 b) X (math-pow a (nth 2 b)) X (math-pow a (nth 3 b)))) X ((and (eq (car-safe a) 'intv) (math-constp a) X (eq (car-safe b) 'intv) (math-constp b) X (or (and (not (Math-negp (nth 2 a))) X (not (Math-negp (nth 2 b)))) X (and (Math-posp (nth 2 a)) X (not (Math-posp (nth 3 b)))))) X (let ((lo (math-pow a (nth 2 a))) X (hi (math-pow a (nth 3 a)))) X (math-combine-intervals (nth 2 lo) (and (memq (nth 1 a) '(2 3)) X (memq (nth 1 lo) '(2 3))) X (nth 3 lo) (and (memq (nth 1 a) '(2 3)) X (memq (nth 1 lo) '(1 3))) X (nth 2 hi) (and (memq (nth 1 a) '(1 3)) X (memq (nth 1 hi) '(2 3))) X (nth 3 hi) (and (memq (nth 1 a) '(1 3)) X (memq (nth 1 hi) '(1 3)))))) X ((and (eq (car-safe a) 'mod) (eq (car-safe b) 'mod) X (equal (nth 2 a) (nth 2 b))) X (math-make-mod (math-pow-mod (nth 1 a) (nth 1 b) (nth 2 a)) X (nth 2 a))) X ((and (eq (car-safe a) 'mod) (Math-anglep b)) X (math-make-mod (math-pow-mod (nth 1 a) b (nth 2 a)) (nth 2 a))) X ((and (eq (car-safe b) 'mod) (Math-anglep a)) X (math-make-mod (math-pow-mod a (nth 1 b) (nth 2 b)) (nth 2 b))) X ((not (Math-numberp a)) X (math-reject-arg a 'numberp)) X (t X (math-reject-arg b 'numberp))) 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 (if (Math-negp b) X (math-div-mod 1 (math-pow-mod a (Math-integer-neg b) m) m) X (if (eq m 1) X 0 X (math-pow-mod-step a b m))) X (math-mod (math-pow a b) m)) X) X X(defun math-pow-mod-step (a n m) ; [I I I I] X (math-working "pow" a) X (let ((val (cond X ((eq n 0) 1) X ((eq n 1) a) X (t X (let ((rest (math-pow-mod-step X (math-imod (math-mul a a) m) X (math-div2 n) X m))) X (if (math-evenp n) X rest X (math-mod (math-mul a rest) m))))))) X (math-working "pow" val) X val) X) X X(defvar math-power-of-2-cache (list 1 2 4 8 16 32 64 128 256 512 1024)) X(defvar math-big-power-of-2-cache nil) X(defun math-power-of-2 (n) ; [I I] [Public] X (if (and (natnump n) (<= n 100)) X (or (nth n math-power-of-2-cache) X (let* ((i (length math-power-of-2-cache)) X (val (nth (1- i) math-power-of-2-cache))) X (while (<= i n) X (setq val (math-mul val 2) X math-power-of-2-cache (nconc math-power-of-2-cache X (list val)) X i (1+ i))) X val)) X (let ((found (assq n math-big-power-of-2-cache))) X (if found X (cdr found) X (let ((po2 (math-ipow 2 n))) X (setq math-big-power-of-2-cache X (cons (cons n po2) math-big-power-of-2-cache)) X po2)))) X) X X(defun math-integer-log2 (n) ; [I I] [Public] X (let ((i 0) X (p math-power-of-2-cache) X val) X (while (and p (Math-natnum-lessp (setq val (car p)) n)) X (setq p (cdr p) X i (1+ i))) X (if p X (and (equal val n) X i) X (while (Math-natnum-lessp X (prog1 X (setq val (math-mul val 2)) X (setq math-power-of-2-cache (nconc math-power-of-2-cache X (list val)))) X n) X (setq i (1+ i))) X (and (equal val n) X i))) X) X X X X;;; Compute the integer square-root floor(sqrt(A)). A > 0. [I I] [Public] X;;; This method takes advantage of the fact that Newton's method starting 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 (let ((len (length a))) X (if (= (% len 2) 0) X (let* ((top (nthcdr (- len 2) a))) X (math-isqrt-bignum-iter X a X (math-scale-bignum-3 X (math-bignum-big X (1+ (math-isqrt-small X (+ (* (nth 1 top) 1000) (car top))))) X (1- (/ len 2))))) X (let* ((top (nth (1- len) a))) X (math-isqrt-bignum-iter X a X (math-scale-bignum-3 X (list (1+ (math-isqrt-small top))) X (/ len 2)))))) X) X X(defun math-isqrt-bignum-iter (a guess) ; [l L l] X (math-working "isqrt" (cons 'bigpos guess)) X (let* ((q (math-div-bignum a guess)) X (s (math-add-bignum (car q) guess)) X (g2 (math-div2-bignum s)) X (comp (math-compare-bignum g2 guess))) X (if (< comp 0) X (math-isqrt-bignum-iter a g2) X (cons (and (= comp 0) X (math-zerop-bignum (cdr q)) X (= (% (car s) 2) 0)) X guess))) X) X X(defun math-scale-bignum-3 (a n) ; [L L S] X (while (> n 0) X (setq a (cons 0 a) X n (1- n))) X a X) X X(defun math-isqrt-small (a) ; A > 0. [S S] X (let ((g (cond ((>= a 10000) 1000) X ((>= a 100) 100) X (t 10))) X g2) X (while (< (setq g2 (/ (+ g (/ a g)) 2)) g) X (setq g g2)) X g) X) X X X(defun math-inexact-result () X (and calc-symbolic-mode X (signal 'inexact-result nil)) X) X X X;;; Compute the square root of a number. X;;; [T N] if possible, else [F N] if possible, else [C N]. [Public] X(defun math-sqrt (a) X (or X (and (Math-zerop a) a) X (and (Math-negp a) X (math-imaginary (math-sqrt (math-neg a)))) X (and (integerp a) X (let ((sqrt (math-isqrt-small a))) X (if (= (* sqrt sqrt) a) X sqrt X (math-sqrt-float (math-float a) (math-float sqrt))))) X (and (eq (car-safe a) 'bigpos) X (let* ((res (math-isqrt-bignum (cdr a))) X (sqrt (math-normalize (cons 'bigpos (cdr res))))) X (if (car res) X sqrt X (math-sqrt-float (math-float a) (math-float sqrt))))) X (and (eq (car-safe a) 'frac) X (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a))))) X (num-sqrt (math-normalize (cons 'bigpos (cdr num-res)))) X (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a))))) X (den-sqrt (math-normalize (cons 'bigpos (cdr den-res))))) X (if (and (car num-res) (car den-res)) X (list 'frac num-sqrt den-sqrt) X (math-sqrt-float (math-float a) X (math-div (math-float num-sqrt) den-sqrt))))) X (and (eq (car-safe a) 'float) X (if calc-symbolic-mode X (if (= (% (nth 2 a) 2) 0) X (let ((res (math-isqrt-bignum X (cdr (Math-bignum-test (nth 1 a)))))) X (if (car res) X (math-make-float (math-normalize X (cons 'bigpos (cdr res))) X (/ (nth 2 a) 2)) X (signal 'inexact-result nil))) X (signal 'inexact-result nil)) X (math-sqrt-float a))) X (and (eq (car-safe a) 'cplx) X (math-with-extra-prec 2 X (let* ((d (math-abs a)) X (imag (math-sqrt (math-mul (math-sub d (nth 1 a)) X '(float 5 -1))))) X (list 'cplx X (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1))) X (if (math-negp (nth 2 a)) (math-neg imag) imag))))) X (and (eq (car-safe a) 'polar) X (list 'polar X (math-sqrt (nth 1 a)) X (math-mul (nth 2 a) '(float 5 -1)))) X (and (eq (car-safe a) 'sdev) (not (math-negp (nth 1 a))) X (let ((sqrt (math-sqrt (nth 1 a)))) X (math-make-sdev sqrt X (math-div (nth 2 a) (math-mul sqrt 2))))) X (and (eq (car-safe a) 'intv) X (math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a)))) X (and (memq (car-safe a) '(* /)) X (let ((s1 (math-sqrt (nth 1 a))) X (s2 (math-sqrt (nth 2 a)))) X (and (not (and (eq (car-safe s1) 'calcFunc-sqrt) X (eq (car-safe s2) 'calcFunc-sqrt))) X (if (eq (car a) '*) X (math-mul s1 s2) X (math-div s1 s2))))) X (progn X (calc-record-why 'numberp a) X (list 'calcFunc-sqrt a))) X) X(fset 'calcFunc-sqrt (symbol-function 'math-sqrt)) X X(defun math-sqrt-float (a &optional guess) ; [F F F] X (if calc-symbolic-mode X (signal 'inexact-result nil) X (math-with-extra-prec 1 (math-sqrt-raw a guess))) X) X X(defun math-sqrt-raw (a &optional guess) ; [F F F] X (if (not (Math-posp a)) X (math-sqrt a) X (if (null guess) X (let ((ldiff (- (math-numdigs (nth 1 a)) 6))) X (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff))) X (setq guess (math-make-float (math-isqrt-small X (math-scale-int (nth 1 a) (- ldiff))) X (/ (+ (nth 2 a) ldiff) 2))))) X (math-sqrt-float-iter a guess)) X) X X(defun math-sqrt-float-iter (a guess) ; [F F F] X (math-working "sqrt" guess) X (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess)) X '(float 5 -1)))) X (if (math-nearly-equal-float g2 guess) X g2 X (math-sqrt-float-iter a g2))) X) 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(defun math-nearly-zerop-float (a b) X (or (eq (nth 1 a) 0) X (<= (+ (math-numdigs (nth 1 a)) (nth 2 a)) 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(defun math-hypot (a b) X (cond ((Math-zerop a) (math-abs b)) X ((Math-zerop b) (math-abs a)) X ((not (Math-scalarp a)) X (calc-record-why 'scalarp a) X (list 'calcFunc-hypot a b)) X ((not (Math-scalarp b)) X (calc-record-why 'scalarp b) X (list 'calcFunc-hypot a b)) X ((and (Math-realp a) (Math-realp b)) X (math-with-extra-prec 1 X (math-sqrt (math-add (math-sqr a) (math-sqr b))))) X ((eq (car-safe a) 'hms) X (if (eq (car-safe b) 'hms) ; this helps sdev's of hms forms X (math-to-hms (math-hypot (math-from-hms a 'deg) X (math-from-hms b 'deg))) X (math-to-hms (math-hypot (math-from-hms a 'deg) b)))) X ((eq (car-safe b) 'hms) X (math-to-hms (math-hypot a (math-from-hms b 'deg)))) X (t nil)) X) X(fset 'calcFunc-hypot (symbol-function 'math-hypot)) X X(defun calcFunc-sqr (x) X (math-pow x 2) X) X X X X;;; Compute the minimum of two real numbers. [R R R] [Public] X(defun math-min (a b) X (if (and (consp a) (eq (car a) 'intv)) X (if (and (consp b) (eq (car b) 'intv)) X (let ((lo (nth 2 a)) X (lom (memq (nth 1 a) '(2 3))) X (hi (nth 3 a)) X (him (memq (nth 1 a) '(1 3))) X res) X (if (= (setq res (math-compare (nth 2 b) lo)) -1) X (setq lo (nth 2 b) lom (memq (nth 1 b) '(2 3))) X (if (= res 0) X (setq lom (or lom (memq (nth 1 b) '(2 3)))))) X (if (= (setq res (math-compare (nth 3 b) hi)) -1) X (setq hi (nth 3 b) him (memq (nth 1 b) '(1 3))) X (if (= res 0) X (setq him (or him (memq (nth 1 b) '(1 3)))))) X (math-make-intv (+ (if lom 2 0) (if him 1 0)) lo hi)) X (math-min a (list 'intv 3 b b))) X (if (and (consp b) (eq (car b) 'intv)) X (math-min (list 'intv 3 a a) b) X (if (Math-lessp a b) X a X b))) X) X X(defun calcFunc-min (a &rest b) X (if (not (or (Math-anglep a) X (and (eq (car a) 'intv) (math-constp a)))) X (math-reject-arg a 'anglep)) X (math-min-list a b) X) X X(defun math-min-list (a b) X (if b X (if (or (Math-anglep (car b)) X (and (eq (car (car b)) 'intv) (math-constp (car b)))) X (math-min-list (math-min a (car b)) (cdr b)) X (math-reject-arg (car b) 'anglep)) X a) X) X X;;; Compute the maximum of two real numbers. [R R R] [Public] X(defun math-max (a b) X (if (or (and (consp a) (eq (car a) 'intv)) X (and (consp b) (eq (car b) 'intv))) X (math-neg (math-min (math-neg a) (math-neg b))) X (if (Math-lessp a b) X b X a)) X) X X(defun calcFunc-max (a &rest b) X (if (not (or (Math-anglep a) X (and (eq (car a) 'intv) (math-constp a)))) X (math-reject-arg a 'anglep)) X (math-max-list a b) X) X X(defun math-max-list (a b) X (if b X (if (or (Math-anglep (car b)) X (and (eq (car (car b)) 'intv) (math-constp (car b)))) X (math-max-list (math-max a (car b)) (cdr b)) X (math-reject-arg (car b) 'anglep)) X a) X) X X X;;; Compute the absolute value of A. [O O; r r] [Public] X(defun math-abs (a) X (cond ((Math-negp a) X (math-neg a)) X ((Math-anglep a) X a) X ((eq (car a) 'cplx) X (math-hypot (nth 1 a) (nth 2 a))) X ((eq (car a) 'polar) X (nth 1 a)) X ((eq (car a) 'vec) X (if (cdr (cdr (cdr a))) X (math-sqrt (math-abssqr a)) X (if (cdr (cdr a)) X (math-hypot (nth 1 a) (nth 2 a)) X (if (cdr a) X (math-abs (nth 1 a)) X a)))) X ((eq (car a) 'sdev) X (list 'sdev (math-abs (nth 1 a)) (nth 2 a))) X ((and (eq (car a) 'intv) (math-constp a)) X (if (Math-posp a) X a X (let* ((nlo (math-neg (nth 2 a))) X (res (math-compare nlo (nth 3 a)))) X (cond ((= res 1) X (math-make-intv (if (memq (nth 1 a) '(0 1)) 2 3) 0 nlo)) X ((= res 0) X (math-make-intv (if (eq (nth 1 a) 0) 2 3) 0 nlo)) X (t X (math-make-intv (if (memq (nth 1 a) '(0 2)) 2 3) X 0 (nth 3 a))))))) X ((eq (car a) 'calcFunc-abs) X (car a)) X ((math-looks-negp a) X (list 'calcFunc-abs (math-neg a))) X (t (calc-record-why 'numvecp a) X (list 'calcFunc-abs a))) X) X(fset 'calcFunc-abs (symbol-function 'math-abs)) 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 (if (math-messy-integerp (nth 2 a)) X (math-trunc (math-make-mod (nth 1 a) (math-trunc (nth 2 a)))) X (math-make-mod (math-trunc (nth 1 a)) (nth 2 a)))) X ((eq (car a) 'intv) X (math-make-intv 3 X (if (and (Math-negp (nth 2 a)) X (Math-num-integerp (nth 2 a)) X (memq (nth 1 a) '(0 1))) X (math-add (math-trunc (nth 2 a)) 1) X (math-trunc (nth 2 a))) X (if (and (Math-posp (nth 3 a)) X (Math-num-integerp (nth 3 a)) X (memq (nth 1 a) '(0 2))) X (math-add (math-trunc (nth 3 a)) -1) X (math-trunc (nth 3 a))))) X ((math-provably-integerp a) a) X (t (math-reject-arg a 'numberp))) X) X(defun calcFunc-ftrunc (a) X (math-float (math-trunc a)) X) X X(defun math-floor-fancy (a) X (cond ((math-provably-integerp a) a) X ((eq (car a) 'hms) X (if (or (math-posp a) X (and (math-zerop (nth 2 a)) X (math-zerop (nth 3 a)))) X (math-trunc a) X (math-add (math-trunc a) -1))) X ((eq (car a) 'intv) X (math-make-intv 3 X (math-floor (nth 2 a)) X (if (and (Math-num-integerp (nth 3 a)) X (memq (nth 1 a) '(0 2))) X (math-add (math-floor (nth 3 a)) -1) X (math-floor (nth 3 a))))) X (t (math-reject-arg a 'anglep))) X) X(defun calcFunc-ffloor (a) X (math-float (math-floor a)) X) X X;;; Coerce A to be an integer (by truncation toward plus infinity). [I N] X(defun math-ceiling (a) ; [Public] X (cond ((Math-integerp a) a) X ((Math-messy-integerp a) (math-trunc a)) X ((Math-realp a) X (if (Math-posp a) X (math-add (math-trunc a) 1) X (math-trunc a))) X ((math-provably-integerp a) a) X ((eq (car a) 'hms) X (if (or (math-negp a) X (and (math-zerop (nth 2 a)) X (math-zerop (nth 3 a)))) X (math-trunc a) X (math-add (math-trunc a) 1))) X ((eq (car a) 'intv) X (math-make-intv 3 X (if (and (Math-num-integerp (nth 2 a)) X (memq (nth 1 a) '(0 1))) X (math-add (math-floor (nth 2 a)) 1) X (math-ceiling (nth 2 a))) X (math-ceiling (nth 3 a)))) X (t (math-reject-arg a 'anglep))) X) X(fset 'calcFunc-ceil (symbol-function 'math-ceiling)) X(defun calcFunc-fceil (a) X (math-float (math-ceiling a)) X) X X;;; Coerce A to be an integer (by rounding to nearest integer). [I N] [Public] X(defun math-round (a) X (cond ((Math-anglep a) X (if (Math-num-integerp a) X (math-trunc a) X (if (Math-negp a) X (math-neg (math-round (math-neg a))) X (math-floor X (let ((calc-angle-mode 'deg)) ; in case of HMS forms X (math-add a (if (Math-ratp a) X '(frac 1 2) X '(float 5 -1)))))))) X ((math-provably-integerp a) a) X ((eq (car a) 'intv) X (math-floor (math-add a '(frac 1 2)))) X (t (math-reject-arg a 'anglep))) X) X(fset 'calcFunc-round (symbol-function 'math-round)) X(defun calcFunc-fround (a) X (math-float (math-round a)) X) X X X;;; Convert a real value to fractional form. [T R I; T R F] [Public] X(defun math-to-fraction (a &optional tol) X (or tol (setq tol 0)) X (cond ((Math-ratp a) X a) X ((memq (car a) '(cplx polar vec hms sdev intv mod)) X (cons (car a) (mapcar (function X (lambda (x) X (math-to-fraction x tol))) X (cdr a)))) X ((Math-negp a) X (math-neg (math-to-fraction (math-neg a) tol))) X ((not (eq (car a) 'float)) X (math-reject-arg a 'numberp)) X ((integerp tol) X (if (<= tol 0) X (setq tol (+ tol calc-internal-prec))) X (math-to-fraction a (list 'float 5 SHAR_EOF echo "End of part 7" echo "File calc-ext.el is continued in part 8" echo "8" > s2_seq_.tmp exit 0