|
BWAPI
|
00001 // Copyright (c) 2002-2008 Max-Planck-Institute Saarbruecken (Germany) 00002 // 00003 // This file is part of CGAL (www.cgal.org); you can redistribute it and/or 00004 // modify it under the terms of the GNU Lesser General Public License as 00005 // published by the Free Software Foundation; version 2.1 of the License. 00006 // See the file LICENSE.LGPL distributed with CGAL. 00007 // 00008 // Licensees holding a valid commercial license may use this file in 00009 // accordance with the commercial license agreement provided with the software. 00010 // 00011 // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE 00012 // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. 00013 // 00014 // $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.5-branch/Polynomial/include/CGAL/Polynomial/modular_gcd_utcf_dfai.h $ 00015 // $Id: modular_gcd_utcf_dfai.h 47339 2008-12-10 10:22:51Z hemmer $ 00016 // 00017 // 00018 // Author(s) : Dominik Huelse <dominik.huelse@gmx.de> 00019 // Michael Hemmer <hemmer@mpi-inf.mpg.de> 00020 // 00021 // ============================================================================ 00022 00026 #ifndef CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_DFAI_H 00027 #define CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_DFAI_H 1 00028 00029 #include <CGAL/basic.h> 00030 #include <CGAL/Scalar_factor_traits.h> 00031 #include <CGAL/Chinese_remainder_traits.h> 00032 #include <CGAL/Algebraic_extension_traits.h> 00033 #include <CGAL/primes.h> 00034 #include <CGAL/Residue.h> 00035 #include <CGAL/Polynomial/modular_gcd.h> 00036 #include <CGAL/Polynomial/Cached_extended_euclidean_algorithm.h> 00037 #include <CGAL/Polynomial.h> 00038 #include <CGAL/Polynomial/modular_gcd_utils.h> 00039 #include <CGAL/Cache.h> 00040 #include <CGAL/Real_timer.h> 00041 00042 00043 00044 namespace CGAL { 00045 00046 namespace CGALi{ 00047 00048 template <class NT> Polynomial<NT> 00049 gcd_utcf_Integral_domain(Polynomial<NT>,Polynomial<NT>); 00050 00051 00052 template <class NT> 00053 Polynomial< Polynomial<NT> > modular_gcd_utcf_dfai( 00054 const Polynomial< Polynomial<NT> >& FF1 , 00055 const Polynomial< Polynomial<NT> >& FF2 ){ 00056 return CGALi::gcd_utcf_Integral_domain(FF1, FF2); 00057 } 00058 00059 // algorithm just computes Gs using the denominator for 00060 // algebraic integers and checks it using pseudo division. 00061 template <class NT> 00062 Polynomial<NT> modular_gcd_utcf_dfai( 00063 const Polynomial<NT>& FF1_ , 00064 const Polynomial<NT>& FF2_ ){ 00065 00066 // Enforce IEEE double precision and to nearest before using modular arithmetic 00067 CGAL::Protect_FPU_rounding<true> pfr(CGAL_FE_TONEAREST); 00068 00069 typedef Polynomial<NT> Poly; 00070 Poly FF1 = FF1_; 00071 Poly FF2 = FF2_; 00072 00073 typedef Polynomial_traits_d<Poly> PT; 00074 typedef typename PT::Innermost_coefficient_type IC; 00075 00076 typename Coercion_traits<Poly,IC>::Cast ictp; 00077 typename PT::Construct_innermost_coefficient_const_iterator_range range; 00078 typename PT::Innermost_leading_coefficient ilcoeff; 00079 00080 typedef Algebraic_extension_traits<IC> ANT; 00081 typename ANT::Denominator_for_algebraic_integers dfai; 00082 typename ANT::Normalization_factor nfac; 00083 00084 typedef typename CGAL::Scalar_factor_traits<Poly> SFT; 00085 typedef typename SFT::Scalar Scalar; 00086 00087 typedef typename CGAL::Modular_traits<Poly>::Residue_type MPoly; 00088 typedef typename CGAL::Modular_traits<Scalar>::Residue_type MScalar; 00089 00090 typedef Chinese_remainder_traits<Poly> CRT; 00091 typename CRT::Chinese_remainder chinese_remainder; 00092 00093 CGAL::Real_timer timer; 00094 00095 if(FF1.is_zero()){ 00096 if(FF2.is_zero()){ 00097 return Poly(1); 00098 } 00099 else{ 00100 return CGAL::canonicalize(FF2); 00101 } 00102 } 00103 if(FF2.is_zero()){; 00104 return CGAL::canonicalize(FF1); 00105 } 00106 00107 if(FF1.degree() == 0 || FF2.degree() == 0){ 00108 return Poly(1); 00109 } 00110 00111 #ifdef CGAL_MODULAR_GCD_TIMER 00112 timer_init.start(); 00113 #endif 00114 00115 Poly F1 = CGAL::canonicalize(FF1); 00116 Poly F2 = CGAL::canonicalize(FF2); 00117 00118 00119 // This is the most important part of the (dfai) algorithm, it computes the 00120 // multiplictive denominator bound according to the algebraic extension 00121 // This is needed to guarantee a termination of the Chenese Remainder, 00122 // i.e. ensures that Gs can be expressed in terms of algebraic integers 00123 // dfai = denominator for algebraic integers 00124 00125 IC denom; 00126 { 00127 Poly tmp = F1+F2; 00128 denom = dfai(range(tmp).first, range(tmp).second); 00129 } 00130 denom *= nfac(denom); 00131 00132 Scalar denominator = scalar_factor(denom); 00133 //std::cout <<" F1*denom*nafc: " << F1 <<std::endl; 00134 //std::cout <<" F2*denom*nfac: " << F2 <<std::endl; 00135 00136 Scalar f1 = scalar_factor(F1.lcoeff()); // ilcoeff(F1) 00137 Scalar f2 = scalar_factor(F2.lcoeff()); // ilcoeff(F2) 00138 Scalar g_ = scalar_factor(f1,f2) * denominator; 00139 00140 bool solved = false; 00141 int prime_index = -1; 00142 00143 int n = 0; // number of lucky primes 00144 00145 MScalar mg_; 00146 MPoly mF1,mF2,mG_,mQ,mR; 00147 00148 typename CRT::Scalar_type p,q,pq,s,t; 00149 Poly Gs,H1s,H2s,Gs_old; // s =^ star 00150 00151 #ifdef CGAL_MODULAR_GCD_TIMER 00152 timer_init.stop(); 00153 #endif 00154 00155 typedef std::vector<int> Vector; 00156 Vector prs_degrees_old, prs_degrees_new; 00157 int current_prime = -1; 00158 00159 00160 while(!solved){ 00161 do{ 00162 //--------------------------------------- 00163 //choose prime not deviding f1 or f2 00164 MScalar tmp1, tmp2; 00165 do{ 00166 prime_index++; 00167 if(prime_index >= 2000){ 00168 std::cerr<<"primes exhausted"<<std::endl; 00169 // return CGAL::CGALi::gcd_utcf_Integral_domain(FF1,FF2); 00170 current_prime = CGALi::get_next_lower_prime(current_prime); 00171 } 00172 else{ 00173 current_prime = CGALi::primes[prime_index]; 00174 } 00175 CGAL_assertion(current_prime != -1); 00176 CGAL::Residue::set_current_prime(current_prime); 00177 #ifdef CGAL_MODULAR_GCD_TIMER 00178 timer_image.start(); 00179 #endif 00180 tmp1 = modular_image(f1); 00181 tmp2 = modular_image(f2); 00182 #ifdef CGAL_MODULAR_GCD_TIMER 00183 timer_image.stop(); 00184 #endif 00185 } 00186 while(!(( tmp1 != 0 ) 00187 && ( tmp2 != 0 ) 00188 && (denominator%current_prime) != 0 )); 00189 00190 // -------------------------------------- 00191 // invoke gcd for current prime 00192 #ifdef CGAL_MODULAR_GCD_TIMER 00193 timer_image.start(); 00194 #endif 00195 mg_ = CGAL::modular_image(g_); 00196 mF1 = CGAL::modular_image(FF1); 00197 mF2 = CGAL::modular_image(FF2); 00198 #ifdef CGAL_MODULAR_GCD_TIMER 00199 timer_image.stop(); 00200 #endif 00201 00202 #ifdef CGAL_MODULAR_GCD_TIMER 00203 timer_gcd.start(); 00204 #endif 00205 00206 // compute gcd over Field[x] 00207 00208 prs_degrees_new.clear(); 00209 CGAL_precondition(mF1 != MPoly(0)); 00210 CGAL_precondition(mF2 != MPoly(0)); 00211 while (!mF2.is_zero()) { 00212 // MPoly::euclidean_division(mF1, mF2, mQ, mR); 00213 euclidean_division_obstinate(mF1, mF2, mQ, mR); 00214 mF1 = mF2; mF2 = mR; 00215 prs_degrees_new.push_back(mR.degree()); 00216 } 00217 mF1 /= mF1.lcoeff(); 00218 mG_ = mF1; 00219 00220 if(n==0) 00221 prs_degrees_old = prs_degrees_new; 00222 00223 mG_ = mG_ * MPoly(mg_); 00224 00225 #ifdef CGAL_MODULAR_GCD_TIMER 00226 timer_gcd.stop(); 00227 #endif 00228 00229 // return if G is constant 00230 if (mG_ == MPoly(1)) return Poly(1); 00231 // use ordinary algorithm if prs sequence is too short 00232 // if (prs_degrees_new.size() <= 2) 00233 // return CGALi::gcd_utcf_Integral_domain(F1, F2); 00234 // -------------------------------------- 00235 } 00236 // repeat until mG_ degree is less equal the known bound 00237 // this is now tested by the prs degree sequence 00238 while(prs_degrees_old > prs_degrees_new); 00239 // check that everything went fine 00240 if( prs_degrees_old < prs_degrees_new ){ 00241 if( n != 0 ) std::cout << "UNLUCKY PRIME !!"<< std::endl; 00242 00243 // restart chinese remainder 00244 // ignore previous unlucky primes 00245 n=1; 00246 prs_degrees_old = prs_degrees_new; 00247 }else{ 00248 CGAL_postcondition( prs_degrees_old == prs_degrees_new); 00249 n++; // increase number of lucky primes 00250 } 00251 00252 // -------------------------------------- 00253 // try chinese remainder 00254 typename CGAL::Modular_traits<Poly>::Modular_image_representative inv_map; 00255 if(n == 1){ 00256 // init chinese remainder 00257 q = CGAL::Residue::get_current_prime(); // implicit ! 00258 Gs_old = Gs = inv_map(mG_); 00259 }else{ 00260 // continue chinese remainder 00261 p = CGAL::Residue::get_current_prime(); // implicit! 00262 Gs_old = Gs ; 00263 #ifdef CGAL_MODULAR_GCD_TIMER 00264 timer_CR.start(); 00265 #endif 00266 CGALi::Cached_extended_euclidean_algorithm <Scalar, 2> ceea; 00267 ceea(q,p,s,t); 00268 pq =p*q; 00269 chinese_remainder(q,p,pq,s,t,Gs,inv_map(mG_),Gs); 00270 #ifdef CGAL_MODULAR_GCD_TIMER 00271 timer_CR.stop(); 00272 #endif 00273 q=pq; 00274 } 00275 00276 try{ 00277 // to catch error in case the extension is not correct yet. 00278 if( n != 1 && Gs == Gs_old ){ 00279 Poly r1,r2; NT dummy; 00280 #ifdef CGAL_MODULAR_GCD_TIMER 00281 timer_division.start(); 00282 #endif 00283 typedef CGAL::Algebraic_structure_traits< Poly > 00284 ASTE_Poly; 00285 typename ASTE_Poly::Divides divides; 00286 00287 00288 FF1*=ictp(ilcoeff(Gs)*denom); 00289 FF2*=ictp(ilcoeff(Gs)*denom); 00290 bool div1=divides(Gs,FF1,H1s); 00291 bool div2=divides(Gs,FF2,H2s); 00292 00293 00294 //This is the old code: 00295 // Poly::pseudo_division(F1,Gs,H1s,r1,dummy); 00296 // Poly::pseudo_division(F2,Gs,H2s,r2,dummy); 00297 if (div1 && div2){ 00298 solved = true; 00299 // std::cout << "number of primes used : "<< n << std::endl; 00300 } 00301 00302 // if (r1.is_zero() && r2.is_zero()){ 00303 // solved = true; 00304 // std::cout << "number of primes used : "<< n << std::endl; 00305 // } 00306 00307 #ifdef CGAL_MODULAR_GCD_TIMER 00308 timer_division.stop(); 00309 #endif 00310 } 00311 } 00312 catch(...){} 00313 } 00314 00315 return CGAL::canonicalize(Gs); 00316 00317 } 00318 00319 } // namespace CGALi 00320 } // namespace CGAL 00321 00322 #endif // CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_DFAI_H
1.7.6.1
