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 39 for branches/f4grobner


Ignore:
Timestamp:
2015-06-02T00:11:37-07:00 (9 years ago)
Author:
Marek Rychlik
Message:

* empty log message *

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/f4grobner/f4-matrix.lisp

    r36 r39  
    2424(macsyma-module f4-maxima)
    2525
    26 (defun poly-monomials (p)
    27   (mapcar #'term-monom (poly-termlist p)))
     26;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
     27;;
     28;; An implementation of the normal form
     29;;
     30;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    2831
    29 (defun make-f4-matrix (f fl &aux ml))
    30  
    31  
     32(defun normal-form-step (ring fl p r c division-count
     33                         &aux (g (find (poly-lm p) fl
     34                                       :test #'monom-divisible-by-p
     35                                       :key #'poly-lm)))
     36  (cond
     37   (g                                   ;division possible
     38    (incf division-count)
     39    (multiple-value-bind (gcd cg cp)
     40        (funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
     41      (declare (ignore gcd))
     42      (let ((m (monom-div (poly-lm p) (poly-lm g))))
     43        ;; Multiply the equation c*f=sum ai*fi+r+p by cg.
     44        (setf r (scalar-times-poly ring cg r)
     45              c (funcall (ring-mul ring) c cg)
     46              ;; p := cg*p-cp*m*g
     47              p (grobner-op ring cp cg m p g))))
     48    (debug-cgb "/"))
     49   (t                                                   ;no division possible
     50    (push (poly-lt p) (poly-termlist r))                ;move lt(p) to remainder
     51    (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
     52    (pop (poly-termlist p))                             ;remove lt(p) from p
     53    (debug-cgb "+")))
     54  (values p r c division-count))
    3255
     56;; Merge it sometime with poly-pseudo-divide
     57(defun normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
     58  ;; Loop invariant: c*f0=sum ai*fi+r+f, where f0 is the initial value of f
     59  #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
     60  (do ((r (make-poly-zero))
     61       (c (funcall (ring-unit ring)))
     62       (division-count 0))
     63      ((or (poly-zerop f)
     64           ;;(endp fl)
     65           (and top-reduction-only (not (poly-zerop r))))
     66       (progn
     67         (debug-cgb "~&~3T~d reduction~:p" division-count)
     68         (when (poly-zerop r)
     69           (debug-cgb " ---> 0")))
     70       (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
     71       (values f c division-count))
     72    (declare (fixnum division-count)
     73             (type poly r))
     74    (multiple-value-setq (f r c division-count)
     75      (normal-form-step ring fl f r c division-count))))
Note: See TracChangeset for help on using the changeset viewer.