- Timestamp:
- 2015-06-22T09:57:51-07:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/f4grobner/polynomial.lisp
r3122 r3123 367 367 368 368 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 459 370 460 371 (defun spoly (ring-and-order f g … … 498 409 (reduce (ring-gcd ring) (mapcar #'term-coeff (rest (poly-termlist p))) :initial-value (poly-lc p))) 499 410 500 (defun read-infix-form (&key (stream t))501 "Parser of infix expressions with integer/rational coefficients502 The parser will recognize two kinds of polynomial expressions:503 504 - polynomials in fully expanded forms with coefficients505 written in front of symbolic expressions; constants can be optionally506 enclosed in (); for example, the infix form507 X^2-Y^2+(-4/3)*U^2*W^3-5508 parses to509 (+ (- (EXPT X 2) (EXPT Y 2)) (* (- (/ 4 3)) (EXPT U 2) (EXPT W 3)) (- 5))510 511 - lists of polynomials; for example512 [X-Y, X^2+3*Z]513 parses to514 (:[ (- X Y) (+ (EXPT X 2) (* 3 Z)))515 where the first symbol [ marks a list of polynomials.516 517 -other infix expressions, for example518 [(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-string523 (concatenate 'string524 "#I("525 (with-output-to-string (s)526 (loop527 (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 &key534 (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 a539 list of polynomials in variables VARS, over the ring RING. The540 polynomial or list of polynomials is returned, with terms in each541 polynomial ordered according to monomial order ORDER."542 (poly-eval (read-infix-form :stream stream) vars ring order))543 544 (defun string->poly (str vars545 &optional546 (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 the554 returned value is ((MONOM[0] . COEFF[0]) (MONOM[1] . COEFF[1]) ...), where555 MONOM[I] is a list of exponents in the monomial and COEFF[I] is the556 corresponding coefficient in the ring."557 (cond558 ((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 vars564 &optional565 (ring +ring-of-integers+)566 (order #'lex>))567 "Convert a string STR representing a polynomial or polynomial list to568 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))583 411 |#
Note:
See TracChangeset
for help on using the changeset viewer.