|
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_algorithm_M.h $ 00015 // $Id: modular_gcd_utcf_algorithm_M.h 47339 2008-12-10 10:22:51Z hemmer $ 00016 // 00017 // 00018 // Author(s) : Dominik Huelse <dominik.huelse@gmx.de> 00019 // Michael Hemmer <mhemmer@uni-mainz.de> 00020 // 00021 // ============================================================================ 00022 00028 #ifndef CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_ALGORITHM_M_H 00029 #define CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_ALGORITHM_M_H 1 00030 00031 #include <CGAL/basic.h> 00032 #include <CGAL/Residue.h> 00033 #include <CGAL/Polynomial/modular_gcd.h> 00034 #include <CGAL/Polynomial.h> 00035 #include <CGAL/Polynomial/Cached_extended_euclidean_algorithm.h> 00036 #include <CGAL/Scalar_factor_traits.h> 00037 #include <CGAL/Chinese_remainder_traits.h> 00038 #include <CGAL/Cache.h> 00039 #include <CGAL/Real_timer.h> 00040 00041 // algorithm M for integer polynomials, without denominator bound 00042 00043 namespace CGAL { 00044 00045 namespace CGALi{ 00046 template <class NT> Polynomial<NT> gcd_utcf_UFD(Polynomial<NT>,Polynomial<NT>); 00047 00048 00049 template <class NT> 00050 Polynomial< Polynomial<NT> > modular_gcd_utcf_algorithm_M( 00051 const Polynomial< Polynomial<NT> >& FF1 , 00052 const Polynomial< Polynomial<NT> >& FF2 ){ 00053 return gcd_utcf_UFD(FF1, FF2); 00054 } 00055 00056 template <class NT> 00057 Polynomial<NT> modular_gcd_utcf_algorithm_M( 00058 const Polynomial<NT>& FF1 , 00059 const Polynomial<NT>& FF2 ){ 00060 00061 // Enforce IEEE double precision and to nearest before using modular arithmetic 00062 CGAL::Protect_FPU_rounding<true> pfr(CGAL_FE_TONEAREST); 00063 00064 // std::cout << "start modular_gcd_utcf_algorithm_M " << std::endl; 00065 #ifdef CGAL_MODULAR_GCD_TIMER 00066 timer_init.start(); 00067 #endif 00068 typedef Polynomial<NT> Poly; 00069 typedef Polynomial_traits_d<Poly> PT; 00070 00071 typedef typename PT::Innermost_coefficient_type IC; 00072 00073 // will paly the role of content 00074 typedef typename CGAL::Scalar_factor_traits<Poly>::Scalar Scalar; 00075 00076 typedef typename CGAL::Modular_traits<Poly>::Residue_type MPoly; 00077 typedef typename CGAL::Modular_traits<Scalar>::Residue_type MScalar; 00078 00079 typedef Chinese_remainder_traits<Poly> CRT; 00080 typename CRT::Chinese_remainder chinese_remainder; 00081 00082 CGAL::Real_timer timer; 00083 00084 00085 if(FF1.is_zero()){ 00086 if(FF2.is_zero()){ 00087 return Poly(1);// TODO: return 0 for CGAL 00088 } 00089 else{ 00090 // std::cout<<"\nFF1 is zero"<<std::endl; 00091 00092 return CGAL::canonicalize(FF2); 00093 } 00094 } 00095 if(FF2.is_zero()){ 00096 return CGAL::canonicalize(FF1); 00097 } 00098 if(FF1.degree() == 0 || FF2.degree() == 0){ 00099 Poly result; 00100 result = CGAL::gcd(FF1.content(),FF2.content()); 00101 return CGAL::canonicalize(result); 00102 } 00103 00104 Poly F1 = CGAL::canonicalize(FF1); 00105 Poly F2 = CGAL::canonicalize(FF2); 00106 00107 Scalar f1 = scalar_factor(F1.lcoeff()); // ilcoeff(F1) 00108 Scalar f2 = scalar_factor(F2.lcoeff()); // ilcoeff(F2) 00109 Scalar g_ = scalar_factor(f1,f2); 00110 00111 //std::cout <<" g_ : "<< g_ << std::endl; 00112 00113 bool solved = false; 00114 int prime_index = -1; 00115 int n = 0; // number of lucky primes 00116 int degree_F1 = F1.degree(); 00117 int degree_F2 = F2.degree(); 00118 int degree_e = std::min(degree_F1,degree_F2); 00119 00120 MScalar mg_; 00121 MPoly mF1,mF2,mG_; 00122 00123 typename CRT::Scalar_type p,q,pq,s,t; 00124 Poly Gs,H1s,H2s, Gs_old; // s =^ star 00125 #ifdef CGAL_MODULAR_GCD_TIMER 00126 timer_init.stop(); 00127 #endif 00128 00129 while(!solved){ 00130 do{ 00131 //--------------------------------------- 00132 //choose prime not deviding f1 or f2 00133 MScalar tmp1, tmp2; 00134 do{ 00135 int current_prime = -1; 00136 prime_index++; 00137 if(prime_index >= 2000){ 00138 std::cerr<<"primes in the array exhausted"<<std::endl; 00139 current_prime = CGALi::get_next_lower_prime(current_prime); 00140 } 00141 else{ 00142 current_prime = CGALi::primes[prime_index]; 00143 } 00144 CGAL::Residue::set_current_prime(current_prime); 00145 #ifdef CGAL_MODULAR_GCD_TIMER 00146 timer_image.start(); 00147 #endif 00148 tmp1 = CGAL::modular_image(f1); 00149 tmp2 = CGAL::modular_image(f2); 00150 #ifdef CGAL_MODULAR_GCD_TIMER 00151 timer_image.stop(); 00152 #endif 00153 } 00154 while(!(( tmp1 != 0 ) && ( tmp2 != 0 ))); 00155 // -------------------------------------- 00156 // invoke gcd for current prime 00157 #ifdef CGAL_MODULAR_GCD_TIMER 00158 timer_image.start(); 00159 #endif 00160 mg_ = CGAL::modular_image(g_); 00161 mF1 = CGAL::modular_image(F1); 00162 mF2 = CGAL::modular_image(F2); 00163 #ifdef CGAL_MODULAR_GCD_TIMER 00164 timer_image.stop(); 00165 // replace mG_ = CGAL::gcd (mF1,mF2)*MPoly(mg_); for multivariat 00166 timer_gcd.start(); 00167 #endif 00168 mG_ = CGAL::gcd(mF1,mF2)*MPoly(mg_); 00169 #ifdef CGAL_MODULAR_GCD_TIMER 00170 timer_gcd.stop(); 00171 #endif 00172 00173 //mH1 = CGAL::integral_div(mF1,mG_); 00174 //mH2 = CGAL::integral_div(mF2,mG_); 00175 //--------------------------------------- 00176 // return if G is constant 00177 if (mG_ == MPoly(1)) return Poly(1); 00178 // -------------------------------------- 00179 }// repeat until mG_ degree is less equal the known bound 00180 // check prime 00181 while( mG_.degree() > degree_e); 00182 00183 if( mG_.degree() < degree_e ){ 00184 if( n != 0 ) std::cout << "UNLUCKY PRIME !!"<< std::endl; 00185 00186 // restart chinese remainder 00187 // ignore previous unlucky primes 00188 n=1; 00189 degree_e= mG_.degree(); 00190 }else{ 00191 CGAL_postcondition( mG_.degree() == degree_e); 00192 n++; // increase number of lucky primes 00193 } 00194 00195 // -------------------------------------- 00196 // try chinese remainder 00197 00198 // std::cout <<" chinese remainder round :" << n << std::endl; 00199 typename CGAL::Modular_traits<Poly>::Modular_image_representative inv_map; 00200 if(n == 1){ 00201 // init chinese remainder 00202 q = CGAL::Residue::get_current_prime(); // implicit ! 00203 Gs_old = Gs = inv_map(mG_); 00204 00205 //H1s_old = H1s = inv_map(mH1); 00206 //H2s_old = H2s = inv_map(mH2); 00207 }else{ 00208 // continue chinese remainder 00209 00210 p = CGAL::Residue::get_current_prime(); // implicit! 00211 00212 Gs_old = Gs ; 00213 //H1s_old = H1s ; 00214 //H2s_old = H2s ; 00215 #ifdef CGAL_MODULAR_GCD_TIMER 00216 timer_CR.start(); 00217 #endif 00218 // chinese_remainder(q,Gs ,p,inv_map(mG_),pq,Gs); 00219 // cached_extended_euclidean_algorithm(q,p,s,t); 00220 CGALi::Cached_extended_euclidean_algorithm 00221 <typename CRT::Scalar_type, 1> ceea; 00222 ceea(q,p,s,t); 00223 pq =p*q; 00224 chinese_remainder(q,p,pq,s,t,Gs,inv_map(mG_),Gs); 00225 #ifdef CGAL_MODULAR_GCD_TIMER 00226 timer_CR.stop(); 00227 #endif 00228 q=pq; 00229 } 00230 00231 try{ 00232 if( n != 1 && Gs == Gs_old ){ 00233 Poly r1,r2; 00234 #ifdef CGAL_MODULAR_GCD_TIMER 00235 timer_division.start(); 00236 #endif 00237 00238 typedef CGAL::Algebraic_structure_traits< Poly > ASTE_Poly; 00239 typename ASTE_Poly::Divides divides; 00240 00241 bool div1=divides(Gs,g_*F1,H1s); 00242 bool div2=divides(Gs,g_*F2,H2s); 00243 if (div1 && div2){ 00244 solved = true; 00245 } 00246 // this is the old code 00247 // NT dummy; 00248 // Poly::euclidean_division(g_*F1,Gs,H1s,r1); 00249 // Poly::euclidean_division(g_*F2,Gs,H2s,r2); 00250 // if (r1.is_zero() && r2.is_zero()) 00251 // solved = true; 00252 00253 #ifdef CGAL_MODULAR_GCD_TIMER 00254 timer_division.stop(); 00255 #endif 00256 // std::cout << "number of primes used : "<< n << std::endl; 00257 } // end while 00258 00259 }catch(...){} 00260 00261 } 00262 00263 00264 //TODO CGAL: change this to multivariat content 00265 // Scalar scalar_content_f1 = scalar_factor(FF1); 00266 // Scalar scalar_content_f2 = scalar_factor(FF2); 00267 // Scalar scalar_content_gcd = CGAL::gcd(scalar_content_f1,scalar_content_f2); 00268 // Poly result = CGAL::canonicalize(Gs)*Poly(scalar_content_gcd); 00269 // return result; 00270 00271 return CGAL::canonicalize(Gs); 00272 00273 } 00274 00275 } // namespace CGALi 00276 00277 } // namespace CGAL 00278 00279 #endif // CGAL_POLYNOMIAL_MODULAR_GCD_UTCF_ALGORITHM_M_H
1.7.6.1
