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

;; Translate output from parse to a pure list form

;; assuming variables are VARS



(defun alist-form (plist vars)

  "Translates an expression PLIST, which should be a list of polynomials

in variables VARS, to an alist representation of a polynomial.

It returns the alist. See also PARSE-TO-ALIST."

  (cond

   ((endp plist) nil)

   ((eql (first plist) '[)

    (cons '[ (mapcar #'(lambda (x) (alist-form x vars)) (rest plist))))

   (t

    (assert (eql (car plist) '+))

    (alist-form-1 (rest plist) vars))))



(defun alist-form-1 (p vars

                         &aux (ht (make-hash-table

                                   :test #'equal :size 16))

                         stack)

  (dolist (term p)

    (assert (eql (car term) '*))

    (incf  (gethash (powers (cddr term) vars) ht 0) (second term)))

  (maphash #'(lambda (key value) (unless (zerop value)

                                   (push (cons key value) stack))) ht)

  stack)

    

(defun powers (monom vars

                     &aux (tab (pairlis vars (make-list (length vars)

                                                        :initial-element 0))))

  (dolist (e monom (mapcar #'(lambda (v) (cdr (assoc v tab))) vars))

    (assert (equal (first e) '^))

    (assert (integerp (third e)))

    (assert (= (length e) 3))

    (let ((x (assoc (second e) tab)))

      (if (null x) (error "Variable ~a not in the list of variables." 

                         (second e))

        (incf (cdr x) (third e))))))

  





(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 :monom (make-monom :dimension n) :coeff 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-poly-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 :monom (make-monom :dimension n) :coeff 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 :monom (make-monom :dimension (length vars)) :coeff 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))))))





(defun parse-string-to-alist (str vars)

  "Parse string STR and return a polynomial as a sorted association

list of pairs (MONOM . COEFFICIENT). For example:

(parse-string-to-alist \"[x^2-y^2+(-4/3)*u^2*w^3-5,y]\" '(x y u w))

 ([ (((0 0 2 3) . -4/3) ((0 2 0 0) . -1) ((2 0 0 0) . 1)

     ((0 0 0 0) . -5))

    (((0 1 0 0) . 1)))

The functions PARSE-TO-SORTED-ALIST and PARSE-STRING-TO-SORTED-ALIST

sort terms by the predicate defined in the ORDER package."

  (with-input-from-string (stream str)

    (parse-to-alist vars stream)))





(defun parse-to-sorted-alist (vars &optional (order #'lex>) (stream t))

  "Parses streasm STREAM and returns a polynomial represented as 

a sorted alist. For example:

(WITH-INPUT-FROM-STRING (S \"X^2-Y^2+(-4/3)*U^2*W^3-5\")

  (PARSE-TO-SORTED-ALIST '(X Y U W) S))

returns

(((2 0 0 0) . 1) ((0 2 0 0) . -1) ((0 0 2 3) . -4/3) ((0 0 0 0) . -5))

and

(WITH-INPUT-FROM-STRING (S \"X^2-Y^2+(-4/3)*U^2*W^3-5\")

  (PARSE-TO-SORTED-ALIST '(X Y U W) T #'GRLEX>) S)

returns

(((0 0 2 3) . -4/3) ((2 0 0 0) . 1) ((0 2 0 0) . -1) ((0 0 0 0) . -5))"

  (sort-poly (parse-to-alist vars stream) order))



(defun parse-string-to-sorted-alist (str vars &optional (order #'lex>))

  "Parse a string to a sorted alist form, the internal representation

of polynomials used by our system."

  (with-input-from-string (stream str)

    (parse-to-sorted-alist vars order stream)))



(defun sort-poly-1 (p order)

  "Sort the terms of a single polynomial P using an admissible monomial order ORDER.

Returns the sorted polynomial. Destructively modifies P."

  (sort p order :key #'first))



;; Sort a polynomial or polynomial list

(defun sort-poly (poly-or-poly-list &optional (order #'lex>))

  "Sort POLY-OR-POLY-LIST, which could be either a single polynomial

or a list of polynomials in internal alist representation, using

admissible monomial order ORDER. Each polynomial is sorted using

SORT-POLY-1."

  (cond

   ((eql poly-or-poly-list :syntax-error) nil)

   ((null poly-or-poly-list) nil)

   ((eql (car poly-or-poly-list) '[)

    (cons '[ (mapcar #'(lambda (p) (sort-poly-1 p order)) 

                     (rest poly-or-poly-list))))

   (t (sort-poly-1 poly-or-poly-list order))))