- Timestamp:
- 2015-06-05T11:30:24-07:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/f4grobner/grobner.lisp
r66 r69 373 373 374 374 375 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;376 ;;377 ;; Operations in ideal theory378 ;;379 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;380 381 ;; Does the term depend on variable K?382 (defun term-depends-p (term k)383 "Return T if the term TERM depends on variable number K."384 (monom-depends-p (term-monom term) k))385 386 ;; Does the polynomial P depend on variable K?387 (defun poly-depends-p (p k)388 "Return T if the term polynomial P depends on variable number K."389 (some #'(lambda (term) (term-depends-p term k)) (poly-termlist p)))390 391 (defun ring-intersection (plist k)392 "This function assumes that polynomial list PLIST is a Grobner basis393 and it calculates the intersection with the ring R[x[k+1],...,x[n]], i.e.394 it discards polynomials which depend on variables x[0], x[1], ..., x[k]."395 (dotimes (i k plist)396 (setf plist397 (remove-if #'(lambda (p)398 (poly-depends-p p i))399 plist))))400 401 (defun elimination-ideal (ring flist k402 &optional (top-reduction-only $poly_top_reduction_only) (start 0)403 &aux (*monomial-order*404 (or *elimination-order*405 (elimination-order k))))406 (ring-intersection (reduced-grobner ring flist start top-reduction-only) k))407 408 (defun colon-ideal (ring f g &optional (top-reduction-only $poly_top_reduction_only))409 "Returns the reduced Grobner basis of the colon ideal Id(F):Id(G),410 where F and G are two lists of polynomials. The colon ideal I:J is411 defined as the set of polynomials H such that for all polynomials W in412 J the polynomial W*H belongs to I."413 (cond414 ((endp g)415 ;;Id(G) consists of 0 only so W*0=0 belongs to Id(F)416 (if (every #'poly-zerop f)417 (error "First ideal must be non-zero.")418 (list (make-poly419 (list (make-term420 (make-monom (monom-dimension (poly-lm (find-if-not #'poly-zerop f)))421 :initial-element 0)422 (funcall (ring-unit ring))))))))423 ((endp (cdr g))424 (colon-ideal-1 ring f (car g) top-reduction-only))425 (t426 (ideal-intersection ring427 (colon-ideal-1 ring f (car g) top-reduction-only)428 (colon-ideal ring f (rest g) top-reduction-only)429 top-reduction-only))))430 431 (defun colon-ideal-1 (ring f g &optional (top-reduction-only $poly_top_reduction_only))432 "Returns the reduced Grobner basis of the colon ideal Id(F):Id({G}), where433 F is a list of polynomials and G is a polynomial."434 (mapcar #'(lambda (x) (poly-exact-divide ring x g)) (ideal-intersection ring f (list g) top-reduction-only)))435 436 437 (defun ideal-intersection (ring f g &optional (top-reduction-only $poly_top_reduction_only)438 &aux (*monomial-order* (or *elimination-order*439 #'elimination-order-1)))440 (mapcar #'poly-contract441 (ring-intersection442 (reduced-grobner443 ring444 (append (mapcar #'(lambda (p) (poly-extend p (make-monom 1 :initial-element 1))) f)445 (mapcar #'(lambda (p)446 (poly-append (poly-extend (poly-uminus ring p)447 (make-monom 1 :initial-element 1))448 (poly-extend p)))449 g))450 0451 top-reduction-only)452 1)))453 454 (defun poly-lcm (ring f g)455 "Return LCM (least common multiple) of two polynomials F and G.456 The polynomials must be ordered according to monomial order PRED457 and their coefficients must be compatible with the RING structure458 defined in the COEFFICIENT-RING package."459 (cond460 ((poly-zerop f) f)461 ((poly-zerop g) g)462 ((and (endp (cdr (poly-termlist f))) (endp (cdr (poly-termlist g))))463 (let ((m (monom-lcm (poly-lm f) (poly-lm g))))464 (make-poly-from-termlist (list (make-term m (funcall (ring-lcm ring) (poly-lc f) (poly-lc g)))))))465 (t466 (multiple-value-bind (f f-cont)467 (poly-primitive-part ring f)468 (multiple-value-bind (g g-cont)469 (poly-primitive-part ring g)470 (scalar-times-poly471 ring472 (funcall (ring-lcm ring) f-cont g-cont)473 (poly-primitive-part ring (car (ideal-intersection ring (list f) (list g) nil)))))))))474 475 ;; Do two Grobner bases yield the same ideal?476 (defun grobner-equal (ring g1 g2)477 "Returns T if two lists of polynomials G1 and G2, assumed to be Grobner bases,478 generate the same ideal, and NIL otherwise."479 (and (grobner-subsetp ring g1 g2) (grobner-subsetp ring g2 g1)))480 481 (defun grobner-subsetp (ring g1 g2)482 "Returns T if a list of polynomials G1 generates483 an ideal contained in the ideal generated by a polynomial list G2,484 both G1 and G2 assumed to be Grobner bases. Returns NIL otherwise."485 (every #'(lambda (p) (grobner-member ring p g2)) g1))486 487 (defun grobner-member (ring p g)488 "Returns T if a polynomial P belongs to the ideal generated by the489 polynomial list G, which is assumed to be a Grobner basis. Returns NIL otherwise."490 (poly-zerop (normal-form ring p g nil)))491 492 ;; Calculate F : p^inf493 (defun ideal-saturation-1 (ring f p start &optional (top-reduction-only $poly_top_reduction_only)494 &aux (*monomial-order* (or *elimination-order*495 #'elimination-order-1)))496 "Returns the reduced Grobner basis of the saturation of the ideal497 generated by a polynomial list F in the ideal generated by a single498 polynomial P. The saturation ideal is defined as the set of499 polynomials H such for some natural number n (* (EXPT P N) H) is in the ideal500 F. Geometrically, over an algebraically closed field, this is the set501 of polynomials in the ideal generated by F which do not identically502 vanish on the variety of P."503 (mapcar504 #'poly-contract505 (ring-intersection506 (reduced-grobner507 ring508 (saturation-extension-1 ring f p)509 start top-reduction-only)510 1)))511 512 513 514 ;; Calculate F : p1^inf : p2^inf : ... : ps^inf515 (defun ideal-polysaturation-1 (ring f plist start &optional (top-reduction-only $poly_top_reduction_only))516 "Returns the reduced Grobner basis of the ideal obtained by a517 sequence of successive saturations in the polynomials518 of the polynomial list PLIST of the ideal generated by the519 polynomial list F."520 (cond521 ((endp plist) (reduced-grobner ring f start top-reduction-only))522 (t (let ((g (ideal-saturation-1 ring f (car plist) start top-reduction-only)))523 (ideal-polysaturation-1 ring g (rest plist) (length g) top-reduction-only)))))524 525 (defun ideal-saturation (ring f g start &optional (top-reduction-only $poly_top_reduction_only)526 &aux527 (k (length g))528 (*monomial-order* (or *elimination-order*529 (elimination-order k))))530 "Returns the reduced Grobner basis of the saturation of the ideal531 generated by a polynomial list F in the ideal generated a polynomial532 list G. The saturation ideal is defined as the set of polynomials H533 such for some natural number n and some P in the ideal generated by G534 the polynomial P**N * H is in the ideal spanned by F. Geometrically,535 over an algebraically closed field, this is the set of polynomials in536 the ideal generated by F which do not identically vanish on the537 variety of G."538 (mapcar539 #'(lambda (q) (poly-contract q k))540 (ring-intersection541 (reduced-grobner ring542 (polysaturation-extension ring f g)543 start544 top-reduction-only)545 k)))546 547 (defun ideal-polysaturation (ring f ideal-list start &optional (top-reduction-only $poly_top_reduction_only))548 "Returns the reduced Grobner basis of the ideal obtained by a549 successive applications of IDEAL-SATURATION to F and lists of550 polynomials in the list IDEAL-LIST."551 (cond552 ((endp ideal-list) f)553 (t (let ((h (ideal-saturation ring f (car ideal-list) start top-reduction-only)))554 (ideal-polysaturation ring h (rest ideal-list) (length h) top-reduction-only)))))555 556 557 375 558 376 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; … … 640 458 "The ring of coefficients, over which all polynomials 641 459 are assumed to be defined.") 642 643 644 645 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;646 ;;647 ;; Maxima expression parsing648 ;;649 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;650 651 (defun equal-test-p (expr1 expr2)652 (alike1 expr1 expr2))653 654 (defun coerce-maxima-list (expr)655 "convert a maxima list to lisp list."656 (cond657 ((and (consp (car expr)) (eql (caar expr) 'mlist)) (cdr expr))658 (t expr)))659 660 (defun free-of-vars (expr vars) (apply #'$freeof `(,@vars ,expr)))661 662 (defun parse-poly (expr vars &aux (vars (coerce-maxima-list vars)))663 "Convert a maxima polynomial expression EXPR in variables VARS to internal form."664 (labels ((parse (arg) (parse-poly arg vars))665 (parse-list (args) (mapcar #'parse args)))666 (cond667 ((eql expr 0) (make-poly-zero))668 ((member expr vars :test #'equal-test-p)669 (let ((pos (position expr vars :test #'equal-test-p)))670 (make-variable *maxima-ring* (length vars) pos)))671 ((free-of-vars expr vars)672 ;;This means that variable-free CRE and Poisson forms will be converted673 ;;to coefficients intact674 (coerce-coeff *maxima-ring* expr vars))675 (t676 (case (caar expr)677 (mplus (reduce #'(lambda (x y) (poly-add *maxima-ring* x y)) (parse-list (cdr expr))))678 (mminus (poly-uminus *maxima-ring* (parse (cadr expr))))679 (mtimes680 (if (endp (cddr expr)) ;unary681 (parse (cdr expr))682 (reduce #'(lambda (p q) (poly-mul *maxima-ring* p q)) (parse-list (cdr expr)))))683 (mexpt684 (cond685 ((member (cadr expr) vars :test #'equal-test-p)686 ;;Special handling of (expt var pow)687 (let ((pos (position (cadr expr) vars :test #'equal-test-p)))688 (make-variable *maxima-ring* (length vars) pos (caddr expr))))689 ((not (and (integerp (caddr expr)) (plusp (caddr expr))))690 ;; Negative power means division in coefficient ring691 ;; Non-integer power means non-polynomial coefficient692 (mtell "~%Warning: Expression ~%~M~%contains power which is not a positive integer. Parsing as coefficient.~%"693 expr)694 (coerce-coeff *maxima-ring* expr vars))695 (t (poly-expt *maxima-ring* (parse (cadr expr)) (caddr expr)))))696 (mrat (parse ($ratdisrep expr)))697 (mpois (parse ($outofpois expr)))698 (otherwise699 (coerce-coeff *maxima-ring* expr vars)))))))700 701 (defun parse-poly-list (expr vars)702 (case (caar expr)703 (mlist (mapcar #'(lambda (p) (parse-poly p vars)) (cdr expr)))704 (t (merror "Expression ~M is not a list of polynomials in variables ~M."705 expr vars))))706 (defun parse-poly-list-list (poly-list-list vars)707 (mapcar #'(lambda (g) (parse-poly-list g vars)) (coerce-maxima-list poly-list-list)))708 709 460 710 461
Note:
See TracChangeset
for help on using the changeset viewer.