daveg@csvax.cs.caltech.edu (David Gillespie) (10/15/90)
Posting-number: Volume 15, Issue 36
Submitted-by: daveg@csvax.cs.caltech.edu (David Gillespie)
Archive-name: calc-1.05/part09
#!/bin/sh
# this is part 9 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file calc.patch continued
#
CurArch=9
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+ (math-bernoulli-number n)))
X+ )
X+
X+ (defun calcFunc-euler (n &optional x)
X+ (or (math-num-natnump n) (math-reject-arg n 'natnump))
X+ (if x
X+ (let* ((n1 (math-add n 1))
X+ (coefs (math-bernoulli-coefs n1))
X+ (fac (math-div (math-pow 2 n1) n1))
X+ (k -1)
X+ (x1 (math-div (math-add x 1) 2))
X+ (x2 (math-div x 2)))
X+ (if (math-numberp x)
X+ (if (and calc-symbolic-mode (math-floatp x))
X+ (math-inexact-result)
X+ (math-mul fac
X+ (math-sub (math-build-polynomial-expr coefs x1)
X+ (math-build-polynomial-expr coefs x2))))
X+ (math-collect-terms
X+ (math-reduce-vec
X+ 'math-add
X+ (cons 'vec
X+ (mapcar (function
X+ (lambda (c)
X+ (setq k (1+ k))
X+ (math-mul (math-mul fac c)
X+ (math-sub (math-pow x1 k)
X+ (math-pow x2 k)))))
X+ coefs)))
X+ x)))
X+ (math-mul (math-pow 2 n)
X+ (if (consp n)
X+ (progn
X+ (math-inexact-result)
X+ (calcFunc-euler n '(float 5 -1)))
X+ (calcFunc-euler n '(frac 1 2)))))
X+ )
X+
X+ (defun math-bernoulli-coefs (n)
X+ (let* ((coefs (list (calcFunc-bern n)))
X+ (nn (math-trunc n))
X+ (k nn)
X+ (term nn)
X+ coef
X+ (calc-prefer-frac (or (integerp n) calc-prefer-frac)))
X+ (while (>= (setq k (1- k)) 0)
X+ (setq term (math-div term (- nn k))
X+ coef (math-mul term (math-bernoulli-number k))
X+ coefs (cons (if (consp n) (math-float coef) coef) coefs)
X+ term (math-mul term k)))
X+ (nreverse coefs))
X+ )
X+
X+ (defun math-bernoulli-number (n)
X+ (if (= (% n 2) 1)
X+ (if (= n 1)
X+ '(frac -1 2)
X+ 0)
X+ (setq n (/ n 2))
X+ (while (>= n math-bernoulli-cache-size)
X+ (let* ((sum 0)
X+ (nk 1) ; nk = n-k+1
X+ (fact 1) ; fact = (n-k+1)!
X+ ofact
X+ (p math-bernoulli-b-cache)
X+ (calc-prefer-frac t))
X+ (math-working "bernoulli B" (* 2 math-bernoulli-cache-size))
X+ (while p
X+ (setq nk (+ nk 2)
X+ ofact fact
X+ fact (math-mul fact (* nk (1- nk)))
X+ sum (math-add sum (math-div (car p) fact))
X+ p (cdr p)))
X+ (setq ofact (math-mul ofact (1- nk))
X+ sum (math-sub (math-div '(frac 1 2) ofact) sum)
X+ math-bernoulli-b-cache (cons sum math-bernoulli-b-cache)
X+ math-bernoulli-B-cache (cons (math-mul sum ofact)
X+ math-bernoulli-B-cache)
X+ math-bernoulli-cache-size (1+ math-bernoulli-cache-size))))
X+ (nth (- math-bernoulli-cache-size n 1) math-bernoulli-B-cache))
X+ )
X+
X+ ;;; Bn = n! bn
X+ ;;; bn = - sum_k=0^n-1 bk / (n-k+1)!
X+
X+ ;;; A faster method would be to use "tangent numbers", c.f., Concrete
X+ ;;; Mathematics pg. 273.
X+
X+ (setq math-bernoulli-b-cache '( (frac -174611
X+ (bigpos 0 200 291 698 662 857 802))
X+ (frac 43867 (bigpos 0 944 170 217 94 109 5))
X+ (frac -3617 (bigpos 0 880 842 622 670 10))
X+ (frac 1 (bigpos 600 249 724 74))
X+ (frac -691 (bigpos 0 368 674 307 1))
X+ (frac 1 (bigpos 160 900 47))
X+ (frac -1 (bigpos 600 209 1))
X+ (frac 1 30240) (frac -1 720)
X+ (frac 1 12) 1 ))
X+
X+ (setq math-bernoulli-B-cache '( (frac -174611 330) (frac 43867 798)
X+ (frac -3617 510) (frac 7 6) (frac -691 2730)
X+ (frac 5 66) (frac -1 30) (frac 1 42)
X+ (frac -1 30) (frac 1 6) 1 ))
X+
X+ (setq math-bernoulli-cache-size 11)
X+
X+
X+
X+ ;;; Probability distributions.
X+
X+ ;;; Binomial.
X+ (defun calcFunc-utpb (x n p)
X+ (calcFunc-betaI p x (math-add (math-sub n x) 1))
X+ )
X+
X+ (defun calcFunc-ltpb (x n p)
X+ (math-sub 1 (calcFunc-betaI p x (math-add (math-sub n x) 1)))
X+ )
X+
X+ ;;; Chi-square.
X+ (defun calcFunc-utpc (chisq v)
X+ (calcFunc-gammaQ (math-div v 2) (math-div chisq 2))
X+ )
X+
X+ (defun calcFunc-ltpc (chisq v)
X+ (calcFunc-gammaP (math-div v 2) (math-div chisq 2))
X+ )
X+
X+ ;;; F-distribution.
X+ (defun calcFunc-utpf (f v1 v2)
X+ (calcFunc-betaI (math-div v2 (math-add v2 (math-mul v1 f)))
X+ (math-div v2 2)
X+ (math-div v1 2))
X+ )
X+
X+ (defun calcFunc-ltpf (f v1 v2)
X+ (math-sub 1 (calcFunc-utpf f v1 v2))
X+ )
X+
X+ ;;; Normal.
X+ (defun calcFunc-utpn (x mean sdev)
X+ (math-mul (math-add '(float 1 0)
X+ (calcFunc-erf (math-div (math-sub mean x)
X+ (math-mul sdev (math-sqrt-2)))))
X+ '(float 5 -1))
X+ )
X+
X+ (defun calcFunc-ltpn (x mean sdev)
X+ (math-mul (math-add '(float 1 0)
X+ (calcFunc-erf (math-div (math-sub x mean)
X+ (math-mul sdev (math-sqrt-2)))))
X+ '(float 5 -1))
X+ )
X+
X+ ;;; Poisson.
X+ (fset 'calcFunc-utpp (symbol-function 'calcFunc-gammaP))
X+ (fset 'calcFunc-ltpp (symbol-function 'calcFunc-gammaQ))
X+
X+ ;;; Student's t. (As defined in Abramowitz & Stegun and Numerical Recipes.)
X+ (defun calcFunc-utpt (tt v)
X+ (calcFunc-betaI (math-div v (math-add v (math-sqr tt)))
X+ (math-div v 2)
X+ '(float 5 -1))
X+ )
X+
X+ (defun calcFunc-ltpt (tt v)
X+ (math-sub 1 (calcFunc-utpt tt v))
X+ )
X+
X+
X+
X+
X ;;;; [calc-arith.el]
X
X ;;;; Number theory.
X***************
X*** 10112,10134 ****
X
X (defun math-factorial (n) ; [I I] [F F] [Public]
X (let (temp)
X! (cond ((Math-integer-negp n) (list 'calcFunc-fact n))
X! ((Math-zerop n) 1)
X! ((integerp n) (math-factorial-iter (1- n) 2 1))
X ((and (math-messy-integerp n)
X (Math-lessp (setq temp (math-trunc n)) 100))
X! (if (natnump temp)
X! (math-with-extra-prec 1
X! (math-factorial-iter (1- temp) 2 '(float 1 0)))
X! (list 'calcFunc-fact max)))
X! ((math-realp n)
X! (math-with-extra-prec 3
X! (math-gammap1-raw (math-float n))))
X! (t (calc-record-why 'realp n)
X (list 'calcFunc-fact n))))
X )
X (fset 'calcFunc-fact (symbol-function 'math-factorial))
X
X (defun math-factorial-iter (count n f)
X (if (= (% n 5) 1)
X (math-working (format "factorial(%d)" (1- n)) f))
X--- 15575,15636 ----
X
X (defun math-factorial (n) ; [I I] [F F] [Public]
X (let (temp)
X! (cond ((Math-integer-negp n)
X! (math-reject-arg n 'natnump))
X! ((integerp n)
X! (if (<= n 20)
X! (aref '[1 1 2 6 24 120 720 5040 40320 362880
X! (bigpos 800 628 3) (bigpos 800 916 39)
X! (bigpos 600 1 479) (bigpos 800 20 227 6)
X! (bigpos 200 291 178 87) (bigpos 0 368 674 307 1)
X! (bigpos 0 888 789 922 20) (bigpos 0 96 428 687 355)
X! (bigpos 0 728 705 373 402 6)
X! (bigpos 0 832 408 100 645 121)
X! (bigpos 0 640 176 8 902 432 2)] n)
X! (math-factorial-iter (1- n) 2 1)))
X ((and (math-messy-integerp n)
X (Math-lessp (setq temp (math-trunc n)) 100))
X! (math-inexact-result)
X! (if (>= temp 0)
X! (if (<= temp 20)
X! (math-float (math-factorial temp))
X! (math-with-extra-prec 1
X! (math-factorial-iter (1- temp) 2 '(float 1 0))))
X! (math-reject-arg n 'natnump)))
X! ((math-numberp n)
X! (let* ((q (math-quarter-integer n))
X! (tn (and q (math-floor n))))
X! (if (and tn (>= (setq tn (1+ tn)) 0) (< tn 20))
X! (if (= q 2)
X! (math-div
X! (math-mul (math-double-factorial-iter (* 2 tn) 3 1 2)
X! (if calc-symbolic-mode
X! (list 'calcFunc-sqrt '(var pi var-pi))
X! (math-sqrt-pi)))
X! (math-pow 2 tn))
X! (math-inexact-result)
X! (math-div
X! (math-mul (math-double-factorial-iter (* 4 tn) q 1 4)
X! (if (= q 1) (math-gamma-1q) (math-gamma-3q)))
X! (math-pow 4 tn)))
X! (math-inexact-result)
X! (math-with-extra-prec 3
X! (math-gammap1-raw (math-float n)
X! (math-float calc-internal-prec)
X! (math-float (- calc-internal-prec)))))))
X! (t (calc-record-why 'numberp n)
X (list 'calcFunc-fact n))))
X )
X (fset 'calcFunc-fact (symbol-function 'math-factorial))
X
X+ (math-defcache math-gamma-1q nil
X+ (math-with-extra-prec 3
X+ (math-gammap1-raw '(float -75 -2))))
X+
X+ (math-defcache math-gamma-3q nil
X+ (math-with-extra-prec 3
X+ (math-gammap1-raw '(float -25 -2))))
X+
X (defun math-factorial-iter (count n f)
X (if (= (% n 5) 1)
X (math-working (format "factorial(%d)" (1- n)) f))
X***************
X*** 10137,10219 ****
X f)
X )
X
X- (math-defcache math-sqrt-two-pi nil
X- (math-sqrt (math-two-pi)))
X-
X- (defun math-gammap1-raw (x) ; compute gamma(1 + x)
X- (cond ((math-lessp-float x '(float 1 1))
X- (if (math-lessp-float x '(float -10 0))
X- (setq x (math-neg-float
X- (math-div-float
X- (math-pi)
X- (math-mul-float (math-gammap1-raw
X- (math-add-float (math-neg-float x)
X- '(float -1 0)))
X- (math-sin-raw
X- (math-mul (math-pi) x)))))))
X- (let ((xplus1 (math-add-float x '(float 1 0))))
X- (math-div-float (math-gammap1-raw xplus1) xplus1)))
X- (t ; x now >= 10.0
X- (let ((xinv (math-div 1 x))
X- (lnx (math-ln-raw x)))
X- (math-mul (math-sqrt-two-pi)
X- (math-exp-raw
X- (math-gamma-series
X- (math-sub (math-mul (math-add x '(float 5 -1))
X- lnx)
X- x)
X- xinv
X- (math-sqr xinv)
X- 2))))))
X- )
X-
X- (defun calcFunc-gamma (x)
X- (calcFunc-fact (math-add x -1))
X- )
X-
X- (defun math-gamma-series (sum x xinvsqr n)
X- (math-working "gamma" sum)
X- (let* ((bn (math-bernoulli-number n)) ; this will always be a "frac" form.
X- (next (math-add-float
X- sum
X- (math-mul-float (math-div-float (math-float (nth 1 bn))
X- (math-float (* (nth 2 bn)
X- (* n (1- n)))))
X- x))))
X- (if (math-nearly-equal-float sum next)
X- next
X- (if (= n 24)
X- (progn
X- (calc-record-why "Gamma computation stopped early, not all digits may be valid")
X- next)
X- (math-gamma-series next (math-mul-float x xinvsqr) xinvsqr (+ n 2)))))
X- )
X-
X- (defun math-bernoulli-number (n)
X- (if (= n 1)
X- '(frac -1 2)
X- (if (= (% n 2) 1)
X- 0
X- (aref '[ 1 (frac 1 6) (frac -1 30) (frac 1 42) (frac -1 30)
X- (frac 5 66) (frac -691 2730) (frac 7 6) (frac -3617 510)
X- (frac 43867 798) (frac -174611 330) (frac 854513 138)
X- (frac (bigneg 91 364 236) 2730) ]
X- (/ n 2))))
X- )
X- ;;; To come up with more, we could use this rule:
X- ;;; Bn = n! bn
X- ;;; bn = - sum_k=0^n-1 bk / (n-k+1)!
X-
X (defun math-double-factorial (n) ; [I I] [F F] [Public]
X (cond ((Math-integer-negp n) (list 'calcFunc-dfact n))
X ((Math-zerop n) 1)
X! ((integerp n) (math-double-factorial-iter n (+ 2 (% n 2)) 1))
X ((math-messy-integerp n)
X (let ((temp (math-trunc n)))
X (if (natnump temp)
X! (math-with-extra-prec 1
X! (math-double-factorial-iter temp (+ 2 (% temp 2))
X! '(float 1 0)))
X (list 'calcFunc-dfact max))))
X (t (calc-record-why 'natnump n)
X (list 'calcFunc-dfact n)))
X--- 15639,15662 ----
X f)
X )
X
X (defun math-double-factorial (n) ; [I I] [F F] [Public]
X (cond ((Math-integer-negp n) (list 'calcFunc-dfact n))
X ((Math-zerop n) 1)
X! ((integerp n) (math-double-factorial-iter n (+ 2 (% n 2)) 1 2))
X ((math-messy-integerp n)
X (let ((temp (math-trunc n)))
X+ (math-inexact-result)
X (if (natnump temp)
X! (if (Math-lessp temp 200)
X! (math-with-extra-prec 1
X! (math-double-factorial-iter temp (+ 2 (% temp 2))
X! '(float 1 0) 2))
X! (let* ((half (math-div2 temp))
X! (even (math-mul (math-pow 2 half)
X! (math-factorial (math-float half)))))
X! (if (math-evenp temp)
X! even
X! (math-div (math-factorial n) even))))
X (list 'calcFunc-dfact max))))
X (t (calc-record-why 'natnump n)
X (list 'calcFunc-dfact n)))
X***************
X*** 10220,10236 ****
X )
X (fset 'calcFunc-dfact (symbol-function 'math-double-factorial))
X
X! (defun math-double-factorial-iter (max n f)
X! (if (< (% n 10) 2)
X! (math-working (format "dfact(%d)" (- n 2)) f))
X (if (<= n max)
X! (math-double-factorial-iter max (+ n 2) (math-mul n f))
X f)
X )
X
X (defun math-permutations (n m) ; [I I I] [F F F] [Public]
X (cond ((and (integerp n) (integerp m) (<= m n) (>= m 0))
X! (math-factorial-iter n (1+ (- n m)) 1))
X ((or (not (math-num-integerp n))
X (not (math-num-integerp m)))
X (or (math-realp n) (math-reject-arg n 'realp))
X--- 15663,15679 ----
X )
X (fset 'calcFunc-dfact (symbol-function 'math-double-factorial))
X
X! (defun math-double-factorial-iter (max n f step)
X! (if (< (% n 12) step)
X! (math-working (format "dfact(%d)" (- n step)) f))
X (if (<= n max)
X! (math-double-factorial-iter max (+ n step) (math-mul n f) step)
X f)
X )
X
X (defun math-permutations (n m) ; [I I I] [F F F] [Public]
X (cond ((and (integerp n) (integerp m) (<= m n) (>= m 0))
X! (math-factorial-iter m (1+ (- n m)) 1))
X ((or (not (math-num-integerp n))
X (not (math-num-integerp m)))
X (or (math-realp n) (math-reject-arg n 'realp))
X***************
X*** 10237,10246 ****
X (or (math-realp m) (math-reject-arg m 'realp))
X (and (math-num-integerp n) (math-negp n) (math-reject-arg n 'range))
X (and (math-num-integerp m) (math-negp m) (math-reject-arg m 'range))
X! (math-div (math-factorial n) (math-factorial m)))
X (t
X (let ((tn (math-trunc n))
X (tm (math-trunc m)))
X (or (integerp tn) (math-reject-arg tn 'fixnump))
X (or (integerp tm) (math-reject-arg tm 'fixnump))
X (or (and (<= tm tn) (>= tm 0)) (math-reject-arg tm 'range))
X--- 15680,15690 ----
X (or (math-realp m) (math-reject-arg m 'realp))
X (and (math-num-integerp n) (math-negp n) (math-reject-arg n 'range))
X (and (math-num-integerp m) (math-negp m) (math-reject-arg m 'range))
X! (math-div (math-factorial n) (math-factorial (math-sub n m))))
X (t
X (let ((tn (math-trunc n))
X (tm (math-trunc m)))
X+ (math-inexact-result)
X (or (integerp tn) (math-reject-arg tn 'fixnump))
X (or (integerp tm) (math-reject-arg tm 'fixnump))
X (or (and (<= tm tn) (>= tm 0)) (math-reject-arg tm 'range))
X***************
X*** 10274,10279 ****
X--- 15718,15724 ----
X (Math-lessp n m))
X 0)
X (t
X+ (math-inexact-result)
X (let ((tm (math-trunc m)))
X (or (integerp tm) (math-reject-arg tm 'fixnump))
X (if (> tm 100)
X***************
X*** 10282,10288 ****
X (math-factorial (math-float
X (math-sub n m)))))
X (math-with-extra-prec 1
X! (math-choose-float-iter tm n 1 '(float 1 0)))))))
X )
X (fset 'calcFunc-choose (symbol-function 'math-choose))
X
X--- 15727,15733 ----
X (math-factorial (math-float
X (math-sub n m)))))
X (math-with-extra-prec 1
X! (math-choose-float-iter tm n 1 1))))))
X )
X (fset 'calcFunc-choose (symbol-function 'math-choose))
X
X***************
X*** 10305,10310 ****
X--- 15750,15805 ----
X )
X
X
X+ ;;; Stirling numbers.
X+
X+ (defun calcFunc-stir1 (n m)
X+ (math-stirling-number n m 1)
X+ )
X+
X+ (defun calcFunc-stir2 (n m)
X+ (math-stirling-number n m 0)
X+ )
X+
X+ (defun math-stirling-number (n m k)
X+ (or (math-num-natnump n) (math-reject-arg n 'natnump))
X+ (or (math-num-natnump m) (math-reject-arg m 'natnump))
X+ (if (consp n) (setq n (math-trunc n)))
X+ (or (integerp n) (math-reject-arg n 'fixnump))
X+ (if (consp m) (setq m (math-trunc m)))
X+ (or (integerp m) (math-reject-arg m 'fixnump))
X+ (if (< n m)
X+ 0
X+ (let ((cache (aref math-stirling-cache k)))
X+ (while (<= (length cache) n)
X+ (let ((i (1- (length cache)))
X+ row)
X+ (setq cache (vconcat cache (make-vector (length cache) nil)))
X+ (aset math-stirling-cache k cache)
X+ (while (< (setq i (1+ i)) (length cache))
X+ (aset cache i (setq row (make-vector (1+ i) nil)))
X+ (aset row 0 0)
X+ (aset row i 1))))
X+ (if (= k 1)
X+ (math-stirling-1 n m)
X+ (math-stirling-2 n m))))
X+ )
X+ (setq math-stirling-cache (vector [[1]] [[1]]))
X+
X+ (defun math-stirling-1 (n m)
X+ (or (aref (aref cache n) m)
X+ (aset (aref cache n) m
X+ (math-add (math-stirling-1 (1- n) (1- m))
X+ (math-mul (- 1 n) (math-stirling-1 (1- n) m)))))
X+ )
X+
X+ (defun math-stirling-2 (n m)
X+ (or (aref (aref cache n) m)
X+ (aset (aref cache n) m
X+ (math-add (math-stirling-2 (1- n) (1- m))
X+ (math-mul m (math-stirling-2 (1- n) m)))))
X+ )
X+
X+
X ;;; Produce a random digit in the range 0..999.
X ;;; Avoid various pitfalls that may lurk in the built-in (random) function!
X (defun math-random-digit ()
X***************
X*** 10386,10395 ****
X--- 15881,15979 ----
X (if (Math-lessp lo hi)
X (math-add (math-random (math-sub hi lo)) lo)
X (math-reject-arg max "Empty interval")))))
X+ ((eq (car max) 'vec)
X+ (if (cdr max)
X+ (nth (1+ (math-random (1- (length max)))) max)
X+ (math-reject-arg max "Empty list")))
X+ ((and (eq (car max) 'sdev) (math-constp max))
X+ (math-add (math-mul (math-gaussian-float) (nth 2 max)) (nth 1 max)))
X (t (math-reject-arg max 'realp)))
X )
X (fset 'calcFunc-random (symbol-function 'math-random))
X
X+ ;;; Choose N objects at random from the set MAX without duplicates.
X+ (defun math-shuffle (n max)
X+ (or (and (Math-num-integerp n)
X+ (or (natnump (setq n (math-trunc n))) (eq n -1)))
X+ (math-reject-arg n 'integerp))
X+ (cond ((or (math-zerop max)
X+ (math-floatp max)
X+ (eq (car-safe max) 'sdev))
X+ (if (< n 0)
X+ (math-reject-arg n 'natnump)
X+ (math-simple-shuffle n max)))
X+ ((and (<= n 1) (>= n 0))
X+ (math-simple-shuffle n max))
X+ ((eq (car-safe max) 'intv)
X+ (let ((num (math-add (math-sub (nth 3 max) (nth 2 max))
X+ (cdr (assq (nth 1 max)
X+ '((0 . -1) (1 . 0)
X+ (2 . 0) (3 . 1))))))
X+ (min (math-add (nth 2 max) (if (memq (nth 1 max) '(0 1))
X+ 1 0))))
X+ (if (< n 0) (setq n num))
X+ (or (math-posp num) (math-reject-arg max 'range))
X+ (and (math-lessp num n) (math-reject-arg n 'range))
X+ (if (math-lessp n (math-quotient num 3))
X+ (math-simple-shuffle n max)
X+ (if (> (* n 4) (* num 3))
X+ (math-add (math-sub min 1)
X+ (math-shuffle-list n num (math-vec-index num)))
X+ (let ((tot 0)
X+ (m 0)
X+ (vec nil))
X+ (while (< m n)
X+ (if (< (math-random (- num tot)) (- n m))
X+ (setq vec (cons (math-add min tot) vec)
X+ m (1+ m)))
X+ (setq tot (1+ tot)))
X+ (math-shuffle-list n n (cons 'vec vec)))))))
X+ ((eq (car-safe max) 'vec)
X+ (let ((size (1- (length max))))
X+ (if (< n 0) (setq n size))
X+ (if (and (> n (/ size 2)) (<= n size))
X+ (math-shuffle-list n size (copy-sequence max))
X+ (let* ((vals (math-shuffle n (list 'intv 3 1 (1- (length max)))))
X+ (p vals))
X+ (while (setq p (cdr p))
X+ (setcar p (nth (car p) max)))
X+ vals))))
X+ ((math-integerp max)
X+ (if (math-posp max)
X+ (math-shuffle n (list 'intv 2 0 max))
X+ (math-shuffle n (list 'intv 1 max 0))))
X+ (t (math-reject-arg max 'realp)))
X+ )
X+ (fset 'calcFunc-shuffle (symbol-function 'math-shuffle))
X+
X+ (defun math-simple-shuffle (n max)
X+ (let ((vec nil)
X+ val)
X+ (while (>= (setq n (1- n)) 0)
X+ (while (math-member (setq val (math-random max)) vec))
X+ (setq vec (cons val vec)))
X+ (cons 'vec vec))
X+ )
X+
X+ (defun math-shuffle-list (n size vec)
X+ (let ((j size)
X+ k temp
X+ (p vec))
X+ (while (cdr (setq p (cdr p)))
X+ (setq k (math-random j)
X+ j (1- j)
X+ temp (nth k p))
X+ (setcar (nthcdr k p) (car p))
X+ (setcar p temp))
X+ (cons 'vec (nthcdr (- size n -1) vec)))
X+ )
X+
X+ (defun math-member (x list)
X+ (while (and list (not (equal x (car list))))
X+ (setq list (cdr list)))
X+ list
X+ )
X+
X
X ;;; Check if the integer N is prime. [X I]
X ;;; Return (nil) if non-prime,
X***************
X*** 10484,10489 ****
X--- 16068,16081 ----
X )
X (defvar math-prime-test-cache '(-1))
X
X+ (defun calcFunc-prime (n &optional iters)
X+ (or (math-num-integerp n) (math-reject-arg 'integerp n))
X+ (or (not iters) (math-num-integerp iters) (math-reject-arg 'integerp iters))
X+ (if (car (math-prime-test (math-trunc n) (math-trunc (or iters 1))))
X+ 1
X+ 0)
X+ )
X+
X ;;; Theory: summing base-10^6 digits modulo 111111 is "casting out 999999s".
X ;;; Initial probability that N is prime is 1/ln(N) = log10(e)/log10(N).
X ;;; After culling [2,3,5,7,11,13,37], probability of primality is 5.36 x more.
X***************
X*** 10572,10578 ****
X (fset 'calcFunc-moebius (symbol-function 'math-moebius))
X
X
X! (defun math-next-prime (n iters)
X (if (Math-integerp n)
X (if (Math-integer-negp n)
X 2
X--- 16164,16170 ----
X (fset 'calcFunc-moebius (symbol-function 'math-moebius))
X
X
X! (defun math-next-prime (n &optional iters)
X (if (Math-integerp n)
X (if (Math-integer-negp n)
X 2
X***************
X*** 10582,10588 ****
X (setq n (math-add n -1)))
X (let (res)
X (while (not (car (setq res (math-prime-test
X! (setq n (math-add n 2)) iters)))))
X (if (and calc-verbose-nextprime
X (eq (car res) 'maybe))
X (calc-report-prime-test res)))
X--- 16174,16181 ----
X (setq n (math-add n -1)))
X (let (res)
X (while (not (car (setq res (math-prime-test
X! (setq n (math-add n 2))
X! (or iters 1))))))
X (if (and calc-verbose-nextprime
X (eq (car res) 'maybe))
X (calc-report-prime-test res)))
X***************
X*** 10594,10600 ****
X (fset 'calcFunc-nextprime (symbol-function 'math-next-prime))
X (setq calc-verbose-nextprime nil)
X
X! (defun math-prev-prime (n iters)
X (if (Math-integerp n)
X (if (Math-lessp n 4)
X 2
X--- 16187,16193 ----
X (fset 'calcFunc-nextprime (symbol-function 'math-next-prime))
X (setq calc-verbose-nextprime nil)
X
X! (defun math-prev-prime (n &optional iters)
X (if (Math-integerp n)
X (if (Math-lessp n 4)
X 2
X***************
X*** 10602,10608 ****
X (setq n (math-add n 1)))
X (let (res)
X (while (not (car (setq res (math-prime-test
X! (setq n (math-add n -2)) iters)))))
X (if (and calc-verbose-nextprime
X (eq (car res) 'maybe))
X (calc-report-prime-test res)))
X--- 16195,16202 ----
X (setq n (math-add n 1)))
X (let (res)
X (while (not (car (setq res (math-prime-test
X! (setq n (math-add n -2))
X! (or iters 1))))))
X (if (and calc-verbose-nextprime
X (eq (car res) 'maybe))
X (calc-report-prime-test res)))
X***************
X*** 11018,11025 ****
X (setq x (cons (car x)
X (mapcar 'math-evaluate-expr-rec (cdr x)))))
X (if (eq (car-safe x) 'var)
X! (if (and (boundp (nth 2 x))
X! (symbol-value (nth 2 x))
X (not (eq (car-safe (symbol-value (nth 2 x)))
X 'incomplete)))
X (let ((val (symbol-value (nth 2 x))))
X--- 16612,16618 ----
X (setq x (cons (car x)
X (mapcar 'math-evaluate-expr-rec (cdr x)))))
X (if (eq (car-safe x) 'var)
X! (if (and (calc-var-value (nth 2 x))
X (not (eq (car-safe (symbol-value (nth 2 x)))
X 'incomplete)))
X (let ((val (symbol-value (nth 2 x))))
X***************
X*** 11037,11042 ****
X--- 16630,16637 ----
X )
X
X
X+ ;;;; [calc-arith.el]
X+
X ;;; Combine two terms being added, if possible.
X (defun math-combine-sum (a b nega negb scalar-okay)
X (if (and scalar-okay (Math-objvecp a) (Math-objvecp b))
X***************
X*** 11143,11148 ****
X--- 16738,16745 ----
X )
X
X
X+ ;;;; [calc-alg.el]
X+
X
X ;;; True if A comes before B in a canonical ordering of expressions. [P X X]
X (defun math-beforep (a b) ; [Public]
X***************
X*** 11172,11180 ****
X )
X
X
X-
X- ;;;; [calc-alg.el]
X-
X (setq math-living-dangerously nil) ; true if unsafe simplifications are okay.
X
X (defun math-simplify-extended (a)
X--- 16769,16774 ----
X***************
X*** 11181,11186 ****
X--- 16775,16781 ----
X (let ((math-living-dangerously t))
X (math-simplify a))
X )
X+ (fset 'calcFunc-esimplify (symbol-function 'math-simplify-extended))
X
X (defun math-simplify (top-expr)
X (calc-with-default-simplification
X***************
X*** 11191,11196 ****
X--- 16786,16792 ----
X (setq top-expr res))))
X top-expr
X )
X+ (fset 'calcFunc-simplify (symbol-function 'math-simplify))
X
X ;;; The following has a "bug" in that if any recursive simplifications
X ;;; occur only the first handler will be tried; this doesn't really
X***************
X*** 11309,11356 ****
X
X (defun math-simplify-divide ()
X (let ((np (cdr expr))
X! n nn)
X! (setq nn (math-common-constant-factor (nth 2 expr)))
X (if nn
X (progn
X (setq n (math-common-constant-factor (nth 1 expr)))
X! (if (and (consp nn) (eq (nth 1 nn) 1) (not n))
X (progn
X! (setcar (cdr expr) (math-mul (nth 1 expr) nn))
X (setcar (cdr (cdr expr))
X! (math-cancel-common-factor (nth 2 expr) nn)))
X (if (and n (not (eq (setq n (math-frac-gcd n nn)) 1)))
X (progn
X (setcar (cdr expr)
X (math-cancel-common-factor (nth 1 expr) n))
X (setcar (cdr (cdr expr))
X! (math-cancel-common-factor (nth 2 expr) n)))))))
X (while (eq (car-safe (setq n (car np))) '*)
X! (math-simplify-divisor (cdr n) (cdr (cdr expr)))
X (setq np (cdr (cdr n))))
X! (math-simplify-divisor np (cdr (cdr expr)))
X expr)
X )
X
X! (defun math-simplify-divisor (np dp)
X! (let ((n (car np))
X! d dd temp)
X! (while (eq (car-safe (setq d (car dp))) '*)
X! (if (setq temp (math-combine-prod n (nth 1 d) nil t t))
X! (progn
X! (setcar np (setq n temp))
X! (setcar (cdr d) 1)))
X! (setq dp (cdr (cdr d))))
X! (if (setq temp (math-combine-prod n d nil t t))
X! (progn
X! (setcar np (setq n temp))
X! (setcar dp 1))))
X )
X
X (defun math-common-constant-factor (expr)
X (if (Math-primp expr)
X (if (Math-ratp expr)
X! (and (not (memq expr '(0 1)))
X (math-abs expr))
X (if (Math-ratp (setq expr (math-to-simple-fraction expr)))
X (math-common-constant-factor expr)))
X--- 16905,16981 ----
X
X (defun math-simplify-divide ()
X (let ((np (cdr expr))
X! (nover nil)
X! (nn (math-common-constant-factor (nth 2 expr)))
X! n op)
X (if nn
X (progn
X (setq n (math-common-constant-factor (nth 1 expr)))
X! (if (and (eq (car-safe nn) 'frac) (eq (nth 1 nn) 1) (not n))
X (progn
X! (setcar (cdr expr) (math-mul (nth 2 nn) (nth 1 expr)))
X (setcar (cdr (cdr expr))
X! (math-cancel-common-factor (nth 2 expr) nn))
X! (if (and (math-negp nn)
X! (setq op (assq (car expr) calc-tweak-eqn-table)))
X! (setcar expr (nth 1 op))))
X (if (and n (not (eq (setq n (math-frac-gcd n nn)) 1)))
X (progn
X (setcar (cdr expr)
X (math-cancel-common-factor (nth 1 expr) n))
X (setcar (cdr (cdr expr))
X! (math-cancel-common-factor (nth 2 expr) n))
X! (if (and (math-negp n)
X! (setq op (assq (car expr) calc-tweak-eqn-table)))
X! (setcar expr (nth 1 op))))))))
X! (if (eq (car-safe (car np)) '/)
X! (progn
X! (setq np (cdr (nth 1 expr)))
X! (while (eq (car-safe (setq n (car np))) '*)
X! (math-simplify-divisor (cdr n) (cdr (cdr expr)) nil t)
X! (setq np (cdr (cdr n))))
X! (math-simplify-divisor np (cdr (cdr expr)) nil t)
X! (setq nover t
X! np (cdr (cdr (nth 1 expr))))))
X (while (eq (car-safe (setq n (car np))) '*)
X! (math-simplify-divisor (cdr n) (cdr (cdr expr)) nover t)
X (setq np (cdr (cdr n))))
X! (math-simplify-divisor np (cdr (cdr expr)) nover t)
X expr)
X )
X
X! (defun math-simplify-divisor (np dp nover dover)
X! (cond ((eq (car-safe (car dp)) '/)
X! (math-simplify-divisor np (cdr (car dp)) nover dover)
X! (math-simplify-divisor np (cdr (cdr (car dp))) nover (not dover)))
X! ((or (or (eq (car expr) '/)
X! (and (memq (car expr) '(calcFunc-eq calcFunc-neq))
X! math-living-dangerously))
X! (math-numberp (car np)))
X! (let ((n (car np))
X! d dd temp op)
X! (while (eq (car-safe (setq d (car dp))) '*)
X! (if (setq temp (math-combine-prod n (nth 1 d) nover dover t))
X! (progn
X! (and (math-negp (nth 1 d))
X! (setq op (assq (car expr) calc-tweak-eqn-table))
X! (setcar expr (nth 1 op)))
X! (setcar np (setq n (if nover (math-div 1 temp) temp)))
X! (setcar (cdr d) 1)))
X! (setq dp (cdr (cdr d))))
X! (if (setq temp (math-combine-prod n d nover dover t))
X! (progn
X! (and (math-negp (car dp))
X! (setq op (assq (car expr) calc-tweak-eqn-table))
X! (setcar expr (nth 1 op)))
X! (setcar np (setq n (if nover (math-div 1 temp) temp)))
X! (setcar dp 1))))))
X )
X
X (defun math-common-constant-factor (expr)
X (if (Math-primp expr)
X (if (Math-ratp expr)
X! (and (not (memq expr '(0 1 -1)))
X (math-abs expr))
X (if (Math-ratp (setq expr (math-to-simple-fraction expr)))
X (math-common-constant-factor expr)))
X***************
X*** 11374,11388 ****
X )
X
X (defun math-frac-gcd (a b)
X! (if (and (Math-integerp a)
X! (Math-integerp b))
X (math-gcd a b)
X! (or (Math-integerp a) (setq a (list 'frac a 1)))
X! (or (Math-integerp b) (setq b (list 'frac b 1)))
X (math-make-frac (math-gcd (nth 1 a) (nth 1 b))
X (math-gcd (nth 2 a) (nth 2 b))))
X )
X
X (math-defsimplify calcFunc-sin
X (or (and (eq (car-safe (nth 1 expr)) 'calcFunc-arcsin)
X (nth 1 (nth 1 expr)))
X--- 16999,17055 ----
X )
X
X (defun math-frac-gcd (a b)
X! (if (or (and (Math-integerp a)
X! (Math-integerp b))
X! (Math-zerop a)
X! (Math-zerop b))
X (math-gcd a b)
X! (and (Math-integerp a) (setq a (list 'frac a 1)))
X! (and (Math-integerp b) (setq b (list 'frac b 1)))
X (math-make-frac (math-gcd (nth 1 a) (nth 1 b))
X (math-gcd (nth 2 a) (nth 2 b))))
X )
X
X+ (math-defsimplify (calcFunc-eq calcFunc-neq calcFunc-lt
X+ calcFunc-gt calcFunc-leq calcFunc-geq)
X+ (math-simplify-ineq))
X+
X+ (defun math-simplify-ineq ()
X+ (math-simplify-divide)
X+ (let ((np (cdr expr)))
X+ (while (memq (car-safe (setq n (car np))) '(+ -))
X+ (math-simplify-add-term (cdr (cdr n)) (cdr (cdr expr)) (eq (car n) '-))
X+ (setq np (cdr n)))
X+ (math-simplify-add-term np (cdr (cdr expr)) nil)
X+ expr)
X+ )
X+
X+ (defun math-simplify-add-term (np dp minus)
X+ (let ((n (car np))
X+ d dd temp)
X+ (while (memq (car-safe (setq d (car dp))) '(+ -))
X+ (if (setq temp (math-combine-sum n (nth 2 d)
X+ minus (eq (car d) '+) t))
X+ (if (eq (math-looks-negp temp) minus)
X+ (progn
X+ (setcar np (setq n (if minus (math-neg temp) temp)))
X+ (setcar (cdr (cdr d)) 0))
X+ (progn
X+ (setcar np 0)
X+ (setcar (cdr (cdr d)) (setq n (if (eq (car d) '+)
X+ (math-neg temp)
X+ temp))))))
X+ (setq dp (cdr d)))
X+ (if (setq temp (math-combine-sum n d minus t t))
X+ (if (eq (math-looks-negp temp) minus)
X+ (progn
X+ (setcar np (setq n (if minus (math-neg temp) temp)))
X+ (setcar dp 0))
X+ (progn
X+ (setcar np 0)
X+ (setcar dp (setq n (math-neg temp)))))))
X+ )
X+
X (math-defsimplify calcFunc-sin
X (or (and (eq (car-safe (nth 1 expr)) 'calcFunc-arcsin)
X (nth 1 (nth 1 expr)))
X***************
X*** 11453,11459 ****
X (nth 1 (nth 1 expr)))
X (and math-living-dangerously
X (eq (car-safe (nth 1 expr)) 'calcFunc-cos)
X! (math-sub (math-div '(var pi var-pi) 2)
X (nth 1 (nth 1 expr)))))
X )
X
X--- 17120,17126 ----
X (nth 1 (nth 1 expr)))
X (and math-living-dangerously
X (eq (car-safe (nth 1 expr)) 'calcFunc-cos)
X! (math-sub (math-quarter-circle t)
X (nth 1 (nth 1 expr)))))
X )
X
X***************
X*** 11463,11469 ****
X (nth 1 (nth 1 expr)))
X (and math-living-dangerously
X (eq (car-safe (nth 1 expr)) 'calcFunc-sin)
X! (math-sub (math-div '(var pi var-pi) 2)
X (nth 1 (nth 1 expr)))))
X )
X
X--- 17130,17136 ----
X (nth 1 (nth 1 expr)))
X (and math-living-dangerously
X (eq (car-safe (nth 1 expr)) 'calcFunc-sin)
X! (math-sub (math-quarter-circle t)
X (nth 1 (nth 1 expr)))))
X )
X
X***************
X*** 11518,11535 ****
X (list '^ (nth 1 (nth 1 expr)) (math-div 1 4))))))
X )
X
X! (math-defsimplify 'calcFunc-exp
X (and (eq (car-safe (nth 1 expr)) 'calcFunc-ln)
X (nth 1 (nth 1 expr)))
X )
X
X! (math-defsimplify 'calcFunc-ln
X (and math-living-dangerously
X (eq (car-safe (nth 1 expr)) 'calcFunc-exp)
X (nth 1 (nth 1 expr)))
X )
X
X! (math-defsimplify '^
X (math-simplify-pow))
X
X (defun math-simplify-pow ()
X--- 17185,17202 ----
X (list '^ (nth 1 (nth 1 expr)) (math-div 1 4))))))
X )
X
X! (math-defsimplify calcFunc-exp
X (and (eq (car-safe (nth 1 expr)) 'calcFunc-ln)
X (nth 1 (nth 1 expr)))
X )
X
X! (math-defsimplify calcFunc-ln
X (and math-living-dangerously
X (eq (car-safe (nth 1 expr)) 'calcFunc-exp)
X (nth 1 (nth 1 expr)))
X )
X
X! (math-defsimplify ^
X (math-simplify-pow))
X
X (defun math-simplify-pow ()
X***************
X*** 11550,11556 ****
X (nth 1 (nth 2 expr))))
X )
X
X! (math-defsimplify 'calcFunc-log10
X (and math-living-dangerously
X (eq (car-safe (nth 1 expr)) '^)
X (math-equal-int (nth 1 (nth 1 expr)) 10)
X--- 17217,17223 ----
X (nth 1 (nth 2 expr))))
X )
X
X! (math-defsimplify calcFunc-log10
X (and math-living-dangerously
X (eq (car-safe (nth 1 expr)) '^)
X (math-equal-int (nth 1 (nth 1 expr)) 10)
X***************
X*** 11622,11714 ****
X
X
X
X! (defun math-apply-rewrite (expr lhs rhs &optional cond)
X! (let ((matches-found nil))
X! (and (math-match-pattern expr lhs)
X! (or (null cond)
X! (math-is-true (math-simplify (math-replace-variables cond))))
X! (math-replace-variables rhs)))
X! )
X!
X! (defun math-apply-rewrite-rules (expr rules)
X! (let ((r rules)
X! next)
X! (while (and r
X! (or (not (setq next (math-apply-rewrite expr
X! (nth 1 (car r))
X! (nth 2 (car r))
X! (nth 3 (car r)))))
X! (equal expr (setq next (math-normalize next)))))
X! (setq r (cdr r)))
X! (and r
X! next))
X! )
X
X (defun math-rewrite (expr rules &optional many)
X! (setq rules (math-check-rewrite-rules rules))
X! (math-map-tree (function (lambda (x) (math-apply-rewrite-rules x rules)))
X! expr many)
X )
X
X! (defun math-check-rewrite-rules (rules)
X (if (and (eq (car-safe rules) 'var)
X! (boundp (nth 2 rules))
X! (symbol-value (nth 2 rules)))
X! (setq rules (symbol-value (nth 2 rules))))
X! (or (Math-vectorp rules)
X! (error "Rules must be a vector"))
X! (setq rules (if (Math-vectorp (nth 1 rules))
X! (cdr rules)
X! (list rules)))
X! (let ((r rules))
X! (while r
X! (or (and (Math-vectorp (car r))
X! (cdr (cdr (car r)))
X! (not (nthcdr 4 (car r))))
X! (error "Malformed rules vector"))
X! (setq r (cdr r))))
X! rules
X! )
X!
X! (defun math-match-pattern (expr pat)
X! (cond ((Math-primp pat)
X! (or (math-equal expr pat)
X! (and (eq (car-safe pat) 'var)
X! (let ((match (assq (nth 1 pat) matches-found)))
X! (if match
X! (equal expr (nth 1 match))
X! (setq matches-found (cons (list (nth 1 pat)
X! expr)
X! matches-found)))))))
X! ((eq (car pat) 'calcFunc-quote)
X! (equal expr (nth 1 pat)))
X! (t
X! (and (eq (car pat) (car-safe expr))
X! (progn
X! (while (and (setq expr (cdr expr) pat (cdr pat))
X! expr
X! (math-match-pattern (car expr) (car pat))))
X! (and (null expr) (null pat))))))
X )
X
X! (defun math-replace-variables (expr)
X (if (Math-primp expr)
X (if (eq (car-safe expr) 'var)
X! (let ((match (assq (nth 1 expr) matches-found)))
X! (if match
X! (nth 1 match)
X expr))
X expr)
X! (cons (car expr) (mapcar 'math-replace-variables (cdr expr))))
X )
X
X ;;;; [calc-ext.el]
X
X (defun math-is-true (expr)
X (and (Math-realp expr)
X (not (Math-zerop expr)))
X )
X
X
X
X
X--- 17289,18634 ----
X
X
X
X! ;;;; [calc-rewr.el]
X!
X
X (defun math-rewrite (expr rules &optional many)
X! (let ((crules (math-compile-rewrites rules))
X! (heads (math-rewrite-heads expr))
X! result)
X! (math-map-tree (function
X! (lambda (x)
X! (if (setq result (math-apply-rewrites x crules heads))
X! (setq heads (math-rewrite-heads result heads)))
X! result))
X! expr many))
X! )
X!
X! (defun calcFunc-rewrite (expr rules &optional many)
X! (or (null many) (integerp many) (math-reject-arg 'integerp many))
X! (condition-case err
X! (math-rewrite expr rules (or many 1))
X! (error (math-reject-arg rules (nth 1 err))))
X! )
X!
X! (defun calcFunc-match (pat vec)
X! (or (math-vectorp vec) (math-reject-arg vec 'vectorp))
X! (condition-case err
X! (math-match-patterns pat vec nil)
X! (error (math-reject-arg pat (nth 1 err))))
X! )
X!
X! (defun calcFunc-matchnot (pat vec)
X! (or (math-vectorp vec) (math-reject-arg vec 'vectorp))
X! (condition-case err
X! (math-match-patterns pat vec t)
X! (error (math-reject-arg pat (nth 1 err))))
X! )
X!
X! (defun math-match-patterns (pat vec &optional not-flag)
X! (let ((newvec nil)
X! (crules (math-compile-patterns pat)))
X! (while (setq vec (cdr vec))
X! (if (eq (not (math-apply-rewrites (car vec) crules))
X! not-flag)
X! (setq newvec (cons (car vec) newvec))))
X! (cons 'vec (nreverse newvec)))
X! )
X!
X!
X!
X! ;;; A compiled rule set is an a-list of entries whose cars are functors,
X! ;;; and whose cdrs are lists of rules. If there are rules with no
X! ;;; well-defined head functor, they are included on all lists and also
X! ;;; on an extra list whose car is nil.
X! ;;;
X! ;;; Rule list entries take the form (regs prog head), where:
X! ;;;
X! ;;; regs is a vector of match registers.
X! ;;;
X! ;;; prog is a match program (see below).
X! ;;;
X! ;;; head is a rare function name appearing in the rule body (but not the
X! ;;; head of the whole rule), or nil if none.
X! ;;;
X! ;;; A match program is a list of match instructions.
X! ;;;
X! ;;; In the following, "part" is a register number that contains the
X! ;;; subexpression to be operated on.
X! ;;;
X! ;;; Register 0 is the whole expression being matched. The others are
X! ;;; meta-variables in the pattern, temporaries used for matching and
X! ;;; backtracking, and constant expressions.
X! ;;;
X! ;;; (same part reg)
X! ;;; The selected part must be math-equal to the contents of "reg".
X! ;;;
X! ;;; (same-neg part reg)
X! ;;; The selected part must be math-equal to the negative of "reg".
X! ;;;
X! ;;; (integer part)
X! ;;; The selected part must be an integer.
X! ;;;
X! ;;; (real part)
X! ;;; The selected part must be a real.
X! ;;;
X! ;;; (constant part)
X! ;;; The selected part must be a constant.
X! ;;;
X! ;;; (negative part)
X! ;;; The selected part must "look" negative.
X! ;;;
X! ;;; (rel part op reg)
X! ;;; The selected part must satisfy "part op reg", where "op"
X! ;;; is one of the 6 relational ops, and "reg" is a register.
X! ;;;
X! ;;; (mod part modulo value)
X! ;;; The selected part must satisfy "part % modulo = value", where
X! ;;; "modulo" and "value" are constants.
X! ;;;
X! ;;; (func part head reg1 reg2 ... regn)
X! ;;; The selected part must be an n-ary call to function "head".
X! ;;; The arguments are stored in "reg1" through "regn".
X! ;;;
X! ;;; (func-def part head defs reg1 reg2 ... regn)
X! ;;; The selected part must be an n-ary call to function "head".
X! ;;; "Defs" is a list of value/register number pairs for default args.
X! ;;; If a match, assign default values to registers and then skip
X! ;;; immediately over any following "func-def" instructions and
X! ;;; the following "func" instruction. If wrong number of arguments,
X! ;;; proceed to the following "func-def" or "func" instruction.
X! ;;;
X! ;;; (func-opt part head defs reg1)
X! ;;; Like func-def with "n=1", except that if the selected part is
X! ;;; not a call to "head", then the part itself successfully matches
X! ;;; "reg1" (and the defaults are assigned).
X! ;;;
X! ;;; (try part heads mark reg1 [def])
X! ;;; The selected part must be a function of the correct type which is
X! ;;; associative and/or commutative. "Heads" is a list of acceptable
X! ;;; types. An initial assignment of arguments to "reg1" is tried.
X! ;;; If the program later fails, it backtracks to this instruction
X! ;;; and tries other assignments of arguments to "reg1".
X! ;;; If "def" exists and normal matching fails, backtrack and assign
X! ;;; "part" to "reg1", and "def" to "reg2" in the following "try2".
X! ;;; The "mark" is a vector of size 4; only "mark[3]" is initialized.
X! ;;; "mark[0]" points to the argument list; "mark[1]" points to the
X! ;;; current argument; "mark[2]" is 0 if there are two arguments,
X! ;;; 1 if reg1 is matching single arguments, 2 if reg2 is matching
X! ;;; single arguments (a+b+c+d is never split as (a+b)+(c+d)), or
X! ;;; 3 if reg2 is matching "def"; "mark[3]" is 0 if the function must
X! ;;; have two arguments, 1 if phase-2 can be skipped, 2 if full
X! ;;; backtracking is necessary.
X! ;;;
X! ;;; (try2 try reg2)
X! ;;; Every "try" will be followed by a "try2" whose "try" field is
X! ;;; a pointer to the corresponding "try". The arguments which were
X! ;;; not stored in "reg1" by that "try" are now stored in "reg2".
X! ;;;
X! ;;; (apply part reg1 reg2)
X! ;;; The selected part must be a function call. The functor
X! ;;; (as a variable name) is stored in "reg1"; the arguments
X! ;;; (as a vector) are stored in "reg2".
X! ;;;
X! ;;; (cons part reg1 reg2)
X! ;;; The selected part must be a vector. The first element of
X! ;;; the vector is stored in "reg1"; the rest of the vector
X! ;;; (as another vector) is stored in "reg2".
X! ;;;
X! ;;; (select part reg)
X! ;;; If the selected part is a unary call to function "select", its
X! ;;; argument is stored in "reg"; otherwise (provided this is an `a r'
X! ;;; and not a `g r' command) the selected part is stored in "reg".
X! ;;;
X! ;;; (cond expr)
X! ;;; The "expr", with registers substituted, must simplify to
X! ;;; a non-zero value.
X! ;;;
X! ;;; (done rhs)
X! ;;; Rewrite the expression to "rhs", with register substituted.
X! ;;; Normalize; if the result is different from the original
X! ;;; expression, the match has succeeded. This is the last
X! ;;; instruction of every program.
X! ;;;
X!
X! ;;; Pseudo-functions related to rewrites:
X! ;;; In patterns: quote, plain, condition, opt, apply, cons, select
X! ;;; In righthand sides: quote, plain, eval, evalsimp, apply, cons, select
X!
X! ;;; Some optimizations that would be nice to have:
X! ;;;
X! ;;; * Merge registers with disjoint lifetimes.
X! ;;; * Merge constant registers with equivalent values.
X! ;;;
X! ;;; * If an argument of a commutative op math-depends neither on the
X! ;;; rest of the pattern nor on any of the conditions, then no backtracking
X! ;;; should be done for that argument. (This won't apply to very many
X! ;;; cases.)
X! ;;;
X! ;;; * If top functor is "select", and its argument is a unique function,
X! ;;; add the rule to the lists for both "select" and that function.
X! ;;; (Currently rules like this go on the "nil" list.)
X! ;;; Same for "func-opt" functions. (Though not urgent for these.)
X! ;;;
X!
X! ;;; Some additional features to add / things to think about:
X! ;;;
X! ;;; * Figure out what happens to "a +/- b" and "a +/- opt(b)".
X! ;;;
X! ;;; * Same for interval forms.
X! ;;;
X! ;;; * Have a name(v,pat) pattern which matches pat, and gives the
X! ;;; whole match the name v. Beware of circular structures!
X! ;;;
X!
X! (defun math-compile-patterns (pats)
X! (if (and (eq (car-safe pats) 'var)
X! (calc-var-value (nth 2 pats)))
X! (let ((prop (get (nth 2 pats) 'math-pattern-cache)))
X! (or prop
X! (put (nth 2 pats) 'math-pattern-cache (setq prop (list nil))))
X! (or (eq (car prop) (symbol-value (nth 2 pats)))
X! (progn
X! (setcdr prop (math-compile-patterns
X! (symbol-value (nth 2 pats))))
X! (setcar prop (symbol-value (nth 2 pats)))))
X! (cdr prop))
X! (let ((math-rewrite-whole t))
X! (math-compile-rewrites (cons
X! 'vec
X! (mapcar (function (lambda (x)
X! (list 'vec x
X! '(var XXX XXX))))
X! (if (eq (car-safe pats) 'vec)
X! (cdr pats)
X! (list pats)))))))
X )
X+ (setq math-rewrite-whole nil)
X+ (setq math-make-import-list nil)
X
X! (defun math-compile-rewrites (rules)
X (if (and (eq (car-safe rules) 'var)
X! (calc-var-value (nth 2 rules)))
X! (let ((prop (get (nth 2 rules) 'math-rewrite-cache))
X! (math-import-list nil)
X! (math-make-import-list t)
X! p)
X! (or prop
X! (put (nth 2 rules) 'math-rewrite-cache
X! (setq prop (list (list (cons (nth 2 rules) nil))))))
X! (setq p (car prop))
X! (while (and p (eq (symbol-value (car (car p))) (cdr (car p))))
X! (setq p (cdr p)))
X! (or (null p)
X! (progn
X! (message "Compiling rule set %s..." (nth 1 rules))
X! (setcdr prop (math-compile-rewrites
X! (symbol-value (nth 2 rules))))
X! (message "Compiling rule set %s...done" (nth 1 rules))
X! (setcar prop (cons (cons (nth 2 rules)
X! (symbol-value (nth 2 rules)))
X! math-import-list))))
X! (cdr prop))
X! (or (eq (car-safe rules) 'vec)
X! (error "Rewrite rules must be a vector"))
X! (if (assq 'calcFunc-import rules)
X! (let ((pp (setq rules (copy-sequence rules)))
X! p part)
X! (while (setq p (car (cdr pp)))
X! (if (eq (car-safe p) 'calcFunc-import)
X! (progn
X! (setcdr pp (cdr (cdr pp)))
X! (or (and (eq (car-safe (nth 1 p)) 'var)
X! (setq part (calc-var-value (nth 2 (nth 1 p))))
X! (eq (car-safe part) 'vec))
X! (error "Argument of import() must be a rules variable"))
X! (if math-make-import-list
X! (setq math-import-list
X! (cons (cons (nth 2 (nth 1 p))
X! (symbol-value (nth 2 (nth 1 p))))
X! math-import-list)))
X! (while (setq p (cdr (cdr p)))
X! (or (cdr p)
X! (error "import() must have odd number of arguments"))
X! (setq part (math-rwcomp-substitute part
X! (car p) (nth 1 p))))
X! (if (memq (car-safe (nth 1 part)) '(vec calcFunc-import))
X! (setq part (cdr part))
X! (setq part (list part)))
X! (setcdr pp (append part (cdr pp))))
X! (setq pp (cdr pp))))))
X! (if (eq (car-safe (nth 1 rules)) 'vec)
X! (setq rules (cdr rules))
X! (setq rules (list rules)))
X! (let ((rule-set nil)
X! (all-heads nil)
X! (nil-rules nil)
X! (rule-count 0))
X! (while rules
X! ;; (message "Compiling rule %d..." (setq rule-count (1+ rule-count)))
X! (or (and (eq (car-safe (car rules)) 'vec)
X! (cdr (cdr (car rules)))
X! (not (nthcdr 4 (car rules))))
X! (error "Rewrite rule set must be a vector of vectors"))
X! (let ((math-pattern (nth 1 (car rules)))
X! (math-rhs (nth 2 (car rules)))
X! (math-prog nil)
X! (math-num-regs 1)
X! (math-regs (list (list nil 0 nil nil)))
X! (math-conds nil))
X! (let ((cond (nth 3 (car rules))))
X! (if (and (eq (car-safe cond) 'calcFunc-quote)
X! (= (length cond) 2))
X! (setq cond (nth 1 cond)))
X! (while (eq (car-safe cond) 'calcFunc-land)
X! (setq math-conds (cons (nth 2 cond) math-conds)
X! cond (nth 1 cond)))
X! (if cond
X! (setq math-conds (cons cond math-conds))))
X! (math-rwcomp-pattern math-pattern 0)
X! (while math-conds
X! (math-rwcomp-cond-instr (car math-conds))
X! (setq math-conds (cdr math-conds)))
X! (math-rwcomp-instr 'done (math-rwcomp-match-vars math-rhs))
X! (setq math-prog (nreverse math-prog))
X! (let* ((heads (math-rewrite-heads math-pattern))
X! (rule (list (vconcat
X! (nreverse
X! (mapcar (function (lambda (x) (nth 3 x)))
X! math-regs)))
X! math-prog
X! heads))
X! (head (and (not (Math-primp math-pattern))
X! (not (and (eq (car (car math-prog)) 'try)
X! (nth 5 (car math-prog))))
X! (not (memq (car (car math-prog)) '(func-opt
X! apply
X! select)))
X! (if (memq (car (car math-prog)) '(func
X! func-def))
X! (nth 2 (car math-prog))
X! (car math-pattern)))))
X! (let (found)
X! (while heads
X! (if (setq found (assq (car heads) all-heads))
X! (setcdr found (1+ (cdr found)))
X! (setq all-heads (cons (cons (car heads) 1) all-heads)))
X! (setq heads (cdr heads))))
X! (if (eq head '-) (setq head '+))
X! (if (eq head 'calcFunc-cons) (setq head 'vec))
X! (if head
X! (nconc (or (assq head rule-set)
X! (car (setq rule-set (cons (cons head
X! (copy-sequence
X! nil-rules))
X! rule-set))))
X! (list rule))
X! (setq nil-rules (nconc nil-rules (list rule)))
X! (let ((ptr rule-set))
X! (while ptr
X! (nconc (car ptr) (list rule))
X! (setq ptr (cdr ptr))))))
X! (setq rules (cdr rules))))
X! ;; (message "Collating rule heads...")
X! (if nil-rules
X! (setq rule-set (cons (cons nil nil-rules) rule-set)))
X! (setq all-heads (mapcar 'car
X! (sort all-heads (function
X! (lambda (x y)
X! (< (cdr x) (cdr y)))))))
X! (let ((set rule-set)
X! rule heads ptr)
X! (while set
X! (setq rule (cdr (car set)))
X! (while rule
X! (if (consp (setq heads (nth 2 (car rule))))
X! (progn
X! (setq heads (delq (car (car set)) heads)
X! ptr all-heads)
X! (while (and ptr (not (memq (car ptr) heads)))
X! (setq ptr (cdr ptr)))
X! (setcar (nthcdr 2 (car rule)) (car ptr))))
X! (setq rule (cdr rule)))
X! (setq set (cdr set))))
X! (let ((plus (assq '+ rule-set)))
X! (if plus
X! (setq rule-set (cons (cons '- (cdr plus)) rule-set))))
X! rule-set))
X! )
X!
X! (defun math-rewrite-heads (expr &optional more)
X! (let ((heads more))
X! (or (Math-primp expr)
X! (math-rewrite-heads-rec expr))
X! heads)
X! )
X!
X! (defun math-rewrite-heads-rec (expr)
X! (or (memq (car expr) heads)
X! (memq (car expr) '(calcFunc-quote calcFunc-plain calcFunc-opt
X! calcFunc-select calcFunc-apply
X! calcFunc-cons calcFunc-condition))
X! (memq 'algebraic (get (car expr) 'math-rewrite-props))
X! (setq heads (cons (car expr) heads)))
X! (while (setq expr (cdr expr))
X! (or (Math-primp (car expr))
X! (math-rewrite-heads-rec (car expr))))
X )
X
X! (defun math-rwcomp-match-vars (expr)
X (if (Math-primp expr)
X (if (eq (car-safe expr) 'var)
X! (let ((entry (assq (nth 2 expr) math-regs)))
X! (if entry
X! (math-rwcomp-register-expr (nth 1 entry))
X expr))
X expr)
X! (if (and (eq (car expr) 'calcFunc-quote)
X! (= (length expr) 2))
X! (math-rwcomp-match-vars (nth 1 expr))
X! (cons (car expr)
X! (mapcar 'math-rwcomp-match-vars (cdr expr)))))
X! )
X!
X! (defun math-rwcomp-register-expr (num)
X! (let ((entry (nth (1- (- math-num-regs num)) math-regs)))
X! (if (nth 2 entry)
X! (list 'neg (list 'calcFunc-register (nth 1 entry)))
X! (list 'calcFunc-register (nth 1 entry))))
X! )
X!
X! (defun math-rwcomp-substitute (expr old new)
X! (if (and (eq (car-safe old) 'var)
X! (memq (car-safe new) '(var calcFunc-lambda)))
X! (let ((old-func (math-var-to-calcFunc old))
X! (new-func (math-var-to-calcFunc new)))
X! (math-rwcomp-subst-rec expr))
X! (let ((old-func nil))
X! (math-rwcomp-subst-rec expr)))
X! )
X!
X! (defun math-rwcomp-subst-rec (expr)
X! (cond ((equal expr old) new)
X! ((Math-primp expr) expr)
X! (t (if (eq (car expr) old-func)
X! (math-build-call new-func (mapcar 'math-rwcomp-subst-rec
X! (cdr expr)))
X! (cons (car expr)
X! (mapcar 'math-rwcomp-subst-rec (cdr expr))))))
X! )
X!
X! (setq math-rwcomp-tracing nil)
X! (defun math-rwcomp-trace (instr)
X! (if math-rwcomp-tracing (progn (terpri) (princ instr)))
X! instr
X! )
X!
X! (defun math-rwcomp-instr (&rest instr)
X! (setq math-prog (cons (math-rwcomp-trace instr) math-prog))
X! )
X!
X! (defun math-rwcomp-multi-instr (tail &rest instr)
X! (setq math-prog (cons (math-rwcomp-trace (append instr tail))
X! math-prog))
X! )
X!
X! (defun math-rwcomp-cond-instr (expr)
X! (setq expr (math-rwcomp-match-vars expr))
X! (let (op arg)
X! (cond ((and (setq op (cdr (assq (car-safe expr)
X! '( (calcFunc-integer . integer)
X! (calcFunc-real . real)
X! (calcFunc-constant . constant)
X! (calcFunc-negative . negative) ))))
X! (= (length expr) 2)
X! (or (and (eq (car-safe (nth 1 expr)) 'neg)
X! (memq op '(integer real constant))
X! (setq arg (nth 1 (nth 1 expr))))
X! (setq arg (nth 1 expr)))
X! (eq (car-safe (setq arg (nth 1 expr))) 'calcFunc-register))
X! (math-rwcomp-instr op (nth 1 arg)))
X! ((and (assq (car-safe expr) calc-tweak-eqn-table)
X! (= (length expr) 3)
X! (eq (car-safe (nth 1 expr)) 'calcFunc-register))
X! (if (math-constp (nth 2 expr))
X! (let ((reg (math-rwcomp-reg)))
X! (setcar (nthcdr 3 (car math-regs)) (nth 2 expr))
X! (math-rwcomp-instr 'rel (nth 1 (nth 1 expr))
X! (car expr) reg))
X! (if (eq (car (nth 2 expr)) 'calcFunc-register)
X! (math-rwcomp-instr 'rel (nth 1 (nth 1 expr))
X! (car expr) (nth 1 (nth 2 expr)))
X! (math-rwcomp-instr 'cond expr))))
X! ((and (eq (car-safe expr) 'calcFunc-eq)
X! (= (length expr) 3)
X! (eq (car-safe (nth 1 expr)) '%)
X! (eq (car-safe (nth 1 (nth 1 expr))) 'calcFunc-register)
X! (math-constp (nth 2 (nth 1 expr)))
X! (math-constp (nth 2 expr)))
X! (math-rwcomp-instr 'mod (nth 1 (nth 1 (nth 1 expr)))
X! (nth 2 (nth 1 expr)) (nth 2 expr)))
X! (t (math-rwcomp-instr 'cond expr))))
X! )
X!
X! (defun math-rwcomp-same-instr (reg1 reg2 neg)
X! (math-rwcomp-instr (if (eq (eq (nth 2 (math-rwcomp-reg-entry reg1))
X! (nth 2 (math-rwcomp-reg-entry reg2)))
X! neg)
X! 'same-neg
X! 'same)
X! reg1 reg2)
X! )
X!
X! (defun math-rwcomp-reg ()
X! (prog1
X! math-num-regs
X! (setq math-regs (cons (list nil math-num-regs nil nil) math-regs)
X! math-num-regs (1+ math-num-regs)))
X! )
X!
X! (defun math-rwcomp-reg-entry (num)
X! (nth (1- (- math-num-regs num)) math-regs)
X! )
X!
X! (defun math-rwcomp-pattern (expr part)
X! (cond ((or (math-rwcomp-no-vars expr)
X! (and (eq (car expr) 'calcFunc-quote)
SHAR_EOF
echo "End of part 9, continue with part 10"
echo "10" > s2_seq_.tmp
exit 0