;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                              
;;;  Copyright (C) 2015 Marek Rychlik <rychlik@u.arizona.edu>		 
;;;  		       								 
;;;  This program is free software; you can redistribute it and/or modify	 
;;;  it under the terms of the GNU General Public License as published by	 
;;;  the Free Software Foundation; either version 2 of the License, or		 
;;;  (at your option) any later version.					 
;;; 		       								 
;;;  This program is distributed in the hope that it will be useful,		 
;;;  but WITHOUT ANY WARRANTY; without even the implied warranty of		 
;;;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the		 
;;;  GNU General Public License for more details.				 
;;; 		       								 
;;;  You should have received a copy of the GNU General Public License		 
;;;  along with this program; if not, write to the Free Software 		 
;;;  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.	 
;;;										 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package :maxima)

(macsyma-module f4-maxima)

(defun f4-grobner-op (ring c1 c2 m f g)
  "Returns C2*F-C1*M*G, where F and G are polynomials M is a monomial.
Assume that the leading terms will cancel."
  #+grobner-check(funcall (ring-zerop ring)
			  (funcall (ring-sub ring)
				   (funcall (ring-mul ring) c2 (poly-lc f))
				   (funcall (ring-mul ring) c1 (poly-lc g))))
  #+grobner-check(monom-equal-p (poly-lm f) (monom-mul m (poly-lm g)))
  (poly-sub ring
	    (scalar-times-poly ring c2 f)
	    (scalar-times-poly ring c1 (monom-times-poly m g))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; An implementation of the normal form
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defun f4-normal-form-step (ring fl p r c division-count
			 &aux (g (find (poly-lm p) fl
				       :test #'monom-divisible-by-p
				       :key #'poly-lm)))
  (cond
   (g					;division possible
    (incf division-count)
    (multiple-value-bind (gcd cg cp)
	(funcall (ring-ezgcd ring) (poly-lc g) (poly-lc p))
      (declare (ignore gcd))
      (let ((m (monom-div (poly-lm p) (poly-lm g))))
	;; Multiply the equation c*f=sum ai*fi+r+p by cg.
	(setf r (scalar-times-poly ring cg r)
	      c (funcall (ring-mul ring) c cg)
              ;; p := cg*p-cp*m*g
	      p (grobner-op ring cp cg m p g))))
    (debug-cgb "/"))
   (t							;no division possible
    (push (poly-lt p) (poly-termlist r))		;move lt(p) to remainder
    (setf (poly-sugar r) (max (poly-sugar r) (term-sugar (poly-lt p))))
    (pop (poly-termlist p))				;remove lt(p) from p
    (debug-cgb "+")))
  (values p r c division-count))

;; Merge it sometime with poly-pseudo-divide
(defun f4-normal-form (ring f fl &optional (top-reduction-only $poly_top_reduction_only))
  ;; Loop invariant: c*f0=sum ai*fi+r+f, where f0 is the initial value of f
  #+grobner-check(when (null fl) (warn "normal-form: empty divisor list."))
  (do ((r (make-poly-zero))
       (c (funcall (ring-unit ring)))
       (division-count 0))
      ((or (poly-zerop f)
	   ;;(endp fl)
	   (and top-reduction-only (not (poly-zerop r))))
       (progn
	 (debug-cgb "~&~3T~d reduction~:p" division-count)
	 (when (poly-zerop r)
	   (debug-cgb " ---> 0")))
       (setf (poly-termlist f) (nreconc (poly-termlist r) (poly-termlist f)))
       (values f c division-count))
    (declare (fixnum division-count)
	     (type poly r))
    (multiple-value-setq (f r c division-count)
      (f4-normal-form-step ring fl f r c division-count))))