more comments / using Gmpz again

Modular_arithmetic/examples/Modular_arithmetic/modular_filter.cpp

@@ -1,36 +1,24 @@
+/* Modular arithmetic can be used as a filter, in this example modular 
+   arithmetic is used to avoid unnecessary gcd computations of polynomials. 
+   A gcd computation can be very costly due to coefficient growth within the 
+   Euclidean algorithm. 
+   
+   The general idea is that firstly the gcd is computed with respect 
+   to one prime only. If this modular gcd is constant we can (in most cases) 
+   conclude that the actual gcd is constant as well. 
+   
+   For this purpose the example introduces the function may_have_common_factor.
+   Note that there are two versions of this function, namely for the case 
+   that the coefficient type is Modularizable and that it is not. 
+   If the type is not Modularizable the filter is just not applied and the 
+   function returns true. 
+*/
 #include <CGAL/basic.h>
+
+#ifdef CGAL_USE_GMP
+
+#include <CGAL/Gmpz.h>        
 #include <CGAL/Polynomial.h>
-#include <CGAL/Modular_traits.h>
-
-//fwd: try to apply modular filtering 
-template< typename Polynomial >
-Polynomial modular_filtered_gcd(const Polynomial& p1, const Polynomial& p2);
-
-// fwd: function if Polynomial is Modularizable
-template< typename Polynomial >
-bool may_have_common_factor(
-    const Polynomial& p1, const Polynomial& p2, CGAL::Tag_true);
-
-// fwd: function if Polynomial is not Modularizable
-template< typename Polynomial >
-bool may_have_common_factor(
-    const Polynomial& p1, const Polynomial& p2, CGAL::Tag_false);
-
-
-
-template< typename Polynomial >
-Polynomial modular_filtered_gcd(const Polynomial& p1, const Polynomial& p2){
-  typedef CGAL::Modular_traits<Polynomial> MT;
-  typedef typename MT::Is_modularizable Is_modularizable;
-    
-  // Try to avoid actual gcd computation 
-  if (may_have_common_factor(p1,p2, Is_modularizable())){
-    // Compute gcd, since the filter indicates a common factor
-    return CGAL::gcd(p1,p2);
-  }else{
-    return Polynomial(1); // return trivial gcd 
-  }
-}
 
 // Function in case  Polynomial is Modularizable
 template< typename Polynomial >
@@ -44,12 +32,18 @@ bool may_have_common_factor(
   typedef typename MT::Modular_image Modular_image;
   MPolynomial mp1 = Modular_image()(p1);
   MPolynomial mp2 = Modular_image()(p2);
+  
+  // check for unlucky primes, the polynomials should not lose a degree 
+  typename CGAL::Polynomial_traits_d<Polynomial>::Degree  degree;
+  typename CGAL::Polynomial_traits_d<MPolynomial>::Degree mdegree;
+  if ( degree(p1) != mdegree(mp1)) return true;   
+  if ( degree(p2) != mdegree(mp2)) return true; 
+
   // compute gcd for modular images 
   MPolynomial mg  = CGAL::gcd(mp1,mp2);
   
-  typename CGAL::Polynomial_traits_d<MPolynomial>::Degree degree;
   // if the modular gcd is not trivial: return true 
-  if ( degree(mg) > 0 ){
+  if ( mdegree(mg) > 0 ){
     std::cout << "The gcd may be non trivial" << std::endl;
     return true;
   }else{
@@ -66,27 +60,55 @@ bool may_have_common_factor(
   return true; 
 }
 
+template< typename Polynomial >
+Polynomial modular_filtered_gcd(const Polynomial& p1, const Polynomial& p2){
+  typedef CGAL::Modular_traits<Polynomial> MT;
+  typedef typename MT::Is_modularizable Is_modularizable;
+    
+  // Try to avoid actual gcd computation 
+  if (may_have_common_factor(p1,p2, Is_modularizable())){
+    // Compute gcd, since the filter indicates a common factor
+    return CGAL::gcd(p1,p2);
+  }else{
+    typename CGAL::Polynomial_traits_d<Polynomial>::Univariate_content  content;
+    return CGAL::gcd(content(p1),content(p2)); // return trivial gcd 
+  }
+}
+
 int main(){
   CGAL::set_pretty_mode(std::cout);
     
-  typedef long NT; 
+  typedef CGAL::Gmpz NT; 
   typedef CGAL::Polynomial<NT> Poly; 
 
-  Poly  f1(NT(2), NT(7), NT(1)); 
-  Poly  f2(NT(3), NT(1), NT(4));
+  Poly  f1(NT(2), NT(6), NT(4)); 
+  Poly  f2(NT(12), NT(4), NT(8));
   Poly  f3(NT(3), NT(4));
     
   std::cout << "f1        : " << f1 << std::endl;
   std::cout << "f2        : " << f2 << std::endl;
-  std::cout << "f3        : " << f3 << std::endl;
 
+  std::cout << "compute modular filtered gcd(f1,f2): " << std::endl;
+  Poly g1 = modular_filtered_gcd(f1,f2);
+  std::cout << "gcd(f1,f2): " << g1 << std::endl;
+
+  std::cout << std::endl;
   Poly p1 = f1*f3;
   Poly p2 = f2*f3;
     
+  std::cout << "f3        : " << f3 << std::endl;
   std::cout << "p1=f1*f3  : " << p1 << std::endl;
   std::cout << "p2=f2*f3  : " << p2 << std::endl;
 
-  std::cout << "modular filtered gcd: " << std::endl;
-  Poly g = modular_filtered_gcd(p1,p2);
-  std::cout << "gcd(p1,p2): " << g << std::endl;
+  std::cout << "compute modular filtered gcd(p1,p2): " << std::endl;
+  Poly g2 = modular_filtered_gcd(p1,p2);
+  std::cout << "gcd(p1,p2): " << g2 << std::endl;
 }
+
+#else
+
+int main (){
+  std::cout << " This examples needs GMP! " << std::endl; 
+}
+
+#endif