Changeset 3899

Timestamp:

2016-05-29T14:28:38-07:00 (8 years ago)

Author:

Marek Rychlik

Message:

* empty log message *

File:

• branches/f4grobner/polynomial.lisp (modified) (1 diff)

branches/f4grobner/polynomial.lisp

@@ -560 +560 @@
                    (poly-termlist self))))
   
+
+(defun poly-eval (expr vars order) 
+  "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))
+           (p-eval-scalar (p) (poly-eval p '() order))
+           (p-eval-list (plist) (mapcar #'p-eval plist)))
+    (cond 
+      ((eq expr 0) 
+       (make-instance 'poly :dimension (length vars)))
+      ((member expr vars :test #'equalp)
+       (let ((pos (position expr vars :test #'equalp)))
+         (make-poly-variable (length vars) pos)))
+      ((atom expr)
+       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))
+              (reduce #'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-scalar (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
+          expr))))))