diff --git a/Modular_arithmetic/examples/Modular_arithmetic/modular_filter.cpp b/Modular_arithmetic/examples/Modular_arithmetic/modular_filter.cpp index 68ad2196986..d852d50c6b2 100644 --- a/Modular_arithmetic/examples/Modular_arithmetic/modular_filter.cpp +++ b/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 + +#ifdef CGAL_USE_GMP + +#include #include -#include - -//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 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::Degree degree; + typename CGAL::Polynomial_traits_d::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::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 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::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 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