BWAPI: SPAR/AIModule/BWTA/vendors/CGAL/CGAL/Polynomial/sturm_habicht

Go to the documentation of this file.
00001 // Copyright (c) 2008 Max-Planck-Institute Saarbruecken (Germany).
00002 // All rights reserved.
00003 //
00004 // This file is part of CGAL (www.cgal.org); you can redistribute it and/or
00005 // modify it under the terms of the GNU Lesser General Public License as
00006 // published by the Free Software Foundation; version 2.1 of the License.
00007 // See the file LICENSE.LGPL distributed with CGAL.
00008 //
00009 // Licensees holding a valid commercial license may use this file in
00010 // accordance with the commercial license agreement provided with the software.
00011 //
00012 // This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
00013 // WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00014 //
00015 // $URL: svn+ssh://afabri@scm.gforge.inria.fr/svn/cgal/trunk/Polynomial/include/CGAL/Polynomial.h $
00016 // $Id: Polynomial.h 46502 2008-10-28 08:36:59Z hemmer $
00017 //
00018 //
00019 // Author(s)     : Michael Kerber <mkerber@mpi-inf.mpg.de>
00020 //
00021 // ============================================================================
00022 
00023 #ifndef CGAL_POLYNOMIAL_STURM_HABICHT
00024 #define CGAL_POLYNOMIAL_STURM_HABICHT 1
00025 
00026 #include <vector>
00027 #include <algorithm>
00028 #include <CGAL/Polynomial/bezout_matrix.h>
00029 #include <CGAL/Polynomial/subresultants.h>
00030 
00031 CGAL_BEGIN_NAMESPACE
00032 
00033 namespace CGALi {
00034 
00048   template<typename Polynomial_traits_d,typename OutputIterator>
00049     OutputIterator prs_principal_sturm_habicht_sequence
00050   (typename Polynomial_traits_d::Polynomial_d P,
00051      OutputIterator out) {
00052       
00053     typedef typename Polynomial_traits_d::Coefficient_type NT;
00054     typename Polynomial_traits_d::Get_coefficient coeff;
00055     typename Polynomial_traits_d::Degree degree;
00056 
00057     std::vector<typename Polynomial_traits_d::Polynomial_d> stha;
00058     CGAL::CGALi::sturm_habicht_sequence<Polynomial_traits_d>
00059         (P,std::back_inserter(stha));
00060     for(int i=0; i<static_cast<int>(stha.size()); i++) {
00061       int d = degree(stha[i]);
00062       CGAL_assertion(d<=i);
00063       if(d<i) {
00064         *out++ = NT(0);
00065       } else {
00066         *out++ = coeff(stha[i],i);
00067       }
00068     }
00069     return out;
00070 
00071   } 
00072 
00086   template<typename Polynomial_traits_d,typename OutputIterator>
00087     OutputIterator bezout_principal_sturm_habicht_sequence
00088       (typename Polynomial_traits_d::Polynomial_d P,
00089        OutputIterator out) {
00090 
00091     typedef typename Polynomial_traits_d::Polynomial_d Polynomial;
00092     typedef typename Polynomial_traits_d::Coefficient_type NT;
00093     typename Polynomial_traits_d::Leading_coefficient lcoeff;
00094     typename Polynomial_traits_d::Differentiate diff;
00095 
00096     Polynomial Px = diff(P);
00097 
00098     CGAL::CGALi::Simple_matrix<NT> M 
00099         = CGAL::CGALi::polynomial_subresultant_matrix<Polynomial_traits_d>
00100             (P,Px,1);
00101     int n = static_cast<int>(M.row_dimension());
00102     for(int i=0; i<n; i++) {
00103       if((n-1-i)%4==0 || (n-1-i)%4==1) {
00104         *out++ = -M[n-1-i][n-1-i];
00105       } else {
00106         *out++ = M[n-1-i][n-1-i];
00107       }
00108     }
00109     *out++=lcoeff(Px);
00110     *out++=lcoeff(P);
00111     
00112     return out;
00113 
00114   } 
00115 
00116 
00120   template<typename Polynomial_traits_d,
00121            typename OutputIterator1,
00122            typename OutputIterator2> 
00123     void prs_first_two_sturm_habicht_coefficients
00124       (typename Polynomial_traits_d::Polynomial_d P,
00125        OutputIterator1 pstha,
00126        OutputIterator2 copstha) {
00127     
00128     typedef typename Polynomial_traits_d::Polynomial_d Polynomial;
00129     typedef typename Polynomial_traits_d::Coefficient_type NT;
00130     typename Polynomial_traits_d::Get_coefficient coeff;
00131     typename Polynomial_traits_d::Degree degree;
00132 
00133     std::vector<Polynomial> stha;
00134     int n = degree(P);
00135 
00136     sturm_habicht_sequence<Polynomial_traits_d>(P,std::back_inserter(stha));
00137     CGAL_assertion(static_cast<int>(stha.size())==n+1);
00138     for(int i=0;i<=n;i++) {
00139       int d = degree(stha[i]);
00140       CGAL_assertion(d<=i);
00141       if(d<i) {
00142         *pstha++ = NT(0);
00143       } else {
00144         *pstha++ = coeff(stha[i],i);
00145       }
00146     }
00147     for(int i=1;i<=n;i++) {
00148       int d = degree(stha[i]);
00149       CGAL_assertion(d<=i);
00150       if(d<i-1) {
00151         *copstha++ = NT(0);
00152       } else {
00153         *copstha++ = coeff(stha[i],i-1);
00154       }
00155     }
00156     return;
00157   }
00158 
00161   template<typename Polynomial_traits_d,
00162            typename OutputIterator1,
00163            typename OutputIterator2> 
00164     void bezout_first_two_sturm_habicht_coefficients
00165       (typename Polynomial_traits_d::Polynomial_d P,
00166        OutputIterator1 pstha,
00167        OutputIterator2 copstha) {
00168 
00169     typedef typename Polynomial_traits_d::Polynomial_d Polynomial;
00170     typedef typename Polynomial_traits_d::Coefficient_type NT;
00171     typename Polynomial_traits_d::Get_coefficient coeff;
00172     typename Polynomial_traits_d::Degree degree;
00173     typename Polynomial_traits_d::Leading_coefficient lcoeff;
00174     typename Polynomial_traits_d::Differentiate diff;
00175 
00176     Polynomial Px=diff(P);
00177     CGAL::CGALi::Simple_matrix<NT> M 
00178         = CGAL::CGALi::polynomial_subresultant_matrix<Polynomial_traits_d>
00179             (P,Px,2);
00180     int n = static_cast<int>(M.row_dimension());
00181     for(int i=0; i<n; i++) {
00182       if((n-1-i)%4==0 || (n-1-i)%4==1) {
00183         *pstha++ = -M[n-1-i][n-1-i];
00184       } else {
00185         *pstha++ = M[n-1-i][n-1-i];
00186       }
00187     }
00188     *pstha++ = lcoeff(Px);
00189     *pstha++ = lcoeff(P);
00190     for(int i=1; i<n; i++) {
00191       if(n-i-1%4==0 || n-i-1%4==1) {
00192         *copstha++ = -M[n-i-1][n-i];
00193       } else {
00194         *copstha++ = M[n-1-i][n-i];
00195       }
00196     }
00197     *copstha++ = coeff(Px,degree(Px)-1);
00198     *copstha++ = coeff(P,degree(P)-1);
00199   }
00200 
00201   
00202     // the general function for CGAL::Integral_domain_without_division_tag
00203   template < typename Polynomial_traits_d,
00204              typename OutputIterator1, 
00205              typename OutputIterator2 > inline
00206       void 
00207       first_two_sturm_habicht_coefficients_
00208         (typename Polynomial_traits_d::Polynomial_d A, 
00209          OutputIterator1 pstha,
00210          OutputIterator2 copstha,
00211          CGAL::Integral_domain_without_division_tag){
00212 
00213         bezout_first_two_sturm_habicht_coefficients<Polynomial_traits_d>
00214             (A,pstha,copstha);
00215   
00216     }
00217 
00218     // the general function for CGAL::Integral_domain_tag
00219     template < typename Polynomial_traits_d,
00220                typename OutputIterator1, 
00221                typename OutputIterator2 > inline
00222       void 
00223       first_two_sturm_habicht_coefficients_
00224         (typename Polynomial_traits_d::Polynomial_d A, 
00225          OutputIterator1 pstha,
00226          OutputIterator2 copstha,
00227          CGAL::Integral_domain_tag) {
00228 
00229         return prs_first_two_sturm_habicht_coefficients<Polynomial_traits_d>
00230             (A,pstha,copstha);
00231   
00232     }
00233     
00234     template < typename Polynomial_traits_d,
00235                typename OutputIterator1, 
00236                typename OutputIterator2 > inline
00237       void 
00238       first_two_sturm_habicht_coefficients_
00239         (typename Polynomial_traits_d::Polynomial_d A,
00240          OutputIterator1 pstha,
00241          OutputIterator2 copstha) {
00242 
00243         typedef typename Polynomial_traits_d::Coefficient_type NT;
00244 
00245         typedef typename 
00246             CGAL::Algebraic_structure_traits<NT>::Algebraic_category 
00247             Algebraic_category;
00248         first_two_sturm_habicht_coefficients_<Polynomial_traits_d>
00249             (A,pstha,copstha,Algebraic_category());
00250     }
00251 
00252     // the general function for CGAL::Integral_domain_without_division_tag
00253     template <typename Polynomial_traits_d,typename OutputIterator> inline 
00254       OutputIterator 
00255       principal_sturm_habicht_sequence_
00256         (typename Polynomial_traits_d::Polynomial_d A, 
00257          OutputIterator out,
00258          CGAL::Integral_domain_without_division_tag) {
00259       
00260         return bezout_principal_sturm_habicht_sequence<Polynomial_traits_d>
00261             (A,out);
00262     }
00263       
00264     // the specialization for CGAL::Integral_domain_tag
00265     template <typename Polynomial_traits_d,typename OutputIterator> inline
00266       OutputIterator
00267       principal_sturm_habicht_sequence_
00268         (typename Polynomial_traits_d::Polynomial_d A, 
00269          OutputIterator out,
00270          CGAL::Integral_domain_tag) {
00271     
00272       return prs_principal_sturm_habicht_sequence<Polynomial_traits_d>(A,out);
00273     }
00274 
00275     template <typename Polynomial_traits_d,typename OutputIterator> inline
00276       OutputIterator principal_sturm_habicht_sequence_
00277         (typename Polynomial_traits_d::Polynomial_d A,
00278          OutputIterator out) {
00279         
00280       typedef typename Polynomial_traits_d::Coefficient_type NT;
00281         
00282       typedef typename 
00283           CGAL::Algebraic_structure_traits<NT>::Algebraic_category 
00284           Algebraic_category;
00285       return principal_sturm_habicht_sequence_<Polynomial_traits_d>
00286           (A,out,Algebraic_category());  
00287     }
00288 
00289 
00296   template < typename Polynomial_traits_d,
00297              typename OutputIterator1, 
00298              typename OutputIterator2 > inline
00299     void first_two_sturm_habicht_coefficients
00300       (typename Polynomial_traits_d::Polynomial_d A, 
00301        OutputIterator1 pstha,
00302        OutputIterator2 copstha){
00303       
00304       return CGAL::CGALi::first_two_sturm_habicht_coefficients_
00305           <Polynomial_traits_d> (A,pstha,copstha);
00306   }
00307 
00308 
00309     
00310 
00315   template <typename Polynomial_traits_d,typename OutputIterator> inline
00316     OutputIterator
00317     principal_sturm_habicht_sequence
00318       (typename Polynomial_traits_d::Polynomial_d A, 
00319        OutputIterator out){
00320       
00321       return CGAL::CGALi::principal_sturm_habicht_sequence_
00322           <Polynomial_traits_d>(A,out);
00323   }
00324   
00325 
00326 
00330   template<typename Polynomial_traits_d,typename OutputIterator> OutputIterator
00331     sturm_habicht_sequence(typename Polynomial_traits_d::Polynomial_d P, 
00332                            OutputIterator out) {
00333     typedef typename Polynomial_traits_d::Polynomial_d Polynomial;
00334     typename Polynomial_traits_d::Degree degree;
00335     typename Polynomial_traits_d::Differentiate diff;
00336 
00337     int p = degree(P);
00338 
00339     Polynomial P_x = diff(P);
00340 
00341     std::vector<Polynomial> stha;
00342     
00343     CGAL::CGALi::polynomial_subresultants<Polynomial_traits_d>
00344         (P,P_x,std::back_inserter(stha));
00345     stha.push_back(P);
00346 
00347     CGAL_assertion(static_cast<int>(stha.size())==p+1);
00348 
00349     for(int i=0;i<=p; i++) {
00350       if((p-i)%4==0 || (p-i)%4==1) {
00351         *out++ = stha[i];
00352       } else {
00353         *out++ = -stha[i];
00354       }
00355     }
00356     
00357     return out;
00358   }
00359 
00363 template<typename Polynomial_traits_d,
00364     typename OutputIterator1,
00365     typename OutputIterator2,
00366     typename OutputIterator3> 
00367     OutputIterator1
00368     sturm_habicht_sequence_with_cofactors
00369       (typename Polynomial_traits_d::Polynomial_d P,
00370        OutputIterator1 stha_out,
00371        OutputIterator2 cof_out,
00372        OutputIterator3 cofx_out) {
00373 
00374     typedef typename Polynomial_traits_d::Polynomial_d Polynomial;
00375     typename Polynomial_traits_d::Degree degree;
00376     typename Polynomial_traits_d::Differentiate diff;
00377     typename Polynomial_traits_d::Construct_polynomial construct;
00378 
00379     int p = degree(P);
00380 
00381     Polynomial P_x = diff(P);
00382 
00383     std::vector<Polynomial> stha,co_f,co_fx;
00384     
00385     CGAL::CGALi::polynomial_subresultants_with_cofactors<Polynomial_traits_d>
00386         (P,P_x,
00387          std::back_inserter(stha),
00388          std::back_inserter(co_f),
00389          std::back_inserter(co_fx));
00390 
00391     stha.push_back(P);
00392     co_f.push_back(construct(1));
00393     co_fx.push_back(construct(0));
00394 
00395     CGAL_assertion(static_cast<int>(stha.size())==p+1);
00396 
00397     for(int i=0;i<=p; i++) {
00398       if((p-i)%4==0 || (p-i)%4==1) {
00399         *stha_out++ = stha[i];
00400         *cof_out++ = co_f[i];
00401         *cofx_out++ = co_fx[i];
00402       } else {
00403         *stha_out++ = -stha[i];
00404         *cof_out++ = -co_f[i];
00405         *cofx_out++ = -co_fx[i];
00406       }
00407     }
00408     
00409     return stha_out;
00410   }
00411 
00412 } // namespace CGALi
00413 
00414 
00415 
00416       
00417 
00422   template<typename InputIterator>
00423     int number_of_real_roots(InputIterator start,InputIterator end) {
00424     if(start==end) {
00425       return 0;
00426     }
00427     int m = 0;
00428 
00429     CGAL::Sign last_non_zero=CGAL::ZERO; //marks the starting point
00430     CGAL::Sign curr_sign;
00431     int k;
00432     InputIterator el=start;
00433     
00434     //std::cout << "Sign of." << (*el) << std::endl;
00435 
00436     curr_sign=CGAL::sign(*el);
00437 
00438     while(curr_sign==CGAL::ZERO && el!=end) {
00439       el++;
00440       curr_sign=CGAL::sign(*el);
00441     }
00442     if(el==end) return 0;
00443 
00444     last_non_zero=curr_sign;
00445     k=0;
00446     el++;
00447 
00448     while(el!=end) {
00449       curr_sign=CGAL::sign(*el);
00450 
00451       el++;
00452 
00453       if(curr_sign==CGAL::ZERO) {
00454         k++;
00455       }
00456       else {
00457         if(k%2==0) { // k is even
00458           k=k/2;
00459           int pm_one = (curr_sign==last_non_zero ? 1 : -1);
00460           pm_one = (k%2==1) ? -pm_one : pm_one;
00461           m+=pm_one;
00462         }
00463         k=0;
00464         last_non_zero=curr_sign;
00465       }
00466           
00467     }
00468     return m;
00469   }
00470 
00474   template<typename Polynomial_d>
00475     int number_of_real_roots(Polynomial_d f) {
00476       
00477       typedef CGAL::Polynomial_traits_d<Polynomial_d> Poly_traits_d;
00478       typedef typename Poly_traits_d::Coefficient_type Coeff;
00479       std::vector<Coeff> stha;
00480       typename Poly_traits_d::Principal_sturm_habicht_sequence()
00481           (f,std::back_inserter(stha));
00482       return CGAL::number_of_real_roots(stha.begin(),stha.end());
00483   }
00484 
00485 CGAL_END_NAMESPACE
00486 
00487 #endif // CGAL_POLYNOMIAL_STURM_HABICHT