From b6aa6d8644f99d7b68b9f4ddd3103d5562d6efc0 Mon Sep 17 00:00:00 2001
From: Michael Hemmer <hemmer@mpi-inf.mpg.de>
Date: Fri, 26 Sep 2008 14:52:49 +0000
Subject: [PATCH] apply modular filter to gcd functors

---
 .../Polynomial/Algebraic_structure_traits.h   | 60 ++++++++++++-------
 Polynomial/include/CGAL/Polynomial_traits_d.h | 10 +++-
 2 files changed, 45 insertions(+), 25 deletions(-)

diff --git a/Polynomial/include/CGAL/Polynomial/Algebraic_structure_traits.h b/Polynomial/include/CGAL/Polynomial/Algebraic_structure_traits.h
index bdfafcb3de0..eadcbd40531 100644
--- a/Polynomial/include/CGAL/Polynomial/Algebraic_structure_traits.h
+++ b/Polynomial/include/CGAL/Polynomial/Algebraic_structure_traits.h
@@ -190,29 +190,43 @@ class Polynomial_algebraic_structure_traits_base< POLY, Unique_factorization_dom
   public:
     typedef Unique_factorization_domain_tag Algebraic_category;
     
-    class Gcd 
-      : public std::binary_function< POLY, POLY, POLY > {
-      public:
-        POLY operator()( const POLY& x, const POLY& y ) const {
-            typedef Algebraic_structure_traits<POLY> AST;
-            typename AST::Integral_division idiv;
-            typename AST::Unit_part upart; 
-
-          // First: the extreme cases and negative sign corrections.          
-          if (x == POLY(0)) {
-              if (y == POLY(0))  
-                  return POLY(0);
-              return idiv(y,upart(y));
-          }
-          if (y == POLY(0))
-              return idiv(x,upart(x));
-          
-          return CGALi::gcd(x,y);
-        }
-        
-        CGAL_IMPLICIT_INTEROPERABLE_BINARY_OPERATOR( POLY )
-        
-    };
+  class Gcd 
+    : public std::binary_function< POLY, POLY, POLY > {
+    typedef typename Polynomial_traits_d<POLY>::Multivariate_content Mcontent;
+    typedef typename Mcontent::result_type ICoeff; 
+    
+    ICoeff gcd_help(const ICoeff& x, const ICoeff& y, Field_tag) const {
+      return ICoeff(1);
+    }
+    ICoeff gcd_help(const ICoeff& x, const ICoeff& y, 
+        Unique_factorization_domain_tag) const {
+      return CGAL::gcd(x,y);
+    }
+  public:
+    POLY operator()( const POLY& x, const POLY& y ) const {
+      typedef Algebraic_structure_traits<POLY> AST;
+      typename AST::Integral_division idiv;
+      typename AST::Unit_part upart; 
+      
+      // First: the extreme cases and negative sign corrections.          
+      if (CGAL::is_zero(x)) {
+        if (CGAL::is_zero(y))
+          return POLY(0);
+        return idiv(y,upart(y));
+      }
+      if (CGAL::is_zero(y))
+        return idiv(x,upart(x));
+      
+      if (CGALi::may_have_common_factor(x,y))
+        return CGALi::gcd(x,y);
+      else{        
+        typename Algebraic_structure_traits<ICoeff>::Algebraic_category category;  
+        return POLY(gcd_help(Mcontent()(x),Mcontent()(y), category));
+      }
+    }
+    
+    CGAL_IMPLICIT_INTEROPERABLE_BINARY_OPERATOR( POLY )  
+  };
 };
 
 // Clone this for a EuclideanRing
diff --git a/Polynomial/include/CGAL/Polynomial_traits_d.h b/Polynomial/include/CGAL/Polynomial_traits_d.h
index 2e04f2e1ecd..51718c0417d 100644
--- a/Polynomial/include/CGAL/Polynomial_traits_d.h
+++ b/Polynomial/include/CGAL/Polynomial_traits_d.h
@@ -1111,9 +1111,15 @@ public:
     :public std::binary_function<Polynomial_d, Polynomial_d, Polynomial_d> {
     Polynomial_d
     operator()(const Polynomial_d& p, const Polynomial_d& q) const {
-      if (CGAL::is_zero(p) && CGAL::is_zero(q)) 
+      if (CGAL::is_zero(p) && CGAL::is_zero(q)){
         return Polynomial_d(0);
-      return CGALi::gcd_utcf(p,q);
+      }
+      // apply modular filter first
+      if (CGALi::may_have_common_factor(p,q)){
+        return CGALi::gcd_utcf(p,q);
+      }else{
+        return Polynomial_d(1);
+      }
     }
   };