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 65


Ignore:
Timestamp:
2015-06-05T11:27:08-07:00 (10 years ago)
Author:
Marek Rychlik
Message:

* empty log message *

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/f4grobner/grobner.lisp

    r56 r65  
    293293  `(when $poly_grobner_debug (format *terminal-io* ,@args)))
    294294
    295 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    296 ;;
    297 ;; An implementation of Grobner basis
    298 ;;
    299 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    300 
    301 
    302 
    303 
    304 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    305 ;;
    306 ;; An implementation of the division algorithm
    307 ;;
    308 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    309 
    310 (defun grobner-op (ring c1 c2 m f g)
    311   "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial.
    312 Assume that the leading terms will cancel."
    313   #+grobner-check(funcall (ring-zerop ring)
    314                           (funcall (ring-sub ring)
    315                                    (funcall (ring-mul ring) c2 (poly-lc f))
    316                                    (funcall (ring-mul ring) c1 (poly-lc g))))
    317   #+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
    318   ;; Note that we can drop the leading terms of f ang g
    319   (poly-sub ring
    320             (scalar-times-poly-1 ring c2 f)
    321             (scalar-times-poly-1 ring c1 (monom-times-poly m g))))
    322 
    323 (defun poly-pseudo-divide (ring f fl)
    324   "Pseudo-divide a polynomial F by the list of polynomials FL. Return
    325 multiple values. The first value is a list of quotients A.  The second
    326 value is the remainder R. The third argument is a scalar coefficient
    327 C, such that C*F can be divided by FL within the ring of coefficients,
    328 which is not necessarily a field. Finally, the fourth value is an
    329 integer count of the number of reductions performed.  The resulting
    330 objects satisfy the equation: C*F= sum A[i]*FL[i] + R."
    331   (declare (type poly f) (list fl))
    332   (do ((r (make-poly-zero))
    333        (c (funcall (ring-unit ring)))
    334        (a (make-list (length fl) :initial-element (make-poly-zero)))
    335        (division-count 0)
    336        (p f))
    337       ((poly-zerop p)
    338        (debug-cgb "~&~3T~d reduction~:p" division-count)
    339        (when (poly-zerop r) (debug-cgb " ---> 0"))
    340        (values (mapcar #'poly-nreverse a) (poly-nreverse r) c division-count))
    341     (declare (fixnum division-count))
    342     (do ((fl fl (rest fl))                              ;scan list of divisors
    343          (b a (rest b)))
    344         ((cond
    345           ((endp fl)                                    ;no division occurred
    346            (push (poly-lt p) (poly-termlist r))         ;move lt(p) to remainder
    347            (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
    348            (pop (poly-termlist p))                      ;remove lt(p) from p
    349            t)
    350           ((monom-divides-p (poly-lm (car fl)) (poly-lm p)) ;division occurred
    351            (incf division-count)
    352            (multiple-value-bind (gcd c1 c2)
    353                (funcall (ring-ezgcd ring) (poly-lc (car fl)) (poly-lc p))
    354              (declare (ignore gcd))
    355              (let ((m (monom-div (poly-lm p) (poly-lm (car fl)))))
    356                ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
    357                (mapl #'(lambda (x)
    358                          (setf (car x) (scalar-times-poly ring c1 (car x))))
    359                      a)
    360                (setf r (scalar-times-poly ring c1 r)
    361                      c (funcall (ring-mul ring) c c1)
    362                      p (grobner-op ring c2 c1 m p (car fl)))
    363                (push (make-term m c2) (poly-termlist (car b))))
    364              t)))))))
    365 
    366 (defun poly-exact-divide (ring f g)
    367   "Divide a polynomial F by another polynomial G. Assume that exact division
    368 with no remainder is possible. Returns the quotient."
    369   (declare (type poly f g))
    370   (multiple-value-bind (quot rem coeff division-count)
    371       (poly-pseudo-divide ring f (list g))
    372     (declare (ignore division-count coeff)
    373              (list quot)
    374              (type poly rem)
    375              (type fixnum division-count))
    376     (unless (poly-zerop rem) (error "Exact division failed."))
    377     (car quot)))
    378 
    379 
    380 
    381 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    382 ;;
    383 ;; An implementation of the normal form
    384 ;;
    385 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    386 
    387 (defun normal-form-step (ring fl p r c division-count
    388                          &aux (g (find (poly-lm p) fl
    389                                        :test #'monom-divisible-by-p
    390                                        :key #'poly-lm)))
    391   (cond
    392    (g                                   ;division possible
    393     (incf division-count)
    394     (multiple-value-bind (gcd cg cp)
    395         (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
    396       (declare (ignore gcd))
    397       (let ((m (monom-div (poly-lm p) (poly-lm g))))
    398         ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
    399         (setf r (scalar-times-poly ring cg r)
    400               c (funcall (ring-mul ring) c cg)
    401               ;; p := cg*p-cp*m*g
    402               p (grobner-op ring cp cg m p g))))
    403     (debug-cgb "/"))
    404    (t                                                   ;no division possible
    405     (push (poly-lt p) (poly-termlist r))                ;move lt(p) to remainder
    406     (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
    407     (pop (poly-termlist p))                             ;remove lt(p) from p
    408     (debug-cgb "+")))
    409   (values p r c division-count))
    410 
    411 ;; Merge it sometime with poly-pseudo-divide
    412 (defun normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
    413   ;; Loop invariant: c*f0=sum ai*fi+r+f, where f0 is the initial value of f
    414   #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
    415   (do ((r (make-poly-zero))
    416        (c (funcall (ring-unit ring)))
    417        (division-count 0))
    418       ((or (poly-zerop f)
    419            ;;(endp fl)
    420            (and top-reduction-only (not (poly-zerop r))))
    421        (progn
    422          (debug-cgb "~&~3T~d reduction~:p" division-count)
    423          (when (poly-zerop r)
    424            (debug-cgb " ---> 0")))
    425        (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
    426        (values f c division-count))
    427     (declare (fixnum division-count)
    428              (type poly r))
    429     (multiple-value-setq (f r c division-count)
    430       (normal-form-step ring fl f r c division-count))))
     295
    431296
    432297
     
    440305;;
    441306;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    442 
    443 (defun buchberger-criterion (ring g)
    444   "Returns T if G is a Grobner basis, by using the Buchberger
    445 criterion: for every two polynomials h1 and h2 in G the S-polynomial
    446 S(h1,h2) reduces to 0 modulo G."
    447   (every
    448    #'poly-zerop
    449    (makelist (normal-form ring (spoly ring (elt g i) (elt g j)) g nil)
    450              (i 0 (- (length g) 2))
    451              (j (1+ i) (1- (length g))))))
    452307
    453308(defun grobner-test (ring g f)
     
    467322  t)
    468323
    469 
    470 
    471 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    472 ;;
    473 ;; Pair queue implementation
    474 ;;
    475 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    476 
    477 (defun sugar-pair-key (p q &aux (lcm (monom-lcm (poly-lm p) (poly-lm q)))
    478                                 (d (monom-sugar lcm)))
    479   "Returns list (S LCM-TOTAL-DEGREE) where S is the sugar of the S-polynomial of
    480 polynomials P and Q, and LCM-TOTAL-DEGREE is the degree of is LCM(LM(P),LM(Q))."
    481   (declare (type poly p q) (type monom lcm) (type fixnum d))
    482   (cons (max
    483          (+  (- d (monom-sugar (poly-lm p))) (poly-sugar p))
    484          (+  (- d (monom-sugar (poly-lm q))) (poly-sugar q)))
    485         lcm))
    486 
    487 (defstruct (pair
    488             (:constructor make-pair (first second
    489                                            &aux
    490                                            (sugar (car (sugar-pair-key first second)))
    491                                            (division-data nil))))
    492   (first nil :type poly)
    493   (second nil :type poly)
    494   (sugar 0 :type fixnum)
    495   (division-data nil :type list))
    496  
    497 ;;(defun pair-sugar (pair &aux (p (pair-first pair)) (q (pair-second pair)))
    498 ;;  (car (sugar-pair-key p q)))
    499 
    500 (defun sugar-order (x y)
    501   "Pair order based on sugar, ties broken by normal strategy."
    502   (declare (type cons x y))
    503   (or (< (car x) (car y))
    504       (and (= (car x) (car y))
    505            (< (monom-total-degree (cdr x))
    506               (monom-total-degree (cdr y))))))
    507 
    508 (defvar *pair-key-function* #'sugar-pair-key
    509   "Function that, given two polynomials as argument, computed the key
    510 in the pair queue.")
    511 
    512 (defvar *pair-order* #'sugar-order
    513   "Function that orders the keys of pairs.")
    514 
    515 (defun make-pair-queue ()
    516   "Constructs a priority queue for critical pairs."
    517   (make-priority-queue
    518    :element-type 'pair
    519    :element-key #'(lambda (pair) (funcall *pair-key-function* (pair-first pair) (pair-second pair)))
    520    :test *pair-order*))
    521 
    522 (defun pair-queue-initialize (pq f start
    523                               &aux
    524                               (s (1- (length f)))
    525                               (b (nconc (makelist (make-pair (elt f i) (elt f j))
    526                                                  (i 0 (1- start)) (j start s))
    527                                         (makelist (make-pair (elt f i) (elt f j))
    528                                                  (i start (1- s)) (j (1+ i) s)))))
    529   "Initializes the priority for critical pairs. F is the initial list of polynomials.
    530 START is the first position beyond the elements which form a partial
    531 grobner basis, i.e. satisfy the Buchberger criterion."
    532   (declare (type priority-queue pq) (type fixnum start))
    533   (dolist (pair b pq)
    534     (priority-queue-insert pq pair)))
    535 
    536 (defun pair-queue-insert (b pair)
    537   (priority-queue-insert b pair))
    538 
    539 (defun pair-queue-remove (b)
    540   (priority-queue-remove b))
    541 
    542 (defun pair-queue-size (b)
    543   (priority-queue-size b))
    544 
    545 (defun pair-queue-empty-p (b)
    546   (priority-queue-empty-p b))
    547 
    548 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    549 ;;
    550 ;; Buchberger Algorithm Implementation
    551 ;;
    552 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    553 
    554 (defun buchberger (ring f start &optional (top-reduction-only $poly_top_reduction_only))
    555   "An implementation of the Buchberger algorithm. Return Grobner basis
    556 of the ideal generated by the polynomial list F.  Polynomials 0 to
    557 START-1 are assumed to be a Grobner basis already, so that certain
    558 critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
    559 reduction will be preformed. This function assumes that all polynomials
    560 in F are non-zero."
    561   (declare (type fixnum start))
    562   (when (endp f) (return-from buchberger f)) ;cut startup costs
    563   (debug-cgb "~&GROBNER BASIS - BUCHBERGER ALGORITHM")
    564   (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
    565   #+grobner-check  (when (plusp start)
    566                      (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
    567   ;;Initialize critical pairs
    568   (let ((b (pair-queue-initialize (make-pair-queue)
    569                                   f start))
    570         (b-done (make-hash-table :test #'equal)))
    571     (declare (type priority-queue b) (type hash-table b-done))
    572     (dotimes (i (1- start))
    573       (do ((j (1+ i) (1+ j))) ((>= j start))
    574         (setf (gethash (list (elt f i) (elt f j)) b-done) t)))
    575     (do ()
    576         ((pair-queue-empty-p b)
    577          #+grobner-check(grobner-test ring f f)
    578          (debug-cgb "~&GROBNER END")
    579          f)
    580       (let ((pair (pair-queue-remove b)))
    581         (declare (type pair pair))
    582         (cond
    583           ((criterion-1 pair) nil)
    584           ((criterion-2 pair b-done f) nil)
    585           (t
    586            (let ((sp (normal-form ring (spoly ring (pair-first pair)
    587                                               (pair-second pair))
    588                                   f top-reduction-only)))
    589              (declare (type poly sp))
    590              (cond
    591                ((poly-zerop sp)
    592                 nil)
    593                (t
    594                 (setf sp (poly-primitive-part ring sp)
    595                       f (nconc f (list sp)))
    596                 ;; Add new critical pairs
    597                 (dolist (h f)
    598                   (pair-queue-insert b (make-pair h sp)))
    599                 (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
    600                            (pair-sugar pair) (length f) (pair-queue-size b)
    601                            (hash-table-count b-done)))))))
    602         (setf (gethash (list (pair-first pair) (pair-second pair)) b-done)
    603               t)))))
    604 
    605 (defun parallel-buchberger (ring f start &optional (top-reduction-only $poly_top_reduction_only))
    606   "An implementation of the Buchberger algorithm. Return Grobner basis
    607 of the ideal generated by the polynomial list F.  Polynomials 0 to
    608 START-1 are assumed to be a Grobner basis already, so that certain
    609 critical pairs will not be examined. If TOP-REDUCTION-ONLY set, top
    610 reduction will be preformed."
    611   (declare (ignore top-reduction-only)
    612            (type fixnum start))
    613   (when (endp f) (return-from parallel-buchberger f)) ;cut startup costs
    614   (debug-cgb "~&GROBNER BASIS - PARALLEL-BUCHBERGER ALGORITHM")
    615   (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
    616   #+grobner-check  (when (plusp start)
    617                      (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
    618   ;;Initialize critical pairs
    619   (let ((b (pair-queue-initialize (make-pair-queue) f start))
    620         (b-done (make-hash-table :test #'equal)))
    621     (declare (type priority-queue b)
    622              (type hash-table b-done))
    623     (dotimes (i (1- start))
    624       (do ((j (1+ i) (1+ j))) ((>= j start))
    625         (declare (type fixnum j))
    626         (setf (gethash (list (elt f i) (elt f j)) b-done) t)))
    627     (do ()
    628         ((pair-queue-empty-p b)
    629          #+grobner-check(grobner-test ring f f)
    630          (debug-cgb "~&GROBNER END")
    631          f)
    632       (let ((pair (pair-queue-remove b)))
    633         (when (null (pair-division-data pair))
    634           (setf (pair-division-data pair) (list (spoly ring
    635                                                        (pair-first pair)
    636                                                        (pair-second pair))
    637                                                 (make-poly-zero)
    638                                                 (funcall (ring-unit ring))
    639                                                 0)))
    640         (cond
    641           ((criterion-1 pair) nil)
    642           ((criterion-2 pair b-done f) nil)
    643           (t
    644            (let* ((dd (pair-division-data pair))
    645                   (p (first dd))
    646                   (sp (second dd))
    647                   (c (third dd))
    648                   (division-count (fourth dd)))
    649              (cond
    650                ((poly-zerop p)          ;normal form completed
    651                 (debug-cgb "~&~3T~d reduction~:p" division-count)
    652                 (cond
    653                   ((poly-zerop sp)
    654                    (debug-cgb " ---> 0")
    655                    nil)
    656                   (t
    657                    (setf sp (poly-nreverse sp)
    658                          sp (poly-primitive-part ring sp)
    659                          f (nconc f (list sp)))
    660                    ;; Add new critical pairs
    661                    (dolist (h f)
    662                      (pair-queue-insert b (make-pair h sp)))
    663                    (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d; Pairs done: ~d;"
    664                               (pair-sugar pair) (length f) (pair-queue-size b)
    665                               (hash-table-count b-done))))
    666                 (setf (gethash (list (pair-first pair) (pair-second pair))
    667                                b-done) t))
    668                (t                               ;normal form not complete
    669                 (do ()
    670                     ((cond
    671                        ((> (poly-sugar sp) (pair-sugar pair))
    672                         (debug-cgb "(~a)?" (poly-sugar sp))
    673                         t)
    674                        ((poly-zerop p)
    675                         (debug-cgb ".")
    676                         t)
    677                        (t nil))
    678                      (setf (first dd) p
    679                            (second dd) sp
    680                            (third dd) c
    681                            (fourth dd) division-count
    682                            (pair-sugar pair) (poly-sugar sp))
    683                      (pair-queue-insert b pair))
    684                   (multiple-value-setq (p sp c division-count)
    685                     (normal-form-step ring f p sp c division-count))))))))))))
    686 
    687 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    688 ;;
    689 ;; Grobner Criteria
    690 ;;
    691 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    692 
    693 (defun criterion-1 (pair)
    694   "Returns T if the leading monomials of the two polynomials
    695 in G pointed to by the integers in PAIR have disjoint (relatively prime)
    696 monomials. This test is known as the first Buchberger criterion."
    697   (declare (type pair pair))
    698   (let ((f (pair-first pair))
    699         (g (pair-second pair)))
    700     (when (monom-rel-prime-p (poly-lm f) (poly-lm g))
    701       (debug-cgb ":1")
    702       (return-from criterion-1 t))))
    703 
    704 (defun criterion-2 (pair b-done partial-basis
    705                     &aux (f (pair-first pair)) (g (pair-second pair))
    706                          (place :before))
    707   "Returns T if the leading monomial of some element P of
    708 PARTIAL-BASIS divides the LCM of the leading monomials of the two
    709 polynomials in the polynomial list PARTIAL-BASIS, and P paired with
    710 each of the polynomials pointed to by the the PAIR has already been
    711 treated, as indicated by the absence in the hash table B-done."
    712   (declare (type pair pair) (type hash-table b-done)
    713            (type poly f g))
    714   ;; In the code below we assume that pairs are ordered as follows:
    715   ;; if PAIR is (I J) then I appears before J in the PARTIAL-BASIS.
    716   ;; We traverse the list PARTIAL-BASIS and keep track of where we
    717   ;; are, so that we can produce the pairs in the correct order
    718   ;; when we check whether they have been processed, i.e they
    719   ;; appear in the hash table B-done
    720   (dolist (h partial-basis nil)
    721     (cond
    722      ((eq h f)
    723       #+grobner-check(assert (eq place :before))
    724       (setf place :in-the-middle))
    725      ((eq h g)
    726       #+grobner-check(assert (eq place :in-the-middle))
    727       (setf place :after))
    728      ((and (monom-divides-monom-lcm-p (poly-lm h) (poly-lm f) (poly-lm g))
    729            (gethash (case place
    730                       (:before (list h f))
    731                       ((:in-the-middle :after) (list f h)))
    732                     b-done)
    733            (gethash (case place
    734                       ((:before :in-the-middle) (list h g))
    735                       (:after (list g h)))
    736                     b-done))
    737       (debug-cgb ":2")
    738       (return-from criterion-2 t)))))
    739 
    740 
    741 
    742 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    743 ;;
    744 ;; An implementation of the algorithm of Gebauer and Moeller, as
    745 ;; described in the book of Becker-Weispfenning, p. 232
    746 ;;
    747 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    748 
    749 (defun gebauer-moeller (ring f start &optional (top-reduction-only $poly_top_reduction_only))
    750   "Compute Grobner basis by using the algorithm of Gebauer and
    751 Moeller.  This algorithm is described as BUCHBERGERNEW2 in the book by
    752 Becker-Weispfenning entitled ``Grobner Bases''. This function assumes
    753 that all polynomials in F are non-zero."
    754   (declare (ignore top-reduction-only)
    755            (type fixnum start))
    756   (cond
    757    ((endp f) (return-from gebauer-moeller nil))
    758    ((endp (cdr f))
    759     (return-from gebauer-moeller (list (poly-primitive-part ring (car f))))))
    760    (debug-cgb "~&GROBNER BASIS - GEBAUER MOELLER ALGORITHM")
    761    (when (plusp start) (debug-cgb "~&INCREMENTAL:~d done" start))
    762   #+grobner-check  (when (plusp start)
    763                      (grobner-test ring (subseq f 0 start) (subseq f 0 start)))
    764   (let ((b (make-pair-queue))
    765         (g (subseq f 0 start))
    766         (f1 (subseq f start)))
    767     (do () ((endp f1))
    768       (multiple-value-setq (g b)
    769         (gebauer-moeller-update g b (poly-primitive-part ring (pop f1)))))
    770     (do () ((pair-queue-empty-p b))
    771       (let* ((pair (pair-queue-remove b))
    772              (g1 (pair-first pair))
    773              (g2 (pair-second pair))
    774              (h (normal-form ring (spoly ring g1 g2)
    775                              g
    776                              nil #| Always fully reduce! |#
    777                              )))
    778         (unless (poly-zerop h)
    779           (setf h (poly-primitive-part ring h))
    780           (multiple-value-setq (g b)
    781             (gebauer-moeller-update g b h))
    782           (debug-cgb "~&Sugar: ~d Polynomials: ~d; Pairs left: ~d~%"
    783                      (pair-sugar pair) (length g) (pair-queue-size b))
    784           )))
    785     #+grobner-check(grobner-test ring g f)
    786     (debug-cgb "~&GROBNER END")
    787     g))
    788 
    789 (defun gebauer-moeller-update (g b h
    790                  &aux
    791                  c d e
    792                  (b-new (make-pair-queue))
    793                  g-new)
    794   "An implementation of the auxillary UPDATE algorithm used by the
    795 Gebauer-Moeller algorithm. G is a list of polynomials, B is a list of
    796 critical pairs and H is a new polynomial which possibly will be added
    797 to G. The naming conventions used are very close to the one used in
    798 the book of Becker-Weispfenning."
    799   (declare
    800    #+allegro (dynamic-extent b)
    801    (type poly h)
    802    (type priority-queue b))
    803   (setf c g d nil)
    804   (do () ((endp c))
    805     (let ((g1 (pop c)))
    806       (declare (type poly g1))
    807       (when (or (monom-rel-prime-p (poly-lm h) (poly-lm g1))
    808                 (and
    809                  (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
    810                                          (poly-lm h) (poly-lm g2)
    811                                          (poly-lm h) (poly-lm g1)))
    812                          c)
    813                  (notany #'(lambda (g2) (monom-lcm-divides-monom-lcm-p
    814                                          (poly-lm h) (poly-lm g2)
    815                                          (poly-lm h) (poly-lm g1)))
    816                          d)))
    817         (push g1 d))))
    818   (setf e nil)
    819   (do () ((endp d))
    820     (let ((g1 (pop d)))
    821       (declare (type poly g1))
    822       (unless (monom-rel-prime-p (poly-lm h) (poly-lm g1))
    823         (push g1 e))))
    824   (do () ((pair-queue-empty-p b))
    825     (let* ((pair (pair-queue-remove b))
    826            (g1 (pair-first pair))
    827            (g2 (pair-second pair)))
    828       (declare (type pair pair)
    829                (type poly g1 g2))
    830       (when (or (not (monom-divides-monom-lcm-p
    831                       (poly-lm h)
    832                       (poly-lm g1) (poly-lm g2)))
    833                 (monom-lcm-equal-monom-lcm-p
    834                  (poly-lm g1) (poly-lm h)
    835                  (poly-lm g1) (poly-lm g2))
    836                 (monom-lcm-equal-monom-lcm-p
    837                  (poly-lm h) (poly-lm g2)
    838                  (poly-lm g1) (poly-lm g2)))
    839         (pair-queue-insert b-new (make-pair g1 g2)))))
    840   (dolist (g3 e)
    841     (pair-queue-insert b-new (make-pair h g3)))
    842   (setf g-new nil)
    843   (do () ((endp g))
    844     (let ((g1 (pop g)))
    845       (declare (type poly g1))
    846       (unless (monom-divides-p (poly-lm h) (poly-lm g1))
    847         (push g1 g-new))))
    848   (push h g-new)
    849   (values g-new b-new))
    850 
    851 
    852 
    853 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    854 ;;
    855 ;; Standard postprocessing of Grobner bases
    856 ;;
    857 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    858 
    859 (defun reduction (ring plist)
    860   "Reduce a list of polynomials PLIST, so that non of the terms in any of
    861 the polynomials is divisible by a leading monomial of another
    862 polynomial.  Return the reduced list."
    863   (do ((q plist)
    864        (found t))
    865       ((not found)
    866        (mapcar #'(lambda (x) (poly-primitive-part ring x)) q))
    867     ;;Find p in Q such that p is reducible mod Q\{p}
    868     (setf found nil)
    869     (dolist (x q)
    870       (let ((q1 (remove x q)))
    871         (multiple-value-bind (h c div-count)
    872             (normal-form ring x q1 nil #| not a top reduction! |# )
    873           (declare (ignore c))
    874           (unless (zerop div-count)
    875             (setf found t q q1)
    876             (unless (poly-zerop h)
    877               (setf q (nconc q1 (list h))))
    878             (return)))))))
    879 
    880 (defun minimization (p)
    881   "Returns a sublist of the polynomial list P spanning the same
    882 monomial ideal as P but minimal, i.e. no leading monomial
    883 of a polynomial in the sublist divides the leading monomial
    884 of another polynomial."
    885   (do ((q p)
    886        (found t))
    887       ((not found) q)
    888     ;;Find p in Q such that lm(p) is in LM(Q\{p})
    889     (setf found nil
    890           q (dolist (x q q)
    891               (let ((q1 (remove x q)))
    892                 (when (member-if #'(lambda (p) (monom-divides-p (poly-lm x) (poly-lm p))) q1)
    893                   (setf found t)
    894                   (return q1)))))))
    895 
    896 (defun poly-normalize (ring p &aux (c (poly-lc p)))
    897   "Divide a polynomial by its leading coefficient. It assumes
    898 that the division is possible, which may not always be the
    899 case in rings which are not fields. The exact division operator
    900 is assumed to be provided by the RING structure of the
    901 COEFFICIENT-RING package."
    902   (mapc #'(lambda (term)
    903             (setf (term-coeff term) (funcall (ring-div ring) (term-coeff term) c)))
    904         (poly-termlist p))
    905   p)
    906 
    907 (defun poly-normalize-list (ring plist)
    908   "Divide every polynomial in a list PLIST by its leading coefficient. "
    909   (mapcar #'(lambda (x) (poly-normalize ring x)) plist))
    910 
    911 
    912 
    913 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    914 ;;
    915 ;; Algorithm and Pair heuristic selection
    916 ;;
    917 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    918324
    919325(defun find-grobner-function (algorithm)
Note: See TracChangeset for help on using the changeset viewer.