close Warning: Can't synchronize with repository "(default)" (The repository directory has changed, you should resynchronize the repository with: trac-admin $ENV repository resync '(default)'). Look in the Trac log for more information.

Changeset 3123 for branches


Ignore:
Timestamp:
2015-06-22T09:57:51-07:00 (10 years ago)
Author:
Marek Rychlik
Message:

* empty log message *

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/f4grobner/polynomial.lisp

    r3122 r3123  
    367367
    368368
    369 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    370 ;;
    371 ;; Evaluation of polynomial (prefix) expressions
    372 ;;
    373 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    374 
    375 (defun coerce-coeff (ring expr vars)
    376   "Coerce an element of the coefficient ring to a constant polynomial."
    377   ;; Modular arithmetic handler by rat
    378   (declare (type ring ring))
    379   (make-poly-from-termlist (list (make-term :monom (make-monom :dimension (length vars))
    380                                             :coeff (funcall (ring-parse ring) expr)))
    381                            0))
    382 
    383 (defun poly-eval (expr vars
    384                   &optional
    385                     (ring +ring-of-integers+)
    386                     (order #'lex>)
    387                     (list-marker :[)
    388                   &aux
    389                     (ring-and-order (make-ring-and-order :ring ring :order order)))
    390   "Evaluate Lisp form EXPR to a polynomial or a list of polynomials in
    391 variables VARS. Return the resulting polynomial or list of
    392 polynomials.  Standard arithmetical operators in form EXPR are
    393 replaced with their analogues in the ring of polynomials, and the
    394 resulting expression is evaluated, resulting in a polynomial or a list
    395 of polynomials in internal form. A similar operation in another computer
    396 algebra system could be called 'expand' or so."
    397   (declare (type ring ring))
    398   (labels ((p-eval (arg) (poly-eval arg vars ring order))
    399            (p-eval-scalar (arg) (poly-eval-scalar arg))
    400            (p-eval-list (args) (mapcar #'p-eval args))
    401            (p-add (x y) (poly-add ring-and-order x y)))
    402     (cond
    403       ((null expr) (error "Empty expression"))
    404       ((eql expr 0) (make-poly-zero))
    405       ((member expr vars :test #'equalp)
    406        (let ((pos (position expr vars :test #'equalp)))
    407          (make-poly-variable ring (length vars) pos)))
    408       ((atom expr)
    409        (coerce-coeff ring expr vars))
    410       ((eq (car expr) list-marker)
    411        (cons list-marker (p-eval-list (cdr expr))))
    412       (t
    413        (case (car expr)
    414          (+ (reduce #'p-add (p-eval-list (cdr expr))))
    415          (- (case (length expr)
    416               (1 (make-poly-zero))
    417               (2 (poly-uminus ring (p-eval (cadr expr))))
    418               (3 (poly-sub ring-and-order (p-eval (cadr expr)) (p-eval (caddr expr))))
    419               (otherwise (poly-sub ring-and-order (p-eval (cadr expr))
    420                                    (reduce #'p-add (p-eval-list (cddr expr)))))))
    421          (*
    422           (if (endp (cddr expr))                ;unary
    423               (p-eval (cdr expr))
    424               (reduce #'(lambda (p q) (poly-mul ring-and-order p q)) (p-eval-list (cdr expr)))))
    425          (/
    426           ;; A polynomial can be divided by a scalar
    427           (cond
    428             ((endp (cddr expr))
    429              ;; A special case (/ ?), the inverse
    430              (coerce-coeff ring (apply (ring-div ring) (cdr expr)) vars))
    431             (t
    432              (let ((num (p-eval (cadr expr)))
    433                    (denom-inverse (apply (ring-div ring)
    434                                          (cons (funcall (ring-unit ring))
    435                                                (mapcar #'p-eval-scalar (cddr expr))))))
    436                (scalar-times-poly ring denom-inverse num)))))
    437          (expt
    438           (cond
    439             ((member (cadr expr) vars :test #'equalp)
    440              ;;Special handling of (expt var pow)
    441              (let ((pos (position (cadr expr) vars :test #'equalp)))
    442                (make-poly-variable ring (length vars) pos (caddr expr))))
    443             ((not (and (integerp (caddr expr)) (plusp (caddr expr))))
    444              ;; Negative power means division in coefficient ring
    445              ;; Non-integer power means non-polynomial coefficient
    446              (coerce-coeff ring expr vars))
    447             (t (poly-expt ring-and-order (p-eval (cadr expr)) (caddr expr)))))
    448          (otherwise
    449           (coerce-coeff ring expr vars)))))))
    450 
    451 (defun poly-eval-scalar (expr
    452                          &optional
    453                            (ring +ring-of-integers+)
    454                          &aux
    455                            (order #'lex>))
    456   "Evaluate a scalar expression EXPR in ring RING."
    457   (declare (type ring ring))
    458   (poly-lc (poly-eval expr nil ring order)))
     369
    459370
    460371(defun spoly (ring-and-order f g
     
    498409  (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p)))
    499410
    500 (defun read-infix-form (&key (stream t))
    501   "Parser of infix expressions with integer/rational coefficients
    502 The parser will recognize two kinds of polynomial expressions:
    503 
    504 - polynomials in fully expanded forms with coefficients
    505   written in front of symbolic expressions; constants can be optionally
    506   enclosed in (); for example, the infix form
    507         X^2-Y^2+(-4/3)*U^2*W^3-5
    508   parses to
    509         (+ (- (EXPT X 2) (EXPT Y 2)) (* (- (/ 4 3)) (EXPT U 2) (EXPT W 3)) (- 5))
    510 
    511 - lists of polynomials; for example
    512         [X-Y, X^2+3*Z]         
    513   parses to
    514           (:[ (- X Y) (+ (EXPT X 2) (* 3 Z)))
    515   where the first symbol [ marks a list of polynomials.
    516 
    517 -other infix expressions, for example
    518         [(X-Y)*(X+Y)/Z,(X+1)^2]
    519 parses to:
    520         (:[ (/ (* (- X Y) (+ X Y)) Z) (EXPT (+ X 1) 2))
    521 Currently this function is implemented using M. Kantrowitz's INFIX package."
    522   (read-from-string
    523    (concatenate 'string
    524                 "#I("
    525                 (with-output-to-string (s)
    526                   (loop
    527                      (multiple-value-bind (line eof)
    528                          (read-line stream t)
    529                        (format s "~A" line)
    530                        (when eof (return)))))
    531                 ")")))
    532 
    533 (defun read-poly (vars &key
    534                          (stream t)
    535                          (ring +ring-of-integers+)
    536                          (order #'lex>))
    537   "Reads an expression in prefix form from a stream STREAM.
    538 The expression read from the strem should represent a polynomial or a
    539 list of polynomials in variables VARS, over the ring RING.  The
    540 polynomial or list of polynomials is returned, with terms in each
    541 polynomial ordered according to monomial order ORDER."
    542   (poly-eval (read-infix-form :stream stream) vars ring order))
    543 
    544 (defun string->poly (str vars
    545                      &optional
    546                        (ring +ring-of-integers+)
    547                        (order #'lex>))
    548   "Converts a string STR to a polynomial in variables VARS."
    549   (with-input-from-string (s str)
    550     (read-poly vars :stream s :ring ring :order order)))
    551 
    552 (defun poly->alist (p)
    553   "Convert a polynomial P to an association list. Thus, the format of the
    554 returned value is  ((MONOM[0] . COEFF[0]) (MONOM[1] . COEFF[1]) ...), where
    555 MONOM[I] is a list of exponents in the monomial and COEFF[I] is the
    556 corresponding coefficient in the ring."
    557   (cond
    558     ((poly-p p)
    559      (mapcar #'term->cons (poly-termlist p)))
    560     ((and (consp p) (eq (car p) :[))
    561      (cons :[ (mapcar #'poly->alist (cdr p))))))
    562 
    563 (defun string->alist (str vars
    564                       &optional
    565                         (ring +ring-of-integers+)
    566                         (order #'lex>))
    567   "Convert a string STR representing a polynomial or polynomial list to
    568 an association list (... (MONOM . COEFF) ...)."
    569   (poly->alist (string->poly str vars ring order)))
    570 
    571 (defun poly-equal-no-sugar-p (p q)
    572   "Compare polynomials for equality, ignoring sugar."
    573   (declare (type poly p q))
    574   (equalp (poly-termlist p) (poly-termlist q)))
    575 
    576 (defun poly-set-equal-no-sugar-p (p q)
    577   "Compare polynomial sets P and Q for equality, ignoring sugar."
    578   (null (set-exclusive-or  p q :test #'poly-equal-no-sugar-p )))
    579 
    580 (defun poly-list-equal-no-sugar-p (p q)
    581   "Compare polynomial lists P and Q for equality, ignoring sugar."
    582   (every #'poly-equal-no-sugar-p p q))
    583411|#
Note: See TracChangeset for help on using the changeset viewer.