source: branches/f4grobner/.junk/poly-eval.lisp@ 1058

(defun convert-number (number-or-poly n)
  "Returns NUMBER-OR-POLY, if it is a polynomial. If it is a number,
it converts it to the constant monomial in N variables. If the result
is a number then convert it to a polynomial in N variables."
  (if (numberp number-or-poly)
      (make-poly-from-termlist (list (make-term (make-monom :dimension n) number-or-poly)))
      number-or-poly))

(defun $poly+ (ring-and-order p q n)
  "Add two polynomials P and Q, where each polynomial is either a
numeric constant or a polynomial in internal representation. If the
result is a number then convert it to a polynomial in N variables."
  (poly-add ring-and-order (convert-number p n) (convert-number q n)))

(defun $poly- (ring-and-order p q n)
  "Subtract two polynomials P and Q, where each polynomial is either a
numeric constant or a polynomial in internal representation. If the
result is a number then convert it to a polynomial in N variables."
  (poly-sub ring-and-order (convert-number p n) (convert-number q n)))

(defun $minus-poly (ring p n)
  "Negation of P is a polynomial is either a numeric constant or a
polynomial in internal representation. If the result is a number then
convert it to a polynomial in N variables."
  (poly-uminus ring (convert-number p n)))

(defun $poly* (ring-and-order p q n)
  "Multiply two polynomials P and Q, where each polynomial is either a
numeric constant or a polynomial in internal representation. If the
result is a number then convert it to a polynomial in N variables."
  (poly-mul ring-and-order (convert-number p n) (convert-number q n)))

(defun $poly/ (ring p q)
  "Divide a polynomials P which is either a numeric constant or a
polynomial in internal representation, by a number Q."
  (if (numberp p)
      (common-lisp:/ p q)
      (scalar-times-poly ring (common-lisp:/ q) p)))
    
(defun $poly-expt (ring-and-order p l n)
  "Raise polynomial P, which is a polynomial in internal
representation or a numeric constant, to power L. If P is a number,
convert the result to a polynomial in N variables."
  (poly-expt ring-and-order (convert-number p n) l))


(defun variable-basis (ring n &aux (basis (make-list n)))
  "Generate a list of polynomials X[i], i=0,1,...,N-1."
  (dotimes (i n basis)
    (setf (elt basis i) (make-variable ring n i))))


(defun poly-eval-1 (expr vars &optional (ring *ring-of-integers*) (order #'lex>)
                    &aux
                      (ring-and-order (make-ring-and-order :ring ring :order order))
                      (n (length vars))
                      (basis (variable-basis ring (length vars))))
  "Evaluate an expression EXPR as polynomial by substituting operators
+ - * expt with corresponding polynomial operators and variables VARS
with the corresponding polynomials in internal form.  We use special
versions of binary operators $poly+, $poly-, $minus-poly, $poly* and
$poly-expt which work like the corresponding functions in the POLY
package, but accept scalars as arguments as well. The result is a
polynomial in internal form.  This operation is somewhat similar to
the function EXPAND in CAS."
  (cond
    ((numberp expr)
     (cond
       ((zerop expr) NIL)
       (t (make-poly-from-termlist (list (make-term (make-monom :dimension n) expr))))))
    ((symbolp expr)
     (nth (position expr vars) basis))
    ((consp expr)
     (case (car expr)
       (expt
        (if (= (length expr) 3)
            ($poly-expt ring-and-order
                        (poly-eval-1 (cadr expr) vars ring order)
                        (caddr expr)
                        n)
            (error "Too many arguments to EXPT")))
       (/
        (if (and (= (length expr) 3)
                 (numberp (caddr expr)))
            ($poly/ ring (cadr expr) (caddr expr))
            (error "The second argument to / must be a number")))
       (otherwise
        (let ((r (mapcar
                  #'(lambda (e) (poly-eval-1 e vars ring order))
                  (cdr expr))))
          (ecase (car expr)
            (+ (reduce #'(lambda (p q) ($poly+ ring-and-order p q n)) r))
            (-
             (if (endp (cdr r))
                 ($minus-poly ring (car r) n)
                 ($poly- ring-and-order
                         (car r)
                         (reduce #'(lambda (p q) ($poly+ ring-and-order p q n)) (cdr r))
                         n)))
            (*
             (reduce #'(lambda (p q) ($poly* ring-and-order p q n)) r))
            )))))))


             
(defun poly-eval (expr vars &optional (order #'lex>) (ring *ring-of-integers*))
  "Evaluate an expression EXPR, which should be a polynomial
expression or a list of polynomial expressions (a list of expressions
marked by prepending keyword :[ to it) given in Lisp prefix notation,
in variables VARS, which should be a list of symbols. The result of
the evaluation is a polynomial or a list of polynomials (marked by
prepending symbol '[) in the internal alist form. This evaluator is
used by the PARSE package to convert input from strings directly to
internal form."
  (cond
   ((numberp expr)
    (unless (zerop expr)
      (make-poly-from-termlist 
       (list (make-term (make-monom :dimension (length vars)) expr)))))
   ((or (symbolp expr) (not (eq (car expr) :[)))
    (poly-eval-1 expr vars ring order))
   (t (cons '[ (mapcar #'(lambda (p) (poly-eval-1 p vars ring order)) (rest expr))))))