source: branches/f4grobner/polynomial-eval.lisp

;;; -*-  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.  
;;;                                                                              
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package "POLYNOMIAL")

(defvar *coefficient-class* 'integer-ring
  "The default class in which coefficients are created from
NUMBER tokens.")

(defun poly-eval (expr vars order &optional (coefficient-class *coefficient-class*))
  "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."
  (labels ((p-eval (p) (poly-eval p vars order coefficient-class))
           (p-eval-list (plist) (mapcar #'p-eval plist)))
    (cond 
      ((eq expr 0) 
       (make-instance 'poly))
      ((member expr vars :test #'equalp)
       (let ((pos (position expr vars :test #'equalp)))
         (make-poly-variable (length vars) pos)))
      ((numberp expr)
       (make-poly-constant (length vars) (make-instance coefficient-class :value expr)))
      ((eq (car expr) +list-marker+)
       (cons +list-marker+ (p-eval-list (cdr expr))))
      (t
       (case (car expr)
         (+ (reduce #'add (p-eval-list (cdr expr))))
         (- (apply #'subtract (p-eval-list (cdr expr))))
         (*
          (if (endp (cddr expr))        ;unary
              (p-eval (cadr expr))
              (apply #'multiply (p-eval-list (cdr expr)))))
         (/ 
          ;; A polynomial can be divided by a scalar
          (cond 
            ((endp (cddr expr))
             ;; A special case (/ ?), the inverse
             (divide (cadr expr)))
            (t
             (let ((num (p-eval (cadr expr)))
                   (denom-inverse (apply #'divide (mapcar #'p-eval (cddr expr)))))
               (multiply 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 (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
             expr)
            (t (universal-expt (p-eval (cadr expr)) (caddr expr)))))
         (otherwise
          (error "Cannot evaluate as polynomial: ~A" expr)))))))