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 69 for branches


Ignore:
Timestamp:
2015-06-05T11:30:24-07:00 (9 years ago)
Author:
Marek Rychlik
Message:

* empty log message *

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/f4grobner/grobner.lisp

    r66 r69  
    373373
    374374
    375 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    376 ;;
    377 ;; Operations in ideal theory
    378 ;;
    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 basis
    393 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 plist
    397           (remove-if #'(lambda (p)
    398                          (poly-depends-p p i))
    399                      plist))))
    400 
    401 (defun elimination-ideal (ring flist k
    402                           &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 is
    411 defined as the set of polynomials H such that for all polynomials W in
    412 J the polynomial W*H belongs to I."
    413   (cond
    414     ((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-poly
    419                 (list (make-term
    420                        (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     (t
    426      (ideal-intersection ring
    427                          (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}), where
    433 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-contract
    441           (ring-intersection
    442            (reduced-grobner
    443             ring
    444             (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             0
    451             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 PRED
    457 and their coefficients must be compatible with the RING structure
    458 defined in the COEFFICIENT-RING package."
    459   (cond
    460     ((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     (t
    466      (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-poly
    471           ring
    472           (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 generates
    483 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 the
    489 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^inf
    493 (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 ideal
    497 generated by a polynomial list F in the ideal generated by a single
    498 polynomial P. The saturation ideal is defined as the set of
    499 polynomials H such for some natural number n (* (EXPT P N) H) is in the ideal
    500 F. Geometrically, over an algebraically closed field, this is the set
    501 of polynomials in the ideal generated by F which do not identically
    502 vanish on the variety of P."
    503   (mapcar
    504    #'poly-contract
    505    (ring-intersection
    506     (reduced-grobner
    507      ring
    508      (saturation-extension-1 ring f p)
    509      start top-reduction-only)
    510     1)))
    511 
    512 
    513 
    514 ;; Calculate F : p1^inf : p2^inf : ... : ps^inf
    515 (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 a
    517 sequence of successive saturations in the polynomials
    518 of the polynomial list PLIST of the ideal generated by the
    519 polynomial list F."
    520   (cond
    521    ((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                          &aux
    527                          (k (length g))
    528                          (*monomial-order* (or *elimination-order*
    529                                                (elimination-order k))))
    530   "Returns the reduced Grobner basis of the saturation of the ideal
    531 generated by a polynomial list F in the ideal generated a polynomial
    532 list G. The saturation ideal is defined as the set of polynomials H
    533 such for some natural number n and some P in the ideal generated by G
    534 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 in
    536 the ideal generated by F which do not identically vanish on the
    537 variety of G."
    538   (mapcar
    539    #'(lambda (q) (poly-contract q k))
    540    (ring-intersection
    541     (reduced-grobner ring
    542                      (polysaturation-extension ring f g)
    543                      start
    544                      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 a
    549 successive applications of IDEAL-SATURATION to F and lists of
    550 polynomials in the list IDEAL-LIST."
    551   (cond
    552    ((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 
    557375
    558376;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     
    640458  "The ring of coefficients, over which all polynomials
    641459are assumed to be defined.")
    642 
    643 
    644 
    645 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    646 ;;
    647 ;; Maxima expression parsing
    648 ;;
    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   (cond
    657    ((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     (cond
    667      ((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 converted
    673       ;;to coefficients intact
    674       (coerce-coeff *maxima-ring* expr vars))
    675      (t
    676       (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         (mtimes
    680          (if (endp (cddr expr))         ;unary
    681              (parse (cdr expr))
    682            (reduce #'(lambda (p q) (poly-mul *maxima-ring* p q)) (parse-list (cdr expr)))))
    683         (mexpt
    684          (cond
    685           ((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 ring
    691            ;; Non-integer power means non-polynomial coefficient
    692            (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         (otherwise
    699          (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 
    709460
    710461
Note: See TracChangeset for help on using the changeset viewer.