source: branches/f4grobner/pol.lisp@ 2001

;;; -*-  Mode: Lisp -*- 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                              
;;;  Copyright (C) 1999, 2002, 2009, 2015 Marek Rychlik <rychlik@u.arizona.edu>          
;;;                                                                              
;;;  This program is free software; you can redistribute it and/or modify        
;;;  it under the terms of the GNU General Public License as published by        
;;;  the Free Software Foundation; either version 2 of the License, or           
;;;  (at your option) any later version.                                         
;;;                                                                              
;;;  This program is distributed in the hope that it will be useful,             
;;;  but WITHOUT ANY WARRANTY; without even the implied warranty of              
;;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               
;;;  GNU General Public License for more details.                                
;;;                                                                              
;;;  You should have received a copy of the GNU General Public License           
;;;  along with this program; if not, write to the Free Software                 
;;;  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  
;;;                                                                              
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Polynomials implemented in CLOS
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; A polynomial is an collection of terms. A
;; term has a monomial and a coefficient.
;;
;; A polynomial can be represented by an s-expp
;; (EXPR . VARS) where EXPR is an arithmetical formula
;; recursively built of the arithmetical operations,
;; and VARS are the variables of the polynomial.
;; If a subtree of this s-exp is not an arithmetical
;; operator +, -, *, expt, and is not a member
;; of VARS then it represents a scalar expression
;; which the Lisp reader must know how to convert
;; into an object for which can be multiplied by a variable,
;; subject to commutativity and associativity rules.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defpackage "POL"
  (:use :cl :ring :ring-and-order :monom :order :term :termlist :infix)
  (:export "POLY" 
           "POLY-TERMLIST" 
           "POLY-SUGAR" 
           "POLY-RESET-SUGAR"
           "POLY-LT"
           "MAKE-POLY-FROM-TERMLIST"
           "MAKE-POLY-ZERO"
           "MAKE-POLY-VARIABLE"
           "POLY-UNIT"
           "POLY-LM"
           "POLY-SECOND-LM"
           "POLY-SECOND-LT"
           "POLY-LC"
           "POLY-SECOND-LC"
           "POLY-ZEROP"
           "POLY-LENGTH"
           "SCALAR-TIMES-POLY"
           "SCALAR-TIMES-POLY-1"
           "MONOM-TIMES-POLY"
           "TERM-TIMES-POLY"
           "POLY-ADD"
           "POLY-SUB"
           "POLY-UMINUS"
           "POLY-MUL"
           "POLY-EXPT"
           "POLY-APPEND"
           "POLY-NREVERSE"
           "POLY-REVERSE"
           "POLY-CONTRACT"
           "POLY-EXTEND"
           "POLY-ADD-VARIABLES"
           "POLY-LIST-ADD-VARIABLES"
           "POLY-STANDARD-EXTENSION"
           "SATURATION-EXTENSION"
           "POLYSATURATION-EXTENSION"
           "SATURATION-EXTENSION-1"
           "COERCE-COEFF"
           "POLY-EVAL"
           "POLY-EVAL-SCALAR"
           "SPOLY"
           "POLY-PRIMITIVE-PART"
           "POLY-CONTENT"
           "READ-INFIX-FORM"
           "READ-POLY"
           "STRING->POLY"
           "POLY->ALIST"
           "STRING->ALIST"
           "POLY-EQUAL-NO-SUGAR-P"
           "POLY-SET-EQUAL-NO-SUGAR-P"
           "POLY-LIST-EQUAL-NO-SUGAR-P"
           ))

(in-package :pol)

(proclaim '(optimize (speed 0) (space 0) (safety 3) (debug 3)))

(defclass poly ()
  ((termlist :initarg :termlist :accessor termlist))
  (:default-initargs :termlist nil))

(defun make-poly-from-termlist (termlist)
  (make-instance 'poly :termlist termlist :sugar sugar))

(defun make-poly-zero (&aux (termlist nil) (sugar -1))
  (make-instance 'poly :termlist termlist))

(defun make-poly-variable (ring nvars pos &optional (power 1)
                           &aux
                             (termlist (list
                                        (make-term-variable ring nvars pos power))))
  (make-instance 'poly :termlist termlist))

(defun poly-unit (ring dimension
                  &aux
                    (termlist (termlist-unit ring dimension)))
  (make-instance 'poly :termlist termlist))


(defmethod print-object ((self poly) stream)
  (princ (slot-value self 'termlist)))

(defmethod poly-termlist ((self poly))
  (slot-value self 'termlist))

(defmethod (setf poly-termlist) (new-value (poly poly))
  (setf (slot-value poly 'termlist) new-value))

(defmethod poly-add ((p poly) (q poly)))

(defmethod poly-sub ((p poly) (q poly)))

(defmethod poly-uminus ((self poly)))

(defmethod poly-mul ((p poly) (poly q)))

(defmethod poly-expt ((self poly) n))

(defmethod initialize-instance :after ((self poly) expr vars)
  "Evaluate Lisp form EXPR to a polynomial or a list of polynomials in
variables VARS. Return the resulting polynomial or list of
polynomials.  Standard arithmetical operators in form EXPR are
replaced with their analogues in the ring of polynomials, and the
resulting expression is evaluated, resulting in a polynomial or a list
of polynomials in internal form. A similar operation in another computer
algebra system could be called 'expand' or so."
    (cond
      ((null expr) (error "Empty expression"))
      ((eql expr 0) (make-poly-zero))
      ((member expr vars :test #'equalp)
       (let ((pos (position expr vars :test #'equalp)))
         (make-poly-variable ring (length vars) pos)))
      ((atom expr)
       (scalar->poly ring expr vars))
      ((eq (car expr) list-marker)
       (cons list-marker (p-eval-list (cdr expr))))
      (t
       (case (car expr)
         (+ (reduce #'p-add (p-eval-list (cdr expr))))
         (- (case (length expr)
              (1 (make-poly-zero))
              (2 (poly-uminus ring (p-eval (cadr expr))))
              (3 (poly-sub ring-and-order (p-eval (cadr expr)) (p-eval (caddr expr))))
              (otherwise (poly-sub ring-and-order (p-eval (cadr expr))
                                   (reduce #'p-add (p-eval-list (cddr expr)))))))
         (*
          (if (endp (cddr expr))                ;unary
              (p-eval (cdr expr))
              (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (p-eval-list (cdr expr)))))
         (/ 
          ;; A polynomial can be divided by a scalar
          (cond 
            ((endp (cddr expr))
             ;; A special case (/ ?), the inverse
             (scalar->poly ring (apply (ring-div ring) (cdr expr)) vars))
            (t
             (let ((num (p-eval (cadr expr)))
                   (denom-inverse (apply (ring-div ring)
                                         (cons (funcall (ring-unit ring)) 
                                               (mapcar #'p-eval-scalar (cddr expr))))))
               (scalar-times-poly ring denom-inverse num)))))
         (expt
          (cond
            ((member (cadr expr) vars :test #'equalp)
             ;;Special handling of (expt var pow)
             (let ((pos (position (cadr expr) vars :test #'equalp)))
               (make-poly-variable ring (length vars) pos (caddr expr))))
            ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
             ;; Negative power means division in coefficient ring
             ;; Non-integer power means non-polynomial coefficient
             (scalar->poly ring expr vars))
            (t (poly-expt ring-and-order (p-eval (cadr expr)) (caddr expr)))))
         (otherwise
          (scalar->poly ring expr vars)))))))

(defun poly-eval-scalar (expr 
                         &optional
                           (ring +ring-of-integers+)
                         &aux 
                           (order #'lex>))
  "Evaluate a scalar expression EXPR in ring RING."
  (declare (type ring ring))
  (poly-lc (poly-eval expr nil ring order)))