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


Ignore:
Timestamp:
2016-05-31T17:29:55-07:00 (9 years ago)
Author:
Marek Rychlik
Message:

* empty log message *

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/f4grobner/division.lisp

    r4048 r4049  
    2121
    2222(defpackage "DIVISION"
    23   (:use :cl :utils :ring :monom :polynomial :grobner-debug :term :ring-and-order)
     23  (:use :cl :utils :monom :polynomial :grobner-debug)
    2424  (:export "$POLY_TOP_REDUCTION_ONLY"
    2525           "POLY-PSEUDO-DIVIDE"
    … …  
    5555                   (multiply c2 (leading-coefficient f))
    5656                   (multiply c1 (leading-coefficient g))))
    57   #+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
     57  #+grobner-check(universal-equalp (leading-monomial f) (multiply m (leading-monomial g)))
    5858  ;; Note that below we can drop the leading terms of f ang g for the
    5959  ;; purpose of polynomial arithmetic. 
    … …  
    6161  ;; TODO: Make sure that the sugar calculation is correct if leading
    6262  ;; terms are dropped.
    63   (poly-sub ring-and-order
    64             (scalar-times-poly-1 ring c2 f)
    65             (scalar-times-poly-1 ring c1 (monom-times-poly m g))))
    66 
    67 (defun check-loop-invariant (ring-and-order c f a fl r p
     63  (subtract
     64   (multiply c2 f)
     65   (multiply c1 (multiply m g))))
     66
     67(defun check-loop-invariant (c f a fl r p
    6868                             &aux
    69                                (ring (ro-ring ring-and-order))
    7069                               (p-zero (make-poly-zero))
    7170                               (a (mapcar #'poly-reverse a))
    … …  
    8483          c f a fl r p)
    8584  |#
    86   (flet ((p-add (x y) (poly-add ring-and-order x y))
    87          (p-sub (x y) (poly-sub ring-and-order x y))
    88          (p-mul (x y) (poly-mul ring-and-order x y)))
    89     (let* ((prod (inner-product a fl p-add p-mul p-zero))
    90            (succeeded-p
    91             (poly-zerop
    92              (p-sub
    93               (scalar-times-poly ring c f)
    94               (reduce #'p-add (list prod r p))))))
    95       (unless succeeded-p
    96           (error "#### Polynomial division Loop invariant failed ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%"
    97                  c f a fl r p))
    98       succeeded-p)))
     85  (let* ((prod (inner-product a fl #'add #'multiply 0))
     86         (succeeded-p
     87          (universal-zerop
     88           (subtract
     89            (multiply c f)
     90            (reduce #'add (list prod r p))))))
     91    (unless succeeded-p
     92      (error "#### Polynomial division Loop invariant failed ####:~%C=~A~%F=~A~%A=~A~%FL=~A~%R=~A~%P=~A~%"
     93             c f a fl r p))
     94    succeeded-p))
    9995 
    10096
    101 (defun poly-pseudo-divide (ring-and-order f fl
    102                            &aux
    103                              (ring (ro-ring ring-and-order)))
     97(defun poly-pseudo-divide (f fl)
    10498  "Pseudo-divide a polynomial F by the list of polynomials FL. Return
    10599multiple values. The first value is a list of quotients A.  The second
    … …  
    112106  (declare (type poly f) (list fl))
    113107  ;; Loop invariant: c*f=sum ai*fi+r+p, where p must eventually become 0
    114   (do ((r (make-poly-zero))
    115        (c (funcall (ring-unit ring)))
    116        (a (make-list (length fl) :initial-element (make-poly-zero)))
     108  (do ((r 0)
     109       (c 1)
     110       (a (make-list (length fl) :initial-element 0))
    117111       (division-count 0)
    118112       (p f))
    119       ((poly-zerop p)
    120        #+grobner-check(check-loop-invariant ring-and-order c f a fl r p)
     113      ((universal-zerop p)
     114       #+grobner-check(check-loop-invariant c f a fl r p)
    121115       (debug-cgb "~&~3T~d reduction~:p" division-count)
    122        (when (poly-zerop r) (debug-cgb " ---> 0"))
     116       (when (universal-zerop r) (debug-cgb " ---> 0"))
    123117       ;; We obtained the terms in reverse order, so must fix that
    124118       (setf a (mapcar #'poly-nreverse a)
    125119             r (poly-nreverse r))
    126120       ;; Initialize the sugar of the quotients
    127        (mapc #'poly-reset-sugar a)
     121       ;; (mapc #'poly-reset-sugar a) ;; TODO: Sugar is currently unimplemented
    128122       (values a r c division-count))
    129123    (declare (fixnum division-count))
    130124    ;; Check the loop invariant here
    131     #+grobner-check(check-loop-invariant ring-and-order c f a fl r p)
     125    #+grobner-check(check-loop-invariant c f a fl r p)
    132126    (do ((fl fl (rest fl))              ;scan list of divisors
    133127         (b a (rest b)))
    134128        ((cond
    135129           ((endp fl)                           ;no division occurred
    136             (push (poly-lt p) (poly-termlist r)) ;move lt(p) to remainder
    137             (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
     130            (push (leading-term p) (poly-termlist r)) ;move lt(p) to remainder
     131            ;;(setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
    138132            (pop (poly-termlist p))     ;remove lt(p) from p
    139133            t)
    140            ((monom-divides-p (poly-lm (car fl)) (poly-lm p)) ;division occurred
     134           ((monom-divides-p (leading-monomial (car fl)) (leading-monomial p)) ;division occurred
    141135            (incf division-count)
    142136            (multiple-value-bind (gcd c1 c2)
    143                 (funcall (ring-ezgcd ring) (poly-lc (car fl)) (poly-lc p))
     137                (universal-ezgcd (leading-coefficient (car fl)) (leading-coefficient p))
    144138              (declare (ignore gcd))
    145               (let ((m (monom-div (poly-lm p) (poly-lm (car fl)))))
     139              (let ((m (divide (leading-monomial p) (leading-monomial (car fl)))))
    146140                ;; Multiply the equation c*f=sum ai*fi+r+p by c1.
    147141                (mapl #'(lambda (x)
    148                           (setf (car x) (scalar-times-poly ring c1 (car x))))
     142                          (setf (car x) (multiply c1 (car x))))
    149143                      a)
    150                 (setf r (scalar-times-poly ring c1 r)
    151                       c (funcall (ring-mul ring) c c1)
    152                       p (grobner-op ring-and-order c2 c1 m p (car fl)))
     144                (setf r (multiply c1 r)
     145                      c (multiply c c1)
     146                      p (grobner-op c2 c1 m p (car fl)))
    153147                (push (make-term :monom m :coeff c2) (poly-termlist (car b))))
    154148              t))))
    155149      )))
    156150
    157 (defun poly-exact-divide (ring-and-order f g)
     151(defun poly-exact-divide (f g)
    158152  "Divide a polynomial F by another polynomial G. Assume that exact division
    159153with no remainder is possible. Returns the quotient."
    160   (declare (type poly f g) (type ring-and-order ring-and-order))
     154  (declare (type poly f g))
    161155  (multiple-value-bind (quot rem coeff division-count)
    162       (poly-pseudo-divide ring-and-order f (list g))
     156      (poly-pseudo-divide f (list g))
    163157    (declare (ignore division-count coeff)
    164158             (list quot)
    165159             (type poly rem)
    166160             (type fixnum division-count))
    167     (unless (poly-zerop rem) (error "Exact division failed."))
     161    (unless (universal-zerop rem) (error "Exact division failed."))
    168162    (car quot)))
    169163
    … …  
    174168;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    175169
    176 (defun normal-form-step (ring-and-order fl p r c division-count
     170(defun normal-form-step (fl p r c division-count
    177171                         &aux
    178                            (ring (ro-ring ring-and-order))
    179                            (g (find (poly-lm p) fl
     172                           (g (find (leading-monomial p) fl
    180173                                    :test #'monom-divisible-by-p
    181                                     :key #'poly-lm)))
     174                                    :key #'leading-monomial)))
    182175  (cond
    183176   (g                                   ;division possible
    184177    (incf division-count)
    185178    (multiple-value-bind (gcd cg cp)
    186         (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
     179        (universal-ezgcd (leading-coefficient g) (leading-coefficient p))
    187180      (declare (ignore gcd))
    188       (let ((m (monom-div (poly-lm p) (poly-lm g))))
     181      (let ((m (divide (leading-monomial p) (leading-monomial g))))
    189182        ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
    190         (setf r (scalar-times-poly ring cg r)
    191               c (funcall (ring-mul ring) c cg)
     183        (setf r (multiply cg r)
     184              c (multiply c cg)
    192185              ;; p := cg*p-cp*m*g
    193               p (grobner-op ring-and-order cp cg m p g))))
     186              p (grobner-op cp cg m p g))))
    194187    (debug-cgb "/"))
    195188   (t                                                   ;no division possible
    196     (push (poly-lt p) (poly-termlist r))                ;move lt(p) to remainder
    197     (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
     189    (push (leading-term p) (poly-termlist r))           ;move lt(p) to remainder
     190    ;;(setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
    198191    (pop (poly-termlist p))                             ;remove lt(p) from p
    199192    (debug-cgb "+")))
    … …  
    207200;; method would be to test NORMAL-FORM using POLY-PSEUDO-DIVIDE.
    208201;;
    209 (defun normal-form (ring-and-order f fl
    210                     &optional
    211                       (top-reduction-only $poly_top_reduction_only)
    212                       (ring (ro-ring ring-and-order)))
     202(defun normal-form (f fl
     203                    &optional
     204                      (top-reduction-only $poly_top_reduction_only))
    213205  #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
    214   (do ((r (make-poly-zero))
    215        (c (funcall (ring-unit ring)))
     206  (do ((r 0)
     207       (c 1)
    216208       (division-count 0))
    217       ((or (poly-zerop f)
     209      ((or (universal-zerop f)
    218210           ;;(endp fl)
    219            (and top-reduction-only (not (poly-zerop r))))
     211           (and top-reduction-only (not (universal-zerop r))))
    220212       (progn
    221213         (debug-cgb "~&~3T~D reduction~:P" division-count)
    222          (when (poly-zerop r)
     214         (when (universal-zerop r)
    223215           (debug-cgb " ---> 0")))
    224216       (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
    … …  
    227219             (type poly r))
    228220    (multiple-value-setq (f r c division-count)
    229       (normal-form-step ring-and-order fl f r c division-count))))
     221      (normal-form-step fl f r c division-count))))
    230222
    231223(defun buchberger-criterion (ring-and-order g)
Note: See TracChangeset for help on using the changeset viewer.