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

;;; -*-  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)

;; Unless NGROBNER system loaded by ASDF,
;; load the dependencies directly
#-ngrobner
(progn
  (require :utils "utils")
  (require :copy "copy")
  (require :ring "ring")
  (require :integer-ring "integer-ring")
  (require :monom "monom")
  (require :polynomial "polynomial")
  (require :infix "infix")
  (require :symbolic-polynomial "symbolic-polynomial")
  (require :division "division"))

(defpackage #:5am-division
  (:use :cl :it.bese.fiveam :monom :polynomial :infix :symbolic-polynomial :division :integer-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 'integer-ring :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))
    )
  )

(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)))))


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