source: branches/f4grobner/5am-division.lisp@ 4445

;;; -*-  Mode: Lisp -*- 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;                                                                              
;;;  Copyright (C) 1999, 2002, 2009, 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.  
;;;                                                                              
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Run tests using 5am unit testing framework
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; We assume that QuickLisp package manager is installed.
;; See :
;;      https://www.quicklisp.org/beta/
;;

;; The following is unnecessary after running:
;; * (ql:add-to-init-file)
;; at lisp prompt:
;;(load "~/quicklisp/setup")

(ql:quickload :fiveam)

(defpackage #:5am-division
  (:use :cl :it.bese.fiveam :copy :monom :polynomial :infix :symbolic-polynomial :division :ring))

(in-package :5am-division)

(def-suite division-suite 
    :description "Division algorithm suite")

(in-suite division-suite)

;; Manual calculation supporting the test below.
;; We divide X^2 by [X+Y,X-2*Y] with LEX> as order.
;; LM(X^2)=X^2 is divisible by LM(X+Y)=X so the first partial quotient is X.
;; Next, X^2 - X*(X+Y) = -X*Y.
;; LM(-X*Y)=X*Y is divibile by LM(X+Y)=X so the second partial quotient is -Y.
;; Next, -X*Y-(-Y)*(X+Y) = Y^2.
;; LM(Y^2)=Y^2 is not divisible by LM(X+Y)=X or LM(X-2*Y)=X. Hence, division
;; ends. The list of quotients is [X-Y,0]. The remainder is Y^2

(def-fixture division-context ()
  (let* ((f (string->poly "x^2" '(x y)))
         (y-sq (string->poly "y^2" '(x y)))
         (fl (cdr (string->poly "[x+y,x-2*y]" '(x y))))
         (quotients (cdr (string->poly "[x-y,0]" '(x y))))
         (one (make-instance 'rational-field :value 1)))
    (&body)))

(test normal-form
  "Normal form"
  (with-fixture division-context ()
    (is (universal-equalp (multiple-value-list (normal-form f fl)) (list y-sq one 2)))
    (is (universal-equalp (multiple-value-list (poly-pseudo-divide f fl)) (list quotients y-sq one 2)))
    (is-false (buchberger-criterion fl))
    )
  )

;; Maxima
;;poly_normal_form(3*x^3*y+2*z, [x^2*y+x,x-y^2-z^3],[x,y,z]);
;;Result: (-3*z^6)-6*y^2*z^3+2*z-3*y^4
(test normal-form-simple 
  (let ((vars '(x y z)))
    (is (universal-equalp (normal-form (string->poly "3*x^3*y+2*z" vars)
                                       (cdr (string->poly "[x^2*y+x,x-y^2-z^3]" vars)))
                          (string->poly "(-3*z^6)-6*y^2*z^3+2*z-3*y^4" vars)))))

(test normal-form-easy
  "Easy normal form tests"
  (is (universal-zerop (normal-form (string->poly "0" '(x y)) (cdr (string->poly "[x,y]" '(x y))))))
  ;; Maxima equivalent: poly_normal_form(3*x^2*y-x*y-1,[x-y,x+y],[x,y]);
  (is (universal-equalp (normal-form (string->poly "3*x^2*y-x*y-1" '(x y)) (cdr (string->poly "[x-y,x+y]" '(x y))))
                        (string->poly "3*y^3-y^2-1" '(x y))))
  ;; Maxima equivalent: poly_normal_form(3*x^2*y*z-x*y^3-1,[x^2-2*y,x+y*z],[x,y,z]);
  (is (universal-equalp (normal-form (string->poly "3*x^2*y*z-x*y^3-1" '(x y z))
                                     (cdr (string->poly "[x^2-2*y,x+y*z]" '(x y z))))
                        (string->poly "y^4*z+6*y^2*z-1" '(x y z)))))

(def-fixture exact-division-context ()
  (let* ((f (string->poly "x^2-4*y^2" '(x y)))
         (g (string->poly "x+2*y" '(x y)))
         (h (string->poly "x-2*y" '(x y))))
    (&body)))

(test exact-division 
  "Exact division in polynomial ring"
  (with-fixture exact-division-context ()
    (is (universal-equalp (poly-exact-divide f g) h))
    (is (universal-zerop (subtract-from (poly-exact-divide f g) h)))))


;; Check if a set of generators satisfies the Buchberger criterion
;; The example is the Enneper surface ideal. Run this in Maxima to
;; obtain the Grobner basis:
;;   poly_grobner([x-3*u-3*u*v^2+u^3,y-3*v-3*u^2*v+v^3,z-3*u^2+3*v^2],[u,v,x,y,z]);
(def-fixture buchberger-criterion-context ()
  (let ((fl (cdr (string->poly "[x-3*u-3*u*v^2+u^3,y-3*v-3*u^2*v+v^3,z-3*u^2+3*v^2]" '(u v x y z))))
        (gb (cdr (string->poly "[x-3*u*v^2+u^3-3*u,y+v^3-3*u^2*v-3*v,z+3*v^2-3*u^2,
       (-u*z)-3*x+6*u*v^2+9*u,(-v*z)+y-2*v^3-3*v,z^2+6*v^2*z-9*z-9*v*y+9*u*x,
       4*u*v*z-3*u*y+3*v*x,2*u*z^2+6*x*z-18*u*z-9*u*v*y+9*v^2*x,
       (-8*u*z^3)-24*x*z^2+72*u*z^2-36*v^2*x*z+27*u*y^2-27*v*x*y,
       z^3+18*v^2*z^2-18*z^2-54*v*y*z+54*v^2*z+81*z+27*y^2-27*x^2,
       (-4*z^4)+48*z^3-108*v*y*z^2+108*z^2+135*y^2*z+324*v*y*z+108*x^2*z
               -1296*v^2*z-1944*z-243*v^2*y^2-648*y^2+243*v^2*x^2+648*x^2,
       8*v*z^3-9*y*z^2+72*v*z^2+54*v^2*y*z-27*y*z-27*v*y^2+27*v*x^2,
       (-8*v*z^4)+12*y*z^3-96*v*z^3-216*v*z^2-135*v*y^2*z+324*y*z-27*v*x^2*z
                 +81*y^3+81*v*y^2-81*x^2*y-81*v*x^2,
       (-64*v*z^6)+120*y*z^5-1152*v*z^5+288*y*z^4-5184*v*z^4-648*v*y^2*z^3
                  -216*y*z^3+6912*v*z^3+81*y^3*z^2-9720*v*y^2*z^2
                  -1539*x^2*y*z^2+31104*y*z^2+62208*v*z^2+8505*y^3*z
                  +46656*v*y^2*z-8505*x^2*y*z-93312*y*z+729*v*y^4-23328*y^3
                  -1458*v*x^2*y^2-23328*v*y^2+23328*x^2*y+729*v*x^4
                  +23328*v*x^2,
       8*z^6-72*z^5+648*v*y*z^4-648*z^4-945*y^2*z^3+5184*v*y*z^3-189*x^2*z^3
            +5832*z^3+972*y^2*z^2+17496*v*y*z^2-2430*x^2*z^2+8748*v*y^3*z
            -19683*y^2*z+2187*x^2*z-5103*y^4-4374*v*y^3+5832*x^2*y^2
            +4374*v*x^2*y-729*x^4,
       8*z^7-48*z^6+648*v*y*z^5-864*z^5-945*y^2*z^4+5832*v*y*z^4-189*x^2*z^4
            +3888*z^4+81*y^2*z^3+17496*v*y*z^3-2997*x^2*z^3+17496*z^3
            +8748*v*y^3*z^2-16767*y^2*z^2+17496*v*y*z^2-5103*x^2*z^2
            -5103*y^4*z+5832*x^2*y^2*z-6561*y^2*z-729*x^4*z+6561*x^2*z
            -2187*y^4+4374*x^2*y^2-2187*x^4,
       64*z^9-10368*z^7+1296*y^2*z^6-1296*x^2*z^6-34992*y^2*z^5-34992*x^2*z^5
             +419904*z^5+174960*y^2*z^4-174960*x^2*z^4-10935*y^4*z^3
             -56862*x^2*y^2*z^3+314928*y^2*z^3-10935*x^4*z^3+314928*x^2*z^3
             +118098*y^4*z^2-118098*x^4*z^2+59049*y^4*z-118098*x^2*y^2*z
             +59049*x^4*z+19683*y^6-59049*x^2*y^4+59049*x^4*y^2-19683*x^6]" '(u v x y z)))))
    (&body)))
    

(test buchberger-containment
  "Check ideal containment"
  (with-fixture buchberger-criterion-context 
      ()
    (let ((fl-copy (mapcar #'copy-instance fl))
          (gb-copy (mapcar #'copy-instance gb)))
      (loop 
         for i from 0 below (length fl)
         do
           (is (universal-zerop (normal-form (elt fl i) gb))
               "Failed with I=~S~%" I))
      ;; GB should not change in the process
      (is (universal-equalp gb gb-copy))
      ;; FL should not change either
      (is (universal-equalp fl fl-copy)))))

(test buchberger-criterion-with-normal-form
  "Buchberger criterion using normal form"
  (with-fixture buchberger-criterion-context 
      ()
    (let ((fl-copy (mapcar #'copy-instance fl))
          (gb-copy (mapcar #'copy-instance gb)))
      (loop 
         for i from 0 below (length gb)
         do
           (loop 
              for j from (1+ i) below (length gb)
              do
                (is (universal-zerop (normal-form (s-polynomial (elt gb i) (elt gb j)) gb))
                    "Failed with I=~S, J=~S~%" I J)))
      (is (universal-equalp fl fl-copy))
      (is (universal-equalp gb gb-copy)))))

(test buchberger-criterion-with-pseudo-division
  "Buchberger criterion using pseudo-division"
  (with-fixture buchberger-criterion-context ()
    (let ((fl-copy (mapcar #'copy-instance fl))
          (gb-copy (mapcar #'copy-instance gb)))
      (loop 
         for i from 0 below (length gb)
         do
           (loop 
              for j from (1+ i) below (length gb)
              do
                (is (universal-zerop (second (multiple-value-list (poly-pseudo-divide (s-polynomial (elt gb i) (elt gb j)) gb))))
                    "Failed with I=~S, J=~S~%" I J)))
      (is (universal-equalp fl fl-copy))
      (is (universal-equalp gb gb-copy)))))

(test buchberger-criterion-with-grobner-test
  "Buchberger criterion using grobner-test"
  (with-fixture buchberger-criterion-context ()
    (is-true (grobner-test gb fl))
    )
  )


(run! 'division-suite)
(format t "All tests done!~%")