// Copyright (c) 2002-2008 Max-Planck-Institute Saarbruecken (Germany) // // This file is part of CGAL (www.cgal.org); you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public License as // published by the Free Software Foundation; either version 3 of the License, // or (at your option) any later version. // // Licensees holding a valid commercial license may use this file in // accordance with the commercial license agreement provided with the software. // // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // // $URL$ // $Id$ // // // Author(s) : Dominik Huelse // Michael Hemmer // // ============================================================================ /*! \file CGAL/Polynomial/modular_gcd_utcf_dfai.h provides gcd for Polynomials, based on Modular arithmetic. */ #ifndef CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_DFAI_H #define CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_DFAI_H 1 #include #include #include #include #include #include #include #include #include #include #include #include namespace CGAL { namespace internal{ template Polynomial gcd_utcf_Integral_domain(Polynomial,Polynomial); template Polynomial< Polynomial > modular_gcd_utcf_dfai( const Polynomial< Polynomial >& FF1 , const Polynomial< Polynomial >& FF2 ){ return internal::gcd_utcf_Integral_domain(FF1, FF2); } // algorithm just computes Gs using the denominator for // algebraic integers and checks it using pseudo division. template Polynomial modular_gcd_utcf_dfai( const Polynomial& FF1_ , const Polynomial& FF2_ ){ // Enforce IEEE double precision and to nearest before using modular arithmetic CGAL::Protect_FPU_rounding pfr(CGAL_FE_TONEAREST); typedef Polynomial Poly; Poly FF1 = FF1_; Poly FF2 = FF2_; typedef Polynomial_traits_d PT; typedef typename PT::Innermost_coefficient_type IC; typename Coercion_traits::Cast ictp; typename PT::Construct_innermost_coefficient_const_iterator_range range; typename PT::Innermost_leading_coefficient ilcoeff; typedef Algebraic_extension_traits ANT; typename ANT::Denominator_for_algebraic_integers dfai; typename ANT::Normalization_factor nfac; typedef typename CGAL::Scalar_factor_traits SFT; typedef typename SFT::Scalar Scalar; typedef typename CGAL::Modular_traits::Residue_type MPoly; typedef typename CGAL::Modular_traits::Residue_type MScalar; typedef Chinese_remainder_traits CRT; typename CRT::Chinese_remainder chinese_remainder; CGAL::Real_timer timer; if(FF1.is_zero()){ if(FF2.is_zero()){ return Poly(1); } else{ return CGAL::canonicalize(FF2); } } if(FF2.is_zero()){; return CGAL::canonicalize(FF1); } if(FF1.degree() == 0 || FF2.degree() == 0){ return Poly(1); } #ifdef CGAL_MODULAR_GCD_TIMER timer_init.start(); #endif Poly F1 = CGAL::canonicalize(FF1); Poly F2 = CGAL::canonicalize(FF2); // This is the most important part of the (dfai) algorithm, it computes the // multiplictive denominator bound according to the algebraic extension // This is needed to guarantee a termination of the Chenese Remainder, // i.e. ensures that Gs can be expressed in terms of algebraic integers // dfai = denominator for algebraic integers IC denom; { Poly tmp = F1+F2; denom = dfai(range(tmp).first, range(tmp).second); } denom *= nfac(denom); Scalar denominator = scalar_factor(denom); //std::cout <<" F1*denom*nafc: " << F1 < Vector; Vector prs_degrees_old, prs_degrees_new; int current_prime = -1; while(!solved){ do{ //--------------------------------------- //choose prime not deviding f1 or f2 MScalar tmp1, tmp2; do{ prime_index++; if(prime_index >= 2000){ std::cerr<<"primes exhausted"< prs_degrees_new); // check that everything went fine if( prs_degrees_old < prs_degrees_new ){ if( n != 0 ) std::cout << "UNLUCKY PRIME !!"<< std::endl; // restart chinese remainder // ignore previous unlucky primes n=1; prs_degrees_old = prs_degrees_new; }else{ CGAL_postcondition( prs_degrees_old == prs_degrees_new); n++; // increase number of lucky primes } // -------------------------------------- // try chinese remainder typename CGAL::Modular_traits::Modular_image_representative inv_map; if(n == 1){ // init chinese remainder q = CGAL::Residue::get_current_prime(); // implicit ! Gs_old = Gs = inv_map(mG_); }else{ // continue chinese remainder p = CGAL::Residue::get_current_prime(); // implicit! Gs_old = Gs ; #ifdef CGAL_MODULAR_GCD_TIMER timer_CR.start(); #endif internal::Cached_extended_euclidean_algorithm ceea; ceea(q,p,s,t); pq =p*q; chinese_remainder(q,p,pq,s,t,Gs,inv_map(mG_),Gs); #ifdef CGAL_MODULAR_GCD_TIMER timer_CR.stop(); #endif q=pq; } try{ // to catch error in case the extension is not correct yet. if( n != 1 && Gs == Gs_old ){ Poly r1,r2; NT dummy; #ifdef CGAL_MODULAR_GCD_TIMER timer_division.start(); #endif typedef CGAL::Algebraic_structure_traits< Poly > ASTE_Poly; typename ASTE_Poly::Divides divides; FF1*=ictp(ilcoeff(Gs)*denom); FF2*=ictp(ilcoeff(Gs)*denom); bool div1=divides(Gs,FF1,H1s); bool div2=divides(Gs,FF2,H2s); //This is the old code: // Poly::pseudo_division(F1,Gs,H1s,r1,dummy); // Poly::pseudo_division(F2,Gs,H2s,r2,dummy); if (div1 && div2){ solved = true; // std::cout << "number of primes used : "<< n << std::endl; } // if (r1.is_zero() && r2.is_zero()){ // solved = true; // std::cout << "number of primes used : "<< n << std::endl; // } #ifdef CGAL_MODULAR_GCD_TIMER timer_division.stop(); #endif } } catch(...){} } return CGAL::canonicalize(Gs); } } // namespace internal } // namespace CGAL #endif // CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_DFAI_H