polynomial.h

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

  }

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

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


  NT content() const
  /*{\Mop returns the content of |\Mvar| (the gcd of its coefficients).}*/
  { CGAL_assertion( degree()>=0 );
    return gcd_of_range(this->ptr()->coeff.begin(),this->ptr()->coeff.end());
  }


  /*{\Mtext 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<NT>
   operator - <>  (const Polynomial<NT>&);
                          
  friend Polynomial<NT>
    operator + <> (const Polynomial<NT>&, const Polynomial<NT>&);

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

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

  friend 
  Polynomial<NT>  operator / <> 
  (const Polynomial<NT>& p1, const Polynomial<NT>& p2);
  
  /*{\Mbinopfunc 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 |NT[x]|.  The correct division algorithm is chosen
  according to a number type traits class.
  If |Number_type_traits<NT>::Has_gcd == Tag_true| then the division is
  done by \emph{pseudo division} based on a |gcd| operation of |NT|.  If
  |Number_type_traits<NT>::Has_gcd == Tag_false| then the division is done
  by \emph{euclidean division} based on the division operation of the
  field |NT|.

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

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

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

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

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

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


  Polynomial<NT>& operator += (const Polynomial<NT>& 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<NT>& operator -= (const Polynomial<NT>& 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<NT>& operator *= (const Polynomial<NT>& p1)
  { (*this)=(*this)*p1; return (*this); }

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  void minus_offsetmult(const Polynomial<NT>& p, const NT& b, int k)
  { CGAL_assertion(!this->is_shared());
    Polynomial<NT> 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);
  }

};

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

//############ POLYNOMIAL< INT > ###################################
// CLASS TEMPLATE int: 
/*{\Xsubst 
 iterator_traits<Forward_iterator>::value_type#int
<>#
}*/

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

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

 |int[x]| becomes a unique factorization domain, if the number type
 |int| 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 int NT;
  /*{\Xtypemember the component type representing the coefficients.}*/

  typedef Handle_for< Polynomial_rep<int> > Base;
  typedef Polynomial_rep<int> 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<int>(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<int>() ) {}
   
  Polynomial(const int& a0)
  /*{\Xcreate introduces a variable |\Mvar| of type |\Mname| representing
  the constant polynomial $a_0$.}*/
    : Base(Polynomial_rep<int>(a0)) { reduce(); }


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

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



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

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

  //protected: // accessing coefficients internally:
  int& 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 int& 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.}*/

  int eval_at(const int& r) const
  /*{\Xop evaluates the polynomial at |r|.}*/
  { CGAL_assertion( degree()>=0 );
    int 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 int& leading_coeff = this->ptr()->coeff.back();
    if (leading_coeff < int(0)) return (CGAL::NEGATIVE);
    if (leading_coeff > int(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]==int(0); }

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

  int content() const
  /*{\Xop returns the content of |\Mvar| (the gcd of its coefficients).
  \precond Requires |int| to provide a |gcd| 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<int>
    operator + <> (const Polynomial<int>&, const Polynomial<int>&);

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