polynomial.h

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

  { 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)
  { (*this)=(*this)*p1; return (*this); }

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

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


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

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

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

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

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

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

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

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

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

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

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

  void minus_offsetmult(const Polynomial<double>& p, const double& b, int k)
  { CGAL_assertion(!this->is_shared());
    Polynomial<double> 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
|double| entries.  The simple arithmetic operations $+,-$ take time
$O(d*T(double))$, multiplication is quadratic in the maximal degree of the
arguments times $T(double)$, where $T(double)$ is the time for a corresponding
operation on two instances of the ring type.}*/


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



// SPECIALIZE_CLASS(NT,int double) END

template <class NT> double to_double 
  (const Polynomial<NT>&) 
  { return 0; }

template <class NT> bool is_valid 
  (const Polynomial<NT>& p) 
  { return (CGAL::is_valid(p[0])); }

template <class NT> bool is_finite 
  (const Polynomial<NT>& p) 
  { return CGAL::is_finite(p[0]); }

template <class NT> 
inline
Polynomial<NT> operator - (const Polynomial<NT>& p)
{
  CGAL_assertion(p.degree()>=0);
  Polynomial<NT> res(std::make_pair(p.coeffs().begin(),p.coeffs().end())); 
  typename Polynomial<NT>::iterator it, ite=res.coeffs().end();
  for(it=res.coeffs().begin(); it!=ite; ++it) *it = -*it;
  return res;
}



template <class NT> 
Polynomial<NT> operator + (const Polynomial<NT>& p1, 
                            const Polynomial<NT>& p2)
{ 
  typedef typename Polynomial<NT>::size_type size_type;
  CGAL_assertion(p1.degree()>=0 && p2.degree()>=0);
  bool p1d_smaller_p2d = p1.degree() < p2.degree();
  int min,max,i;
  if (p1d_smaller_p2d) { min = p1.degree(); max = p2.degree(); }
  else                 { max = p1.degree(); min = p2.degree(); }
  Polynomial<NT>  p( size_type(max + 1));
  for (i = 0; i <= min; ++i ) p.coeff(i) = p1[i]+p2[i];
  if (p1d_smaller_p2d)  for (; i <= max; ++i ) p.coeff(i)=p2[i];
  else /* p1d >= p2d */ for (; i <= max; ++i ) p.coeff(i)=p1[i];
  p.reduce();
  return p;
}

template <class NT> 
Polynomial<NT> operator - (const Polynomial<NT>& p1, 
                            const Polynomial<NT>& p2)
{ 
  typedef typename Polynomial<NT>::size_type size_type;
  CGAL_assertion(p1.degree()>=0 && p2.degree()>=0);
  bool p1d_smaller_p2d = p1.degree() < p2.degree();
  int min,max,i;
  if (p1d_smaller_p2d) { min = p1.degree(); max = p2.degree(); }
  else                 { max = p1.degree(); min = p2.degree(); }
  Polynomial<NT>  p( size_type(max+1) );
  for (i = 0; i <= min; ++i ) p.coeff(i)=p1[i]-p2[i];
  if (p1d_smaller_p2d)  for (; i <= max; ++i ) p.coeff(i)= -p2[i];
  else /* p1d >= p2d */ for (; i <= max; ++i ) p.coeff(i)=  p1[i];
  p.reduce();
  return p;
}

template <class NT> 
Polynomial<NT> operator * (const Polynomial<NT>& p1, 
                            const Polynomial<NT>& p2)
{
  typedef typename Polynomial<NT>::size_type size_type;
  CGAL_assertion(p1.degree()>=0 && p2.degree()>=0);
  Polynomial<NT>  p( size_type(p1.degree()+p2.degree()+1) ); 
  // initialized with zeros
  for (int i=0; i <= p1.degree(); ++i)
    for (int j=0; j <= p2.degree(); ++j)
      p.coeff(i+j) += (p1[i]*p2[j]); 
  p.reduce();
  return p;
}

template <class NT> inline
Polynomial<NT> operator / (const Polynomial<NT>& p1, 
                           const Polynomial<NT>& p2)
{ 
    typedef Algebraic_structure_traits<NT> AST;
    return divop(p1,p2, typename AST::Algebraic_category());
}


template <class NT> inline
Polynomial<NT> operator % (const Polynomial<NT>& p1,
			   const Polynomial<NT>& p2)
{ 
    typedef Algebraic_structure_traits<NT> AST;
    return modop(p1,p2, typename AST::Algebraic_category());
}


template <class NT> 
Polynomial<NT> divop (const Polynomial<NT>& p1, 
                      const Polynomial<NT>& p2,
                      Integral_domain_without_division_tag)
{ CGAL_assertion(!p2.is_zero());
  if (p1.is_zero()) {
	return 0;
  }
  Polynomial<NT> q;
  Polynomial<NT> r;
  Polynomial<NT>::euclidean_div(p1,p2,q,r);
  CGAL_postcondition( (p2*q+r==p1) );
  return q;
}


template <class NT> 
Polynomial<NT> divop (const Polynomial<NT>& p1, 
                      const Polynomial<NT>& p2,
                      Unique_factorization_domain_tag)
{ CGAL_assertion(!p2.is_zero());
  if (p1.is_zero()) return 0;
  Polynomial<NT> q,r; NT D; 
  Polynomial<NT>::pseudo_div(p1,p2,q,r,D); 
  CGAL_postcondition( (p2*q+r==p1*Polynomial<NT>(D)) );
  return q/=D;
}


template <class NT> 
Polynomial<NT> modop (const Polynomial<NT>& p1, 
                      const Polynomial<NT>& p2,
                      Integral_domain_without_division_tag)
{ CGAL_assertion(!p2.is_zero());
  if (p1.is_zero()) return 0;
  Polynomial<NT> q,r;
  Polynomial<NT>::euclidean_div(p1,p2,q,r);
  CGAL_postcondition( (p2*q+r==p1) );
  return r;
}


template <class NT> 
Polynomial<NT> modop (const Polynomial<NT>& p1, 
                      const Polynomial<NT>& p2,
                      Unique_factorization_domain_tag)
{ CGAL_assertion(!p2.is_zero());
  if (p1.is_zero()) return 0;
  Polynomial<NT> q,r; NT D; 
  Polynomial<NT>::pseudo_div(p1,p2,q,r,D); 
  CGAL_postcondition( (p2*q+r==p1*Polynomial<NT>(D)) );
  return q/=D;
}


template <class NT> 
inline Polynomial<NT> 
gcd(const Polynomial<NT>& p1, const Polynomial<NT>& p2)
{ return Polynomial<NT>::gcd(p1,p2); }


template <class NT> bool operator == 
  (const Polynomial<NT>& p1, const Polynomial<NT>& p2)
  { return ( (p1-p2).sign() == CGAL::ZERO ); }    

template <class NT> bool operator != 
  (const Polynomial<NT>& p1, const Polynomial<NT>& p2)
  { return ( (p1-p2).sign() != CGAL::ZERO ); }    

template <class NT> bool operator <  
  (const Polynomial<NT>& p1, const Polynomial<NT>& p2)
  { return ( (p1-p2).sign() == CGAL::NEGATIVE ); }    

template <class NT> bool operator <= 
  (const Polynomial<NT>& p1, const Polynomial<NT>& p2)
  { return ( (p1-p2).sign() != CGAL::POSITIVE ); }    

template <class NT> bool operator >  
  (const Polynomial<NT>& p1, const Polynomial<NT>& p2)
  { return ( (p1-p2).sign() == CGAL::POSITIVE ); }    

template <class NT> bool operator >= 
  (const Polynomial<NT>& p1, const Polynomial<NT>& p2)
  { return ( (p1-p2).sign() != CGAL::NEGATIVE ); }    

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

//------------------------------------------------------------------
// SPECIALIZE_FUNCTION(NT,int double) START
// SPECIALIZING inline to :

  // lefthand side
  inline    Polynomial<int> operator + 
  (const int& num, const Polynomial<int>& p2)
  { return (Polynomial<int>(num) + p2); }
  inline    Polynomial<int> operator - 
  (const int& num, const Polynomial<int>& p2)
  { return (Polynomial<int>(num) - p2); }
  inline    Polynomial<int> operator * 
  (const int& num, const Polynomial<int>& p2)
  { return (Polynomial<int>(num) * p2); }
  inline    Polynomial<int> operator / 
  (const int& num, const Polynomial<int>& p2)
  { return (Polynomial<int>(num)/p2); }
  inline    Polynomial<int> operator % 
  (const int& num, const Polynomial<int>& p2)
  { return (Polynomial<int>(num)%p2); }

  // righthand side
  inline    Polynomial<int> operator + 
  (const Polynomial<int>& p1, const int& num)
  { return (p1 + Polynomial<int>(num)); }
  inline    Polynomial<int> operator - 
  (const Polynomial<int>& p1, const int& num)
  { return (p1 - Polynomial<int>(num)); }
  inline    Polynomial<int> operator * 
  (const Polynomial<int>& p1, const int& num)
  { return (p1 * Polynomial<int>(num)); }
  inline    Polynomial<int> operator / 
  (const Polynomial<int>& p1, const int& num)
  { return (p1 / Polynomial<int>(num)); }
  inline    Polynomial<int> operator % 
  (const Polynomial<int>& p1, const int& num)
  { return (p1 % Polynomial<int>(num)); }

  // lefthand side
  inline    bool operator ==  
  (const int& num, const Polynomial<int>& p) 
  { return ( (Polynomial<int>(num)-p).sign() == CGAL::ZERO );}
  inline    bool operator != 
  (const int& num, const Polynomial<int>& p) 
  { return ( (Polynomial<int>(num)-p).sign() != CGAL::ZERO );}
  inline    bool operator <  
  (const int& num, const Polynomial<int>& p) 
  { return ( (Polynomial<int>(num)-p).sign() == CGAL::NEGATIVE );}
  inline    bool operator <=  
  (const int& num, const Polynomial<int>& p) 
  { return ( (Polynomial<int>(num)-p).sign() != CGAL::POSITIVE );}
  inline    bool operator >  
  (const int& num, const Polynomial<int>& p) 
  { return ( (Polynomial<int>(num)-p).sign() == CGAL::POSITIVE );}
  inline    bool operator >=  
  (const int& num, const Polynomial<int>& p) 
  { return ( (Polynomial<int>(num)-p).sign() != CGAL::NEGATIVE );}

  // righthand side
  inline    bool operator ==
  (const Polynomial<int>& p, const int& num) 
  { return ( (p-Polynomial<int>(num)).sign() == CGAL::ZERO );}
  inline    bool operator !=
  (const Polynomial<int>& p, const int& num) 
  { return ( (p-Polynomial<int>(num)).sign() != CGAL::ZERO );}
  inline    bool operator < 
  (const Polynomial<int>& p, const int& num) 
  { return ( (p-Polynomial<int>(num)).sign() == CGAL::NEGATIVE );}
  inline    bool operator <= 
  (const Polynomial<int>& p, const int& num) 
  { return ( (p-Polynomial<int>(num)).sign() != CGAL::POSITIVE );}
  inline    bool operator > 
  (const Polynomial<int>& p, const int& num) 
  { return ( (p-Polynomial<int>(num)).sign() == CGAL::POSITIVE );}
  inline    bool operator >=
  (const Polynomial<int>& p, const int& num) 
  { return ( (p-Polynomial<int>(num)).sign() != CGAL::NEGATIVE );}