polynomial.h

很多二维 三维几何计算算法 C++ 类库

// Copyright (c) 2000  Max-Planck-Institute Saarbruecken (Germany).
// All rights reserved.
//
// 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; version 2.1 of the License.
// See the file LICENSE.LGPL distributed with CGAL.
//
// 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: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Nef_2/include/CGAL/Nef_2/Polynomial.h $
// $Id: Polynomial.h 37034 2007-03-12 17:34:47Z hemmer $
// 
//
// Author(s)     : Michael Seel <[e-mail]>
//                 Andreas Fabri <[e-mail]>

#ifndef CGAL_POLYNOMIAL_H
#define CGAL_POLYNOMIAL_H

#include <CGAL/basic.h>
#include <CGAL/kernel_assertions.h>
#include <CGAL/Handle_for.h>
#include <CGAL/number_type_basic.h>
#include <CGAL/number_utils.h>
#include <CGAL/IO/io.h>
#include <cstddef>
#undef CGAL_NEF_DEBUG
#define CGAL_NEF_DEBUG 3
#include <CGAL/Nef_2/debug.h>
#include <vector>



CGAL_BEGIN_NAMESPACE

namespace Nef {

template <class NT> class Polynomial_rep;

// SPECIALIZE_CLASS(NT,int double) START
// CLASS TEMPLATE NT: 
template <class NT> class Polynomial;
// CLASS TEMPLATE int: 
template <> class Polynomial<int> ;
// CLASS TEMPLATE double: 
template <> class Polynomial<double> ;
// SPECIALIZE_CLASS(NT,int double) END



/*{\Mtext \headerline{Range template}}*/


template <class Forward_iterator>
typename std::iterator_traits<Forward_iterator>::value_type 
gcd_of_range(Forward_iterator its, Forward_iterator ite, Unique_factorization_domain_tag)
/*{\Mfunc calculates the greates common divisor of the
set of numbers $\{ |*its|, |*++its|, \ldots, |*it| \}$ of type |NT|,
where |++it == ite| and |NT| is the value type of |Forward_iterator|. 
\precond there exists a pairwise gcd operation |NT gcd(NT,NT)| and 
|its!=ite|.}*/
{ CGAL_assertion(its!=ite);
  typedef typename std::iterator_traits<Forward_iterator>::value_type NT;
  NT res = *its++;
  for(; its!=ite; ++its) res = 
    (*its==0 ? res : CGAL_NTS gcd(res, *its));
  if (res==0) res = 1;
  return res;
}

template <class Forward_iterator>
typename std::iterator_traits<Forward_iterator>::value_type 
gcd_of_range(Forward_iterator its, Forward_iterator ite, Integral_domain_without_division_tag)
/*{\Mfunc calculates the greates common divisor of the
set of numbers $\{ |*its|, |*++its|, \ldots, |*it| \}$ of type |NT|,
where |++it == ite| and |NT| is the value type of |Forward_iterator|. 
\precond |its!=ite|.}*/
{ CGAL_assertion(its!=ite);
  typedef typename std::iterator_traits<Forward_iterator>::value_type NT;
  NT res = *its++;
  for(; its!=ite; ++its) res = 
    (*its==0 ? res : 1);
  if (res==0) res = 1;
  return res;
}

template <class Forward_iterator>
typename std::iterator_traits<Forward_iterator>::value_type 
gcd_of_range(Forward_iterator its, Forward_iterator ite)
{
    
  typedef typename std::iterator_traits<Forward_iterator>::value_type NT;
  typedef CGAL::Algebraic_structure_traits<NT> AST;
  return gcd_of_range(its,ite,typename AST::Algebraic_category());
}


template <class NT>  inline Polynomial<NT>
  operator - (const Polynomial<NT>&);
template <class NT>  Polynomial<NT>
  operator + (const Polynomial<NT>&, const Polynomial<NT>&);
template <class NT>  Polynomial<NT>
  operator - (const Polynomial<NT>&, const Polynomial<NT>&);
template <class NT>  Polynomial<NT>
  operator * (const Polynomial<NT>&, const Polynomial<NT>&);
template <class NT> inline Polynomial<NT>
  operator / (const Polynomial<NT>&, const Polynomial<NT>&);
template <class NT> inline Polynomial<NT>
  operator % (const Polynomial<NT>&, const Polynomial<NT>&);

template<class NT> CGAL::Sign 
  sign(const Polynomial<NT>& p);

template <class NT> double 
  to_double(const Polynomial<NT>& p) ;
template <class NT> bool 
  is_valid(const Polynomial<NT>& p) ;
template <class NT> bool 
  is_finite(const Polynomial<NT>& p) ;

template<class NT>  
  std::ostream& operator << (std::ostream& os, const Polynomial<NT>& p);
template <class NT>  
  std::istream& operator >> (std::istream& is, Polynomial<NT>& p);


  // lefthand side
template<class NT> inline Polynomial<NT> operator + 
  (const NT& num, const Polynomial<NT>& p2);
template<class NT> inline Polynomial<NT> operator - 
  (const NT& num, const Polynomial<NT>& p2);
template<class NT> inline Polynomial<NT> operator * 
  (const NT& num, const Polynomial<NT>& p2);
template<class NT> inline Polynomial<NT> operator / 
  (const NT& num, const Polynomial<NT>& p2);
template<class NT> inline Polynomial<NT> operator % 
  (const NT& num, const Polynomial<NT>& p2);

  // righthand side
template<class NT> inline Polynomial<NT> operator + 
  (const Polynomial<NT>& p1, const NT& num);
template<class NT> inline Polynomial<NT> operator - 
  (const Polynomial<NT>& p1, const NT& num);
template<class NT> inline Polynomial<NT> operator * 
  (const Polynomial<NT>& p1, const NT& num);
template<class NT> inline Polynomial<NT> operator / 
  (const Polynomial<NT>& p1, const NT& num);
template<class NT> inline Polynomial<NT> operator % 
  (const Polynomial<NT>& p1, const NT& num);


  // lefthand side
template<class NT> inline bool operator ==  
  (const NT& num, const Polynomial<NT>& p);
template<class NT> inline bool operator != 
  (const NT& num, const Polynomial<NT>& p);
template<class NT> inline bool operator <  
  (const NT& num, const Polynomial<NT>& p);
template<class NT> inline bool operator <=  
  (const NT& num, const Polynomial<NT>& p);
template<class NT> inline bool operator >  
  (const NT& num, const Polynomial<NT>& p);
template<class NT> inline bool operator >=  
  (const NT& num, const Polynomial<NT>& p);

  // righthand side
template<class NT> inline bool operator ==
  (const Polynomial<NT>& p, const NT& num);
template<class NT> inline bool operator !=
  (const Polynomial<NT>& p, const NT& num);
template<class NT> inline bool operator < 
  (const Polynomial<NT>& p, const NT& num);
template<class NT> inline bool operator <= 
  (const Polynomial<NT>& p, const NT& num);
template<class NT> inline bool operator > 
  (const Polynomial<NT>& p, const NT& num);
template<class NT> inline bool operator >=
  (const Polynomial<NT>& p, const NT& num);



template <class pNT> class Polynomial_rep { 
  typedef pNT NT;

  typedef std::vector<NT> Vector;


  typedef typename Vector::size_type      size_type;
  typedef typename Vector::iterator       iterator;
  typedef typename Vector::const_iterator const_iterator;
  Vector coeff;

  Polynomial_rep() : coeff(0) {}
  Polynomial_rep(const NT& n) : coeff(1) { coeff[0]=n; }
  Polynomial_rep(const NT& n, const NT& m) : coeff(2)
    { coeff[0]=n; coeff[1]=m; }
  Polynomial_rep(const NT& a, const NT& b, const NT& c) : coeff(3)
    { coeff[0]=a; coeff[1]=b; coeff[2]=c; }
  Polynomial_rep(size_type s) : coeff(s,NT(0)) {}


  template <class Forward_iterator>
  Polynomial_rep(std::pair<Forward_iterator,Forward_iterator> poly) 
    : coeff(poly.first,poly.second) {}


  void reduce() 
  { while ( coeff.size()>1 && coeff.back()==NT(0) ) coeff.pop_back(); }

  friend class Polynomial<pNT>;
  friend class Polynomial<int>;
  friend class Polynomial<double>;
  friend std::istream& operator >> <> (std::istream&, Polynomial<NT>&);

};

// SPECIALIZE_CLASS(NT,int double) START
// CLASS TEMPLATE NT: 
/*{\Msubst 
typename iterator_traits<Forward_iterator>::value_type#NT
<>#
}*/

/*{\Manpage{Polynomial}{NT}{Polynomials in one variable}{p}}*/

template <class pNT> class Polynomial : 
  public Handle_for< Polynomial_rep<pNT> >
{
/*  {\Mdefinition An instance |\Mvar| of the data type |\Mname| represents
  a polynomial $p = a_0 + a_1 x + \ldots a_d x^d $ from the ring |NT[x]|. 
  The data type offers standard ring operations and a sign operation which 
  determines the sign for the limit process $x \rightarrow \infty$. 

  |NT[x]| becomes a unique factorization domain, if the number type
  |NT| is either a field type (1) or a unique factorization domain
  (2). In both cases there's a polynomial division operation defined.}
*/

  /*{\Mtypes 5}*/
  public:
  typedef pNT NT;
  
  typedef Handle_for< Polynomial_rep<NT> > Base;
  typedef Polynomial_rep<NT> Rep;
  typedef typename Rep::Vector    Vector;
  typedef typename Rep::size_type size_type;
  typedef typename Rep::iterator  iterator;

  typedef typename Rep::const_iterator const_iterator;
  /*{\Mtypemember a random access iterator for read-only access to the
  coefficient vector.}*/

  //protected:
  void reduce() { this->ptr()->reduce(); }
  Vector& coeffs() { return this->ptr()->coeff; }
  const Vector& coeffs() const { return this->ptr()->coeff; }
  Polynomial(size_type s) : Base( Polynomial_rep<NT>(s) ) {}
  // creates a polynomial of degree s-1


  /*{\Mcreation 3}*/
  public:

  Polynomial()
  /*{\Mcreate introduces a variable |\Mvar| of type |\Mname| of undefined
  value. }*/
    : Base( Polynomial_rep<NT>() ) {}
   
  Polynomial(const NT& a0)
  /*{\Mcreate introduces a variable |\Mvar| of type |\Mname| representing
  the constant polynomial $a_0$.}*/
    : Base(Polynomial_rep<NT>(a0)) { reduce(); }

  Polynomial(NT a0, NT a1) 
  /*{\Mcreate introduces a variable |\Mvar| of type |\Mname| representing
  the polynomial $a_0 + a_1 x$. }*/
    : Base(Polynomial_rep<NT>(a0,a1)) { reduce(); }

  Polynomial(const NT& a0, const NT& a1,const NT& a2)
  /*{\Mcreate introduces a variable |\Mvar| of type |\Mname| representing
  the polynomial $a_0 + a_1 x + a_2 x^2$. }*/
    : Base(Polynomial_rep<NT>(a0,a1,a2)) { reduce(); }

  template <class Forward_iterator>
  Polynomial(std::pair<Forward_iterator, Forward_iterator> poly) 
  /*{\Mcreate introduces a variable |\Mvar| of type |\Mname| representing
  the polynomial whose coefficients are determined by the iterator range,
  i.e. let $(a_0 = |*first|, a_1 = |*++first|, \ldots a_d = |*it|)$, 
  where |++it == last| then |\Mvar| stores the polynomial $a_1 + a_2 x + 
  \ldots a_d x^d$.}*/
    : Base(Polynomial_rep<NT>(poly)) { reduce(); }


  // KILL double START
  Polynomial(double n) : Base(Polynomial_rep<NT>(NT(n))) { reduce(); }
  Polynomial(double n1, double n2) 
    : Base(Polynomial_rep<NT>(NT(n1),NT(n2))) { reduce(); }
  // KILL double END

  // KILL int START
  Polynomial(int n) : Base(Polynomial_rep<NT>(NT(n))) { reduce(); }
  Polynomial(int n1, int n2)
    : Base(Polynomial_rep<NT>(NT(n1),NT(n2))) { reduce(); }
  // KILL int END

  Polynomial(const Polynomial<NT>& p) : Base(p) {}

  //protected: // accessing coefficients internally:
  NT& coeff(unsigned int i) 
  { CGAL_assertion(!this->is_shared() && i<(this->ptr()->coeff.size()));
    return this->ptr()->coeff[i]; 
  }
  public:

  /*{\Moperations 3 3 }*/
  const_iterator begin() const { return this->ptr()->coeff.begin(); }
  /*{\Mop a random access iterator pointing to $a_0$.}*/
  const_iterator end()   const { return this->ptr()->coeff.end(); }
  /*{\Mop a random access iterator pointing beyond $a_d$.}*/

  int degree() const 
  { return this->ptr()->coeff.size()-1; } 
  /*{\Mop the degree of the polynomial.}*/

  const NT& operator[](unsigned int i) const 
  { CGAL_assertion( i<(this->ptr()->coeff.size()) );
    return this->ptr()->coeff[i]; }
  //{\Marrop the coefficient $a_i$ of the polynomial.}

  NT operator()(unsigned int i) const
  { 
    if(i<(this->ptr()->coeff.size()))
      return this->ptr()->coeff[i];
    return 0;
  }

  NT eval_at(const NT& r) const
  /*{\Mop evaluates the polynomial at |r|.}*/
  { CGAL_assertion( degree()>=0 );
    NT res = this->ptr()->coeff[0], x = r;
    for(int i=1; i<=degree(); ++i) 
    { res += this->ptr()->coeff[i]*x; x*=r; }
    return res; 
  }

  CGAL::Sign sign() const
  /*{\Mop returns the sign of the limit process for $x \rightarrow \infty$
  (the sign of the leading coefficient).}*/
  { const NT& leading_coeff = this->ptr()->coeff.back();
    if (leading_coeff < NT(0)) return (CGAL::NEGATIVE);
    if (leading_coeff > NT(0)) return (CGAL::POSITIVE);
    return CGAL::ZERO;