home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
ARM Club 3
/
TheARMClub_PDCD3.iso
/
hensa
/
maths
/
b116_1
/
jacal
/
norm
< prev
next >
Wrap
Text File
|
1993-10-03
|
9KB
|
272 lines
;;; JACAL: Symbolic Mathematics System. -*-scheme-*-
;;; Copyright 1989, 1990, 1991, 1992, 1993 Aubrey Jaffer.
;;; See the file "COPYING" for terms applying to this program.
(proclaim '(optimize (speed 3) (compilation-speed 0)))
(define (vsubst new old e)
(cond ((eq? new old) e)
((number? e) e)
((bunch? e) (map (lambda (e) (vsubst new old e)) e))
((var:> new (car e)) (univ:norm0 new (cdr (promote old e))))
(else (poly:resultant (make-var-eqn new old) e old))))
(define (make-var-eqn new old)
(if (var:> old new)
(list old (list new 0 -1) 1)
(list new (list old 0 -1) 1)))
(define (swapvars x y p)
(vsubst x _$
(vsubst y x
(vsubst _$ y p))))
;;; Makes an expression whose value is the variable VAR in the equation
;;; E or (if E is an expression) E=0
(define (suchthat var e)
(set! e (poly:subst0 $ (licit->poleqn e)))
(extize (normalize (vsubst $ var e))))
;; canonicalizers
(define (normalize x)
(cond ((and math:phases (not (novalue? x)))
(display-diag 'normalizing:)
(newline-diag)
(math:write x *output-grammar*)))
(let ((ans (normalize1 x)))
(cond ((and math:phases (not (novalue? x)))
(display-diag 'yielding:)
(newline-diag)
(math:write ans *output-grammar*)))
ans))
(define (normalize1 x)
(cond ((bunch? x) (map normalize x))
((symbol? x) (eval-error 'normalize-symbol?- x))
((eqn? x)
(poly->eqn (signcan (poly:square-and-num-cont-free
(alg:simplify (eqn->poly x))))))
(else (expr:normalize x))))
(define (expr:normalize p)
(if (expl? p) (set! p (expl->impl p)))
(expr:norm-or-signcan
(poly:square-free-var
(alg:simplify (if (rat? p) (alg:clear-denoms p) p))
$)))
(define (extize p)
(cond ((bunch? p) (eval-error 'cannot-suchthat-a-vector p))
((eqn? p) p)
((expl? p) p)
((rat? p) p)
(else
(set! newextstr (sect:next-string newextstr))
(let ((v (defext (string->var
(if (clambda? p) (string-append "@" newextstr)
newextstr))
p)))
(set! var-news (cons v var-news))
(var->expl v)))))
;(trace normalize normalize1 extize signcan
; expr:norm-or-signcan expr:normalize
; alg:simplify alg:clear-denoms
; poly:square-free-var poly:square-and-num-cont-free)
;; differentials
(define (total-diffn p vars)
(if (null? vars) 0
(poly:+ (poly:* (var->expl (var:differential (car vars)))
(poly:diff p (car vars)))
(total-diffn p (cdr vars)))))
(define (chain-rule v vd)
(if (extrule v)
(total-chain-exts (total-diffn (extrule v) (poly:vars (extrule v)))
(var:funcs (extrule v)))
(let ((functor (sexp->math (car (var:sexp v)))))
(do ((pos 1 (+ 1 pos))
(al (cdr (func-arglist v)) (cdr al))
(sum 0 (app* $1*$2+$3
(apply deferop
(deferop _partial functor pos)
(cdr (func-arglist v)))
(total-differential (car al))
sum)))
((null? al) (vsubst vd $ sum))))))
(define (total-chain-exts drule es)
(if (null? es) drule
(let ((ed (var:differential (car es))))
(total-chain-exts
(poly:resultant drule (chain-rule (car es) ed) ed)
(if (extrule (car es))
(union (cdr es) (alg:exts (extrule (car es))))
(cdr es))))))
(define (total-differential a)
(cond
((bunch? a) (map total-differential a))
((eqn? a) (poly->eqn
(total-diffn (eqn->poly a) (poly:vars (eqn->poly a)))))
(else (let ((aes (chainables a)))
(if (and (null? aes) (expl? a))
(total-diffn a (poly:vars a))
(let ((pa (licit->poleqn a)))
(total-chain-exts
(vsubst $ d$ (poly:resultant
pa (total-diffn pa (poly:vars pa)) $))
aes)))))))
(define (diff a var)
(cond
((bunch? a) (map (lambda (x) (diff x var)) a))
((eqn? a) (poly->eqn (diff (eqn->poly a) var)))
(else (let* ((td (total-differential a))
(vd (var->expl (var:differential var)))
(td1 (app* $1/$2 td vd))
(dpvs '()))
(poly:for-each-var
(lambda (v)
(if (and (not (eq? (car vd) v))
(var:differential? v))
(set! dpvs (adjoin v dpvs))))
td)
(reduce-init (lambda (e x) (poly:coeff e x 0)) td1 dpvs)))))
(define (expls:diff a var)
(cond
((bunch? a) (map (lambda (x) (expl:diff x var)) a))
((eqn? a) (poly->eqn (expl:diff (eqn->poly a) var)))
(else (let ((aes (alg:exts a)))
(if (and (null? aes) (expl? a)
(every (lambda (x) (null? (var:depends x))) (poly:vars a)))
(poly:diff a var)
(math:error 'polydiff 'not-a-polynomial? a))))))
;(trace total-differential total-chain-exts chain-rule total-diffn
; diff expls:diff)
;;;; FINITE DIFFERENCES
;;; shift needs to go through extensions; which will create new
;;; extensions (yucc). It is clear what to do for radicals, but other
;;; extensions will be hard to link up. For instance y: {x|x^5+a+b+9+x}
;;; needs to yield the same number whether a or b is substituted first.
(define (shift p var)
(vsubst var
$2
(poly:resultant (list $2 (list var -1 -1) 1)
p
var)))
(define (unsum p var)
(app* $1-$2 p (shift p (expl->var var))))
(define (poly:fdiffn p vars)
(if (null? vars) 0
(poly:+ (poly:* (var->expl (var:finite-differential (car vars)))
(unsum p (car vars)))
(poly:fdiffn p (cdr vars)))))
(define (total-finite-differential e)
(if (bunch? e)
(map total-finite-differential e)
(poly:fdiffn e (alg:vars e))))
;;;logical operations on licits
;(define (impl:not p)
; (poly:+ (poly:* (licit->poleqn p)
; (var->expl (sexp->var (new-symbol "~")))) -1))
;(define (impl:and p . qs)
; (cond ((bunch? p) (impl:and (append p qs)))))
(define (expl:t? e) (equal? e expl:t))
(define (ncexpt a pow)
(cond ((not (or (integer? pow) (expl:t? pow)))
(deferop _^^ a pow))
((eqns? a) (math:error 'Expt-of-equation?:- a))
((not (bunch? a)) (fcexpt a pow))
((expl:t? pow) (transpose a))
(else (mtrx:expt a pow))))
;;;; Routines for factoring
(define (poly:diff-coefs el n)
(if (null? el)
el
(cons (poly:* n (car el))
(poly:diff-coefs (cdr el) (+ 1 n)))))
(define (poly:diff p var)
(cond ((number? p) 0)
((eq? (car p) var) (univ:norm0 var (poly:diff-coefs (cddr p) 1)))
((var:> var (car p)) 0)
(else (univ:norm0 (car p) (map-no-end-0s
(lambda (x) (poly:diff x var))
(cdr p))))))
(define (poly:diff-all p)
(let ((ans 0))
(do ((vars (poly:vars p) (cdr vars)))
((null? vars) ans)
(set! ans (poly:+ (poly:diff p (car vars)) ans)))))
;;;functions involved with square free.
(define (poly:square-free-var p var)
(poly:/ p (poly:gcd p (poly:diff p var))))
(define (poly:square-and-num-cont-free p)
(if (number? p) (if (zero? p) p 1)
(poly:* (poly:square-and-num-cont-free (univ:cont p))
(poly:/ p (poly:gcd p (poly:diff p (car p)))))))
(define (poly:factorq p) (poly:factor-all p))
(define (poly:factor-var c var)
(poly:factor-split c (poly:diff c var)))
(define (poly:factor-all c)
(poly:factor-split c (poly:diff-all c)))
;;; This algorithm is due to:
;;; Multivariate Polynomial Factorization
;;; David R. Musser
;;; Journal of the Association for Computing Machinery
;;; Vol 22, No. 2, April 1975
(define (poly:factor-split c splitter)
(let ((d '()) (aj '()) (b (poly:gcd c splitter))) ; changed #f's to ()
(do ((b b (poly:/ b d))
(a (poly:/ c b) d))
((number? b)
(if (eqv? 1 b)
(cons a aj) ; nreverse removed
(cons b (cons a aj)))) ; nreverse removed
(set! d (poly:gcd a b))
(set! a (poly:/ a d)) ; 'a' is used as a temporary variable
(if (not (eqv? a 1)) ; keeps extraneous `1's
(set! aj (cons a aj)))))) ; out of the results list
;;; the following algorithm attempts to separate factors in a multivariate
;;; polynomial with major variable. It substitues 0 for each variable
;;; that it finds in turn and takes GCD against the original expression.
;;; It assumes that it's argument is squarefree and contentfree in the
;;; major variable.
(define (univ:split pe varlist)
(cond ((unit? pe) (list))
((null? varlist) (list pe))
((let ((p0 (signcan
(poly:gcd pe (poly:subst0 (car varlist) pe))))
(cvl (cdr varlist)))
(if (unit? p0)
(univ:split pe cvl)
(nconc (univ:split (poly:/ pe p0) cvl)
(univ:split p0 cvl)))))))
(define (univ:split-all poly) (univ:split poly (poly:vars poly)))
(define (factor-free-var p var)
(poly:gcd p (poly:subst0 var p)))
(define (factor:test)
(test (list (list x -1 -2) (list x -1 0 1))
poly:factorq
(list x -1 -4 -3 4 4)))
;;; Copyright 1989, 1990, 1991, 1992, 1993 Aubrey Jaffer.