ARM Club 3/TheARMClub_PDCD3.iso/hensa/maths/b116_1/jacal/poly

Text File | 1993-06-11 | 16KB | 508 lines

;;; JACAL: Symbolic Mathematics System.        -*-scheme-*-
;;; Copyright 1989, 1990, 1991, 1992, 1993 Aubrey Jaffer.
;;; See the file "COPYING" for terms applying to this program.

;;;Functions which operate on polynomials as polynomials
;;;have prefix POLY:
;;;Functions which operate on polynomials with the same major variable
;;;have prefix UNIV:
;;;Functions which operate on coefficients (without major varible)
;;;have prefix COES:

;(proclaim '(optimize (speed 3) (compilation-speed 0)))

(define (one? n) (eqv? 1 n))
;;; POLY:0? should do those normalizations which could lead to colapse
;;; of an expression to 0.
(define (poly:0? p) (eqv? 0 p))
(define (divides? a b) (zero? (remainder b a)))

;;; poly is the internal workhorse data type in the form
;;; of numeric or list (var coeff0 coeff1 ...)
;;; where var is a variable and coeffn is the coefficient of var^n.
;;; coeffn is poly.  The variables currently are arranged
;;; with the reverse alphabetically with z higher order than A.

(define (poly:find-var? poly var)
  (poly:find-var-if? poly (lambda (x) (eqv? var x))))

(define (poly:find-var-if? poly proc)
  (cond ((number? poly) #f)
    ((proc (car poly)))
    (else (some (lambda (x) (poly:find-var-if? x proc)) (cdr poly)))))

;;; This can call proc more than once per var
(define (poly:for-each-var proc poly)
  (cond ((number? poly))
    (else
     (proc (car poly))
     (for-each (lambda (b) (poly:for-each-var proc b))
           (cdr poly)))))

;;;POLY:VARS returns a list of all vars used in POLY
(define (poly:vars poly)
  (let ((elts '()))
    (poly:for-each-var (lambda (v) (set! elts (adjoin v elts))) poly)
    elts))

;;;; the following functions are for internal use on the poly data type

;;; this normalizes short polys.
(define (univ:norm0 var col)
  (cond ((null? col) 0)
    ((null? (cdr col)) (car col))
    (else (cons var col))))

(define (map-no-end-0s proc l)
  (if (null? l)
      l
    (let ((res (proc (car l)))
      (rest (map-no-end-0s proc (cdr l))))
      (if (and (null? rest) (eqv? 0 res))
      rest
    (cons res rest)))))
(define (map2c-no-end-0s proc l1 l2)
  (cond ((null? l1) l2)
    ((null? l2) l1)
    (else
     (let ((res (proc (car l1) (car l2)))
      (rest (map2c-no-end-0s proc (cdr l1) (cdr l2))))
       (if (and (null? rest) (eqv? 0 res))
           rest
         (cons res rest))))))
(define (univ-w/o-lt p)
  (univ:norm0 (car p) (map-no-end-0s identity (butlast (cdr p) 1))))
(define (make-list-length lst len fill)
  (let ((ld (- len (length lst))))
    (cond ((negative? ld) (butlast lst (- ld)))
      (else (append lst (make-list ld fill))))))

(define (ipow-by-squaring x n acc proc)
  (cond ((zero? n) acc)
    ((one? n) (proc acc x))
    (else (ipow-by-squaring (proc x x)
                (quotient n 2)
                (if (even? n) acc (proc acc x))
                proc))))

(define (poly:add-const term p2)
  (cons (car p2) (cons (poly:+ term (cadr p2)) (cddr p2))))
(define (poly:+ p1 p2)
  (cond ((and (number? p1) (number? p2)) (+ p1 p2))
    ((number? p1) (if (zero? p1) p2 (poly:add-const p1 p2)))
    ((number? p2) (if (zero? p2) p1 (poly:add-const p2 p1)))
    ((eq? (car p1) (car p2))
     (univ:norm0 (car p1) (map2c-no-end-0s poly:+ (cdr p1) (cdr p2))))
    ((var:> (car p2) (car p1)) (poly:add-const p1 p2))
    (else (poly:add-const p2 p1))))

(define (univ:* p1 p2)
  (let ((res (make-list (+ (length (cdr p1)) (length (cdr p2)) -1) 0)))
    (do ((rpl res (cdr rpl))
     (a (cdr p1) (cdr a)))
    ((null? a) (cons (car p1) res))
    (do ((b (cdr p2) (cdr b))
         (rp rpl (cdr rp)))
        ((null? b))
        (set-car! rp (poly:+ (poly:* (car a) (car b)) (car rp)))))))
(define (poly:times-const term p2)
  (cons (car p2) (map (lambda (x) (poly:* term x)) (cdr p2))))
(define (poly:* p1 p2)
  (cond ((and (number? p1) (number? p2)) (* p1 p2))
    ((number? p1) (cond ((zero? p1) 0)
                ((one? p1) p2)
                (else (poly:times-const p1 p2))))
    ((number? p2) (cond ((zero? p2) 0)
                ((one? p2) p1)
                (else (poly:times-const p2 p1))))
    ((eq? (car p1) (car p2)) (univ:* p1 p2))
    ((var:> (car p2) (car p1)) (poly:times-const p1 p2))
    (else (poly:times-const p2 p1))))

(define (poly:negate p) (poly:* -1 p))
(define (poly:- p1 p2) (poly:+ p1 (poly:negate p2)))

(define (univ:/? u v)
  (let ((r (list->vector (cdr u)))
    (m (length (cddr u)))
    (n (length (cddr v)))
    (vn (car (last-pair v)))
    (q '()))
    (do ((k (- m n) (- k 1))
     (qk (poly:/? (vector-ref r m) vn)
         (and (> k 0) (poly:/? (vector-ref r (+ n k -1)) vn))))
    ((not qk)
     (and (< k 0)
          (do ((k (- n 2) (- k 1)))
          ((or (< k 0) (not (poly:0? (vector-ref r k))))
           (< k 0)))
          (univ:norm0 (car u) q)))
      (set! q (cons qk q))
      (let ((qk- (poly:negate qk)))
    (do ((j (+ n k -1) (- j 1)))
        ((< j k))
      (vector-set! r j (poly:+ (vector-ref r j)
                   (poly:* (list-ref v (+ (- j k) 1)) qk-))))))))

;;; POLY:/? returns U / V if V divides U, otherwise returns #f
(define (poly:/? u v)
  (cond ((equal? u v) 1)
    ((eqv? 0 u) 0)
    ((eqv? 0 v) #f)
    ((unit? v) (poly:* v u))
    ((and (number? u) (number? v)) (and (divides? v u) (quotient u v)))
    ((number? v) (univ:/? u (const:promote (car u) v)))
    ((number? u) #f)
    ((eq? (car u) (car v)) (univ:/? u v))
    ((var:> (car u) (car v))
     (univ:/? u (const:promote (car u) v)))
    (else #f)))

(define (univ:/ dividend divisor)
  (or (univ:/? dividend divisor)
      (math:error divisor 'does-not-udivide- dividend)))

(define (poly:/ dividend divisor)
  (or (poly:/? dividend divisor)
      (math:error divisor 'does-not-divide- dividend)))

(define (univ:monomial coeff n var)
  (cond ((eq? 0 coeff) 0)
    ((>= 0 n) coeff)
    (else
     (cons var
           ((lambda (x) (set-car! (last-pair x) coeff) x)
        (make-list (+ 1 n) 0))))))

(define (poly:degree p var)
   (cond ((number? p) 0)
     ((eq? var (car p)) (length (cddr p)))
     ((var:> var (car p)) 0)
     (else (reduce-init (lambda (m c) (max m (poly:degree c var)))
                0
                (cdr p)))))

;(define (leading-coeff p var) (poly:coeff p var (poly:degree p var)))
;(define (monic? u var) (one? (leading-coeff u var)))

(define (poly:^ x n)
  (if (number? x)
;      (expt x n)
    (ipow-by-squaring x n 1 *)
    (ipow-by-squaring x n 1 poly:*)))

;;;; Routines used in normalizing IMPL polynomials

(define (leading-number p)
  (if (number? p) p (leading-number (car (last-pair p)))))

;;; This canonicalizes polys with respect to sign by forcing the
;;; numerical coefficient of the a certain term to always be positive.
(define (signcan p)
  (if (negative? (leading-number p)) (poly:negate p) p))

(define (shorter? x y) (< (length x) (length y)))

(define (univ:degree p var)
  (if (or (number? p) (not (eq? (car p) var))) 0 (length (cddr p))))

;;; THE NEXT 3 ROUTINES FOR SUBRESULTANT GCD ASSUME THAT THE ARGUMENTS
;;; ARE POLYNOMIALS WITH THE SAME MAJOR VARIABLE.
;;;  THESE TWO ROUTINES ASSUME THAT THE FIRST ARGUMENT IS OF GREATER
;;;  OR EQUAL ORDER THAN THE SECOND.
;;; These algorithms taken from:
;;; Knuth, D. E.,
;;; The Art Of Computer Programming, Vol. 2: Seminumerical Algorithms,
;;; Addison Wesley, Reading, MA 1969.
;;; Pseudo Remainder

;;; This returns a list of the pseudo quotient and pseudo remainder.
(define (univ:pdiv u v)
  (let* ((r (list->vector (cdr u)))
     (m (length (cddr u)))
     (n (length (cddr v)))
     (vn (car (last-pair v)))
     (q (make-vector (+ (- m n) 1) 1)))
    (do ((tt (- (- m n) 1) (- tt 1))
     (k 1 (+ 1 k))
     (vnp 1))
    ((< tt 0))
    (set! vnp (poly:* vnp vn))
    (vector-set! q k vnp)
    (vector-set! r tt (poly:* (vector-ref r tt) vnp)))
    (do ((k (- m n) (- k 1))
     (rnk 0))
    ((< k 0))
      (set! rnk (poly:negate (vector-ref r (+ n k))))
      (do ((j (+ n k -1) (- j 1)))
      ((< j k))
    (vector-set! r j (poly:+ (poly:* (vector-ref r j) vn)
                 (poly:* (list-ref v (+ (- j k) 1)) rnk)))))
    (list
     (do ((k (- m n) (+ -1 k))
      (end '() (cons (poly:* (vector-ref r (+ n k))
                 (vector-ref q k)) end)))
     ((zero? k) (univ:norm0 (car u) (cons (vector-ref r n) end))))
     (do ((j (- n 1) (- j 1))
      (end '()))
     ((< j 0) (univ:norm0 (car u) end))
     (if (not (and (null? end) (eqv? 0 (vector-ref r j))))
         (set! end (cons (vector-ref r j) end)))))))
(define (poly:pdiv dividend divisor var)
  (let ((pd1 (poly:degree dividend var))
    (pd2 (poly:degree divisor var)))
    (cond ((< pd1 pd2) (list 0 dividend))
      ((zero? (+ pd1 pd2))
       (list (quotient dividend divisor) (remainder dividend divisor)))
      ((zero? pd1) (list 0 dividend))
      ((zero? pd2)
;;; This should work but doesn't.
;;;       (map univ:demote (univ:pdiv (promote var dividend)
;;;                       (const:promote var divisor)))
       (list 0 dividend))
      (else
       (map univ:demote (univ:pdiv (promote var dividend)
                       (promote var divisor)))))))
(define (univ:prem u v)
  (let* ((r (list->vector (cdr u)))
     (m (length (cddr u)))
     (n (length (cddr v)))
     (vn (car (last-pair v))))
    (do ((k (- (- m n) 1) (- k 1))
     (vnp 1))
    ((< k 0))
      (set! vnp (poly:* vnp vn))
      (vector-set! r k (poly:* (vector-ref r k) vnp)))
    (do ((k (- m n) (- k 1))
     (rnk 0))
    ((< k 0))
      (set! rnk (poly:negate (vector-ref r (+ n k))))
      (do ((j (+ n k -1) (- j 1)))
      ((< j k))
    (vector-set! r j (poly:+ (poly:* (vector-ref r j) vn)
                 (poly:* (list-ref v (+ (- j k) 1)) rnk)))))
    (do ((j (- n 1) (- j 1))
     (end '()))
    ((< j 0) (univ:norm0 (car u) end))
      (if (and (null? end) (eqv? 0 (vector-ref r j)))
      #f
      (set! end (cons (vector-ref r j) end))))))

;;; Pseudo Remainder Sequence
(define (univ:prs u v)
  (let ((var (car u))
    (g 1)
    (h 1)
    (delta 0))
    (do ((r (univ:prem u v) (univ:prem u v)))
    ((eqv? 0 (univ:degree r var))
     (if (eqv? 0 r) v r))
      (set! delta (- (univ:degree u var) (univ:degree v var)))
      (set! u v)
      (set! v (univ:/ r (const:promote (car r) (poly:* g (poly:^ h delta)))))
      (set! g (car (last-pair u)))
      (set! h (cond ((one? delta) g)
            ((zero? delta) h)
            (else (poly:/ (poly:^ g delta)
                  (poly:^ h (+ -1 delta)))))))))

(define (univ:gcd u v)
  (let* ((cu (univ:cont u))
     (cv (univ:cont v))
     (c (poly:gcd cu cv))
     (ppu (poly:/ u cu))
     (ppv (poly:/ v cv))
     (ans (if (shorter? ppv ppu)
          (univ:prs ppu ppv)
        (univ:prs ppv ppu))))
    (if (zero? (univ:degree ans (car u)))
    c
      (poly:* c (univ:primpart ans)))))

(define (poly:gcd p1 p2)
  (cond ((equal? p1 p2) p1)
    ((and (number? p1) (number? p2)) (gcd p1 p2))
    ((number? p1) (if (zero? p1)
              p2
            (apply poly:gcd* p1 (cdr p2))))
    ((number? p2) (if (zero? p2)
              p1
            (apply poly:gcd* p2 (cdr p1))))
    ((eq? (car p1) (car p2)) (univ:gcd p1 p2))
    ((var:> (car p2) (car p1)) (apply poly:gcd* p1 (cdr p2)))
    (else (apply poly:gcd* p2 (cdr p1)))))

(define (poly:gcd* . li)
  (let ((nums (remove-if-not number? li)))
    (if (null? nums)
    (reduce poly:gcd li)
    (let ((gnum (reduce gcd nums)))
      (if (= 1 gnum) 1
          (reduce-init poly:gcd gnum (remove-if number? li)))))))

(define (univ:cont p) (apply poly:gcd* (cdr p)))
(define (univ:primpart p) (poly:/ p (univ:cont p)))
(define (poly:num-cont p)
  (if (number? p) p
      (do ((l (cdr p) (cdr l))
       (n (poly:num-cont (cadr p))
          (gcd n (poly:num-cont (cadr l)))))
      ((or (= 1 n) (null? (cdr l))) n))))

(define (list-ref? l n)
  (cond ((null? l) #f)
    ((zero? n) (car l))
    (else (list-ref? (cdr l) (- n 1)))))

(define (univ:coeff p ord) (or (list-ref? (cdr p) ord) 0))
(define (poly:coeff p var ord)
  (cond ((or (number? p) (var:> var (car p)))
     (if (zero? ord) p 0))
    ((eq? var (car p)) (univ:coeff p ord))
    (else
     (univ:norm0 (car p)
             (map-no-end-0s (lambda (c) (poly:coeff c var ord))
                    (cdr p))))))

(define (poly:subst0 old e) (poly:coeff e old 0))

(define const:promote list)

(define (promote var p)
  (if (eq? var (car p))
      p
      (let ((dgr (poly:degree p var)))
    (do ((i dgr (+ -1 i))
         (ol (list (poly:coeff p var dgr))
         (cons (poly:coeff p var (+ -1 i)) ol)))
        ((zero? i) (cons var ol))))))

;;;this is bummed if v has higher priority than any variable in (cdr p)
(define (univ:demote p)
  (if (number? p)
      p
    (let ((v (car p)))
      (if (every (lambda (cof) (or (number? cof) (var:> v (car cof))))
         (cdr p))
      p
    (poly:+ (cadr p)
        (do ((trms (cddr p) (cdr trms))
             (sum 0)
             (mon (list v 0 1) (cons v (cons 0 (cdr mon)))))
            ((null? trms) sum)
            (set! sum (poly:+ sum (poly:* mon (car trms))))))))))

(define (sylvester p1 p2 var)
  (set! p1 (promote var p1))
  (set! p2 (promote var p2))
  (let ((d1 (univ:degree p1 var))
    (d2 (univ:degree p2 var))
    (m (list)))
    (do ((i d1 (+ -1 i))
     (row (nconc (make-list (+ -1 d1) 0) (reverse (cdr p2)))
          (append (cdr row) (list 0))))
    ((<= i 1) (set! m (cons row m)))
    (set! m (cons row m)))
    (do ((i d2 (+ -1 i))
     (row (nconc (make-list (+ -1 d2) 0) (reverse (cdr p1)))
          (append (cdr row) (list 0))))
    ((<= i 1) (set! m (cons row m)))
    (set! m (cons row m)))
    m))

;;; Bareiss's integer preserving gaussian elimination.
;;; Bareiss, E.H.: Sylvester's identity and multistep
;;; integer-preserving Gaussian elimination. Mathematics of
;;; Computation 22, 565-578, 1968.
;;; as related by:
;;; Akritas, A.G.: Exact Algorithms for the Matrix-Triangulation
;;; Subresultant PRS Method.  Computers and Mathematics, 145-155.
;;; Springer Verlag, 1989.
(define (bareiss m)
  4)

(define (poly:resultant p1 p2 var)
  (let ((u1 (promote var p1))
    (u2 (promote var p2)))
    (or (not (zero? (univ:degree u1 var)))
    (not (zero? (univ:degree u2 var)))
    (math:error var 'does-not-appear-in- p1 'or- p2))
    (let ((res (cond ((zero? (univ:degree u1 var)) p1)
             ((zero? (univ:degree u2 var)) p2)
             ((shorter? u1 u2) (univ:prs u2 u1))
             (else (univ:prs u1 u2)))))
      (if (zero? (univ:degree res var)) res
      0))))

(define (poly:elim2 p1 p2 var)
  (let* ((u1 (promote var p1))
     (u2 (promote var p2))
     (pg (poly:gcd (car (last-pair u1)) (car (last-pair u2)))))
    (or (not (zero? (univ:degree u1 var)))
    (not (zero? (univ:degree u2 var)))
    (math:error var 'does-not-appear-in- p1 'or- p2))
    (let* ((res (cond ((zero? (univ:degree u1 var)) p1)
              ((zero? (univ:degree u2 var)) p2)
              ((shorter? u1 u2) (univ:prs u2 u1))
              (else (univ:prs u1 u2))))
       (e (if (zero? (univ:degree res var)) res
          0)))
      (if (or (number? pg)) e
      (let ((q (poly:/ e pg)))
        (if (number? q) e (univ:primpart q)))))))

(define (poly-mod-number poly modulus)
  (if (number? poly)
      (modulo poly modulus)
    (cons (car poly)
      (map-no-end-0s (lambda (x) (poly-mod-number x modulus))
             (cdr poly)))))

(define (poly:prem dividend divisor)
  (if (number? divisor)
      (poly-mod-number dividend divisor)
      (let ((var (car divisor)))
    (if (> (poly:degree divisor var) (poly:degree dividend var))
        dividend
        (univ:demote (univ:prem (promote var dividend)
                    (promote var divisor)))))))

;;;; VERIFICATION TESTS
(define (poly:test)
  (define a (sexp->var 'a))
  (define b (sexp->var 'b))
  (define c (sexp->var 'c))
  (define x (sexp->var 'x))
  (define y (sexp->var 'y))
  (test (list a 0 -2)
    poly:gcd
    (list a 0 -2)
    (list a 0 0 -2))
  (test (list x (list a 0 1) 1)
    poly:gcd
    (list x (list a 0 0 -1) 0 1)
    (list x (list a 0 0 1) (list a 0 2) 1))
  (test (list x 0 (list a 0 1))
    poly:gcd
    (list x 0 (list a 0 0 1))
    (list x 0 0 (list a 0 1)))
  (test (list x (list b 0 0 1) 0 (list b 1 2) (list a 0 1) 1)
    poly:resultant
    (list y (list x (list b 0 1) 0 1) (list x 0 1))
    (list y (list x 1 (list a 0 1)) 0 1)
    y)
  (test (list y (list b 0 0 1) 0 (list b 1 2) (list a 0 1) 1)
    poly:resultant
    (list y (list b 0 1) (list x 0 1) 1)
    (list y (list x 1 0 1) (list a 0 1))
    x)
  (test 1
     poly:gcd
     (list x -5 2 8 -3 -3 1 1)
     (list x 21 -9 -4 5 3))
  (test 1
     poly:gcd
     (list x -5 2 8 -3 -3 0 1 0 1)
     (list x 21 -9 -4 0 5 0 3))
  'done)