polynomial.h

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

  friend Polynomial<int>
    operator * <> (const Polynomial<int>&, const Polynomial<int>&);

  friend 
  Polynomial<int>  operator / <> 
  (const Polynomial<int>& p1, const Polynomial<int>& p2);
  /*{\Xbinopfunc implements polynomial division of |p1| and |p2|. if
  |p1 = p2 * p3| then |p2| is returned. The result is undefined if |p3|
  does not exist in |int[x]|.  The correct division algorithm is chosen
  according to a number type traits class.
  If |Number_type_traits<int>::Has_gcd == Tag_true| then the division is
  done by \emph{pseudo division} based on a |gcd| operation of |int|.  If
  |Number_type_traits<int>::Has_gcd == Tag_false| then the division is done
  by \emph{euclidean division} based on the division operation of the
  field |int|.

  \textbf{Note} that |int=int| quickly leads to overflow
  errors when using this operation.}*/

  /*{\Xtext \headerline{Static member functions}}*/

  static Polynomial<int> gcd
    (const Polynomial<int>& p1, const Polynomial<int>& p2);
  /*{\Xstatic returns the greatest common divisor of |p1| and |p2|.
  \textbf{Note} that |int=int| quickly leads to overflow errors when
  using this operation.  \precond Requires |int| to be a unique
  factorization domain, i.e. to provide a |gcd| operation.}*/

  static void pseudo_div
    (const Polynomial<int>& f, const Polynomial<int>& g, 
     Polynomial<int>& q, Polynomial<int>& r, int& D);
  /*{\Xstatic implements division with remainder on polynomials of 
  the ring |int[x]|: $D*f = g*q + r$.  \precond |int| is a unique
  factorization domain, i.e., there exists a |gcd| operation and an
  integral division operation on |int|.}*/

  static void euclidean_div 
    (const Polynomial<int>& f, const Polynomial<int>& g, 
     Polynomial<int>& q, Polynomial<int>& r);
  /*{\Xstatic implements division with remainder on polynomials of 
  the ring |int[x]|: $f = g*q + r$.  \precond |int| is a field, i.e.,
  there exists a division operation on |int|.  }*/

  friend double to_double <> (const Polynomial<int>& p);


  Polynomial<int>& operator += (const Polynomial<int>& p1)
  { this->copy_on_write();
    int d = (std::min)(degree(),p1.degree()), i;
    for(i=0; i<=d; ++i) coeff(i) += p1[i];
    while (i<=p1.degree()) this->ptr()->coeff.push_back(p1[i++]);
    reduce(); return (*this); }

  Polynomial<int>& operator -= (const Polynomial<int>& p1)
  { this->copy_on_write();
    int d = (std::min)(degree(),p1.degree()), i;
    for(i=0; i<=d; ++i) coeff(i) -= p1[i];
    while (i<=p1.degree()) this->ptr()->coeff.push_back(-p1[i++]);
    reduce(); return (*this); }

  Polynomial<int>& operator *= (const Polynomial<int>& p1)
  { (*this)=(*this)*p1; return (*this); }

  Polynomial<int>& operator /= (const Polynomial<int>& p1)
  { (*this)=(*this)/p1; return (*this); }

  Polynomial<int>& operator %= (const Polynomial<int>& p1)
  { (*this)=(*this)%p1; return (*this); }


  //------------------------------------------------------------------
  // SPECIALIZE_MEMBERS(int double) START
    
  Polynomial<int>& operator += (const int& num)
  { this->copy_on_write();
    coeff(0) += (int)num; return *this; }

  Polynomial<int>& operator -= (const int& num)
  { this->copy_on_write();
    coeff(0) -= (int)num; return *this; }

  Polynomial<int>& operator *= (const int& num)
  { this->copy_on_write();
    for(int i=0; i<=degree(); ++i) coeff(i) *= (int)num; 
    reduce(); return *this; }

  Polynomial<int>& operator /= (const int& num)
  { this->copy_on_write(); CGAL_assertion(num!=0);
    for(int i=0; i<=degree(); ++i) coeff(i) /= (int)num; 
    reduce(); return *this; }

  Polynomial<int>& operator %= (const int& num)
  { this->copy_on_write(); CGAL_assertion(num!=0);
    for(int i=0; i<=degree(); ++i) coeff(i) %= (int)num; 
    reduce(); return *this; }

// SPECIALIZING_MEMBERS FOR const double& 
    
  Polynomial<int>& operator += (const double& num)
  { this->copy_on_write();
    coeff(0) += (int)num; return *this; }

  Polynomial<int>& operator -= (const double& num)
  { this->copy_on_write();
    coeff(0) -= (int)num; return *this; }

  Polynomial<int>& operator *= (const double& num)
  { this->copy_on_write();
    for(int i=0; i<=degree(); ++i) coeff(i) *= (int)num; 
    reduce(); return *this; }

  Polynomial<int>& operator /= (const double& num)
  { this->copy_on_write(); CGAL_assertion(num!=0);
    for(int i=0; i<=degree(); ++i) coeff(i) /= (int)num; 
    reduce(); return *this; }

  Polynomial<int>& operator %= (const double& num)
  { this->copy_on_write(); CGAL_assertion(num!=0);
    for(int i=0; i<=degree(); ++i) coeff(i) %= (int)num; 
    reduce(); return *this; }

  // SPECIALIZE_MEMBERS(int double) END
  //------------------------------------------------------------------

  void minus_offsetmult(const Polynomial<int>& p, const int& b, int k)
  { CGAL_assertion(!is_shared());
    Polynomial<int> s(size_type(p.degree()+k+1)); // zero entries
    for (int i=k; i <= s.degree(); ++i) s.coeff(i) = b*p[i-k];
    operator-=(s);
  }

};

/*{\Ximplementation This data type is implemented as an item type
via a smart pointer scheme. The coefficients are stored in a vector of
|int| entries.  The simple arithmetic operations $+,-$ take time
$O(d*T(int))$, multiplication is quadratic in the maximal degree of the
arguments times $T(int)$, where $T(int)$ is the time for a corresponding
operation on two instances of the ring type.}*/


//############ POLYNOMIAL< INT > ###################################

//############ POLYNOMIAL< DOUBLE > ################################


// CLASS TEMPLATE double: 
/*{\Xsubst 
 iterator_traits<Forward_iterator>::value_type#double
<>#
}*/

/*{\Xanpage{Polynomial}{double}{Polynomials in one variable}{p}}*/

template <>
class Polynomial<double> : 
  public Handle_for< Polynomial_rep<double> >
{
/*{\Xdefinition 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 |double[x]|. 
The data type offers standard ring operations and a sign operation which 
determines the sign for the limit process $x \rightarrow \infty$. 

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

  /*{\Xtypes 5}*/
  public:
  typedef double NT;
  /*{\Xtypemember the component type representing the coefficients.}*/

  typedef Handle_for< Polynomial_rep<double> > Base;
  typedef Polynomial_rep<double> Rep;
  typedef  Rep::Vector    Vector;
  typedef  Rep::size_type size_type;
  typedef  Rep::iterator  iterator;

  typedef  Rep::const_iterator const_iterator;
  /*{\Xtypemember 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<double>(s) ) {}
  // creates a polynomial of degree s-1


  /*{\Xcreation 3}*/
  public:

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

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

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

  template <class Forward_iterator>
  Polynomial(std::pair<Forward_iterator, Forward_iterator> poly) 
  /*{\Xcreate 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<double>(poly)) { reduce(); }



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

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

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

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

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

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

  double eval_at(const double& r) const
  /*{\Xop evaluates the polynomial at |r|.}*/
  { CGAL_assertion( degree()>=0 );
    double 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
  /*{\Xop returns the sign of the limit process for $x \rightarrow \infty$
  (the sign of the leading coefficient).}*/
  { const double& leading_coeff = this->ptr()->coeff.back();
    if (leading_coeff < double(0)) return (CGAL::NEGATIVE);
    if (leading_coeff > double(0)) return (CGAL::POSITIVE);
    return CGAL::ZERO;
  }

  bool is_zero() const
  /*{\Xop returns true iff |\Mvar| is the zero polynomial.}*/
  { return degree()==0 && this->ptr()->coeff[0]==double(0); }

  Polynomial<double> abs() const
  /*{\Xop returns |-\Mvar| if |\Mvar.sign()==NEGATIVE| and |\Mvar| 
  otherwise.}*/
  { if ( sign()==CGAL::NEGATIVE ) return -*this; return *this; }


  double content() const
  /*{\Xop returns the content of |\Mvar| (the gcd of its coefficients).
  \precond Requires |double| to provide a |gdc| operation.}*/
  { CGAL_assertion( degree()>=0 );
    return gcd_of_range(this->ptr()->coeff.begin(),this->ptr()->coeff.end());
  }
 

  /*{\Xtext Additionally |\Mname| offers standard arithmetic ring
  opertions like |+,-,*,+=,-=,*=|. By means of the sign operation we can
  also offer comparison predicates as $<,>,\leq,\geq$. Where $p_1 < p_2$
  holds iff $|sign|(p_1 - p_2) < 0$. This data type is fully compliant
  to the requirements of CGAL number types. \setopdims{3cm}{2cm}}*/
                     
  friend Polynomial<double>
    operator + <> (const Polynomial<double>&, const Polynomial<double>&);

  friend Polynomial<double>
    operator - <> (const Polynomial<double>&, const Polynomial<double>&);

  friend Polynomial<double>
    operator * <> (const Polynomial<double>&, const Polynomial<double>&);

  friend 
  Polynomial<double>  operator / <> 
  (const Polynomial<double>& p1, const Polynomial<double>& p2);
  /*{\Xbinopfunc implements polynomial division of |p1| and |p2|. if
  |p1 = p2 * p3| then |p2| is returned. The result is undefined if |p3|
  does not exist in |double[x]|.  The correct division algorithm is chosen
  according to the number type traits class.
  If |Number_type_traits<double>::Has_gcd == Tag_true| then the division is
  done by \emph{pseudo division} based on a |gcd| operation of |double|.  If
  |Number_type_traits<double>::Has_gcd == Tag_false| then the division is done
  by \emph{euclidean division} based on the division operation of the
  field |double|.

  \textbf{Note} that |double=int| quickly leads to overflow
  errors when using this operation.}*/

  /*{\Xtext \headerline{Static member functions}}*/

  static Polynomial<double> gcd
    (const Polynomial<double>& p1, const Polynomial<double>& p2);
  /*{\Xstatic returns the greatest common divisor of |p1| and |p2|.
  \textbf{Note} that |double=int| quickly leads to overflow errors when
  using this operation.  \precond Requires |double| to be a unique
  factorization domain, i.e. to provide a |gdc| operation.}*/

  static void pseudo_div
    (const Polynomial<double>& f, const Polynomial<double>& g, 
     Polynomial<double>& q, Polynomial<double>& r, double& D);
  /*{\Xstatic implements division with remainder on polynomials of 
  the ring |double[x]|: $D*f = g*q + r$.  \precond |double| is a unique
  factorization domain, i.e., there exists a |gcd| operation and an
  integral division operation on |double|.}*/

  static void euclidean_div 
    (const Polynomial<double>& f, const Polynomial<double>& g, 
     Polynomial<double>& q, Polynomial<double>& r);
  /*{\Xstatic implements division with remainder on polynomials of 
  the ring |double[x]|: $f = g*q + r$.  \precond |double| is a field, i.e.,
  there exists a division operation on |double|.  }*/

  friend double to_double <> (const Polynomial<double>& p);


  Polynomial<double>& operator += (const Polynomial<double>& p1)
  { this->copy_on_write();
    int d = (std::min)(degree(),p1.degree()), i;
    for(i=0; i<=d; ++i) coeff(i) += p1[i];
    while (i<=p1.degree()) this->ptr()->coeff.push_back(p1[i++]);
    reduce(); return (*this); }

  Polynomial<double>& operator -= (const Polynomial<double>& p1)