The World of Computer Software/World_Of_Computer_Software-02-387-Vol-3of3.iso/j/jacal1a0.zip/jacal/poly.scm

Text File | 1992-12-23 | 16KB | 500 lines

;;; JACAL: Symbolic Mathematics System.        -*-scheme-*-
;;; Copyright 1989, 1990, 1991, 1992 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))
(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)))))

;;;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))))))
;;; replaces each var in poly with (proc var)
(define (poly_do-vars proc poly)
  (if (number? poly)
      poly
    (cons (proc (car poly))
      (map (lambda (b) (poly_do-vars proc b))
           (cdr poly)))))
(define (polys_do-vars proc polys)
  (bunch_map (lambda (poly) (poly_do-vars proc poly)) polys))
(define (poleqns_do-vars proc licits)
  (bunch_map (lambda (poly) (poly_do-vars proc (licit->poleqn poly))) licits))
;;; 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)))))
(define (polys_for-each-var proc polys)
  (bunch_for-each (lambda (poly) (poly_for-each-var proc poly)) polys))
(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 (> (leading-number p) 0) (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_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)