sturm_root_rep.h

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

    //CGAL_assertion( this_line_should_not_have_been_reached );

    return CGAL::EQUAL;
  }                                         // end-compare_intersecting

  //--------------------------------------------------------------
  void subdivide() const
  {
    if (!refined_) {
      refined_=true;
      //std::cout << "Refining " << *this << " " << std::endl;
    }


    Interval_pair ivl_pair = ivl.split();
    Interval first_half  = ivl_pair.first;
    Interval second_half = ivl_pair.second;

    Sign_at_functor sign_at_f(this);

    Sign s_mid =
      first_half.apply_to_endpoint( sign_at_f, Interval::UPPER );

    if ( s_mid == CGAL::ZERO ) {
      ivl = first_half.endpoint_interval( Interval::UPPER );
      return;
    }

    if ( s_mid == sign_upper() ) {
      ivl = first_half;
    }
    else {
      ivl = second_half;
    }
  }

  

  template<class Child>
  Polynomial compute_simple(Field_tag, const Child&) const
  {
    Polynomial gcd = sseq.exact( sseq.exact_size() - 1 );
    return p_ / gcd;
  }

  /*Polynomial compute_simple(Integral_domain_without_division_tag, const Self&) const
    {
    Polynomial gcd = sseq[sseq.size() - 1];
    return tr_.pseudo_quotient_object()(p_, gcd);
    }*/
  
  Polynomial compute_simple(Field_tag, const Self&) const
  {
    Polynomial gcd = sseq[sseq.size() - 1];
    return tr_.quotient_object()(p_, gcd);
  }

  template<class This>
  Comparison_result
  //    compare_same_interval(const Self& r) const
  compare_same_interval(const This& r) const
  {
    CGAL_precondition( lower_bound() == r.lower_bound() );
    CGAL_precondition( upper_bound() == r.upper_bound() );

    // here I have to do some subdivisions...

    Method_tag mtag;
    Sign_Sturm_sequence seq
      = tr_.sign_Sturm_sequence_object(r.compute_simple(mtag, r),
				       compute_simple(mtag, r) );

    int sign_of_r_at_this;

    sign_of_r_at_this =
      seq.sum_of_signs(lower_bound(), upper_bound());

    if ( sign_of_r_at_this == 0 ) { return CGAL::EQUAL; }

    Sign s_r;
    if ( sign_of_r_at_this < 0 ) {
      s_r = CGAL::NEGATIVE;
    }
    else {
      s_r = CGAL::POSITIVE;
    }

    Sign s_l = sign_lower();

    return ( s_r != s_l ) ? CGAL::SMALLER : CGAL::LARGER;
  }

  //===================
  // SIGNS EVALUATORS
  //===================
  Sign sign_lower() const
  {
    return s_lower;
  }

  Sign sign_upper() const
  {
    return s_upper;
  }

  template<class This>
  Sign sign_at(const This& r, typename Interval::Endpoint b) const
  {
    Sign_at_functor sign_at_f(this);

    return r.ivl.apply_to_endpoint(sign_at_f, b);
  }

public:
  //==================================
  // POSITIVE AND NEGATIVE INFINITY
  //==================================
  static Self infinity(const Traits& tr = Traits()) {
    return Self(-1, 0, tr);
  }

  //================
  // CONSTRUCTORS
  //================
protected:
  Sturm_root_rep(int i, int, const Traits& tr)
    : idx(i), ivl(), p_(), sseq(), s_lower(CGAL::ZERO),
      s_upper(CGAL::ZERO), multiplicity_(0), tr_(tr), refined_(false)  {}

public:
  Sturm_root_rep(const Traits& tr = Traits())
    : idx(-1), ivl(), p_(), sseq(), s_lower(CGAL::ZERO),
      s_upper(CGAL::ZERO), multiplicity_(0), tr_(tr), refined_(false)  {}
  //-------
  template<class T>
  Sturm_root_rep(const T& a, const Traits& tr = Traits())
    : idx(1), ivl(a), p_(), sseq(), s_lower(CGAL::ZERO),
      s_upper(CGAL::ZERO), multiplicity_(0), tr_(tr), refined_(false)  {}
  //-------
  Sturm_root_rep(const Interval& ivl,
		 const Polynomial& p,
		 const Standard_sequence& sseq, int idx,
		 const Traits& tr)
    : idx(idx), ivl(ivl), p_(p), sseq(sseq),
      s_lower(CGAL::ZERO), s_upper(CGAL::ZERO), multiplicity_(0),
      tr_(tr), refined_(false) {
    //++sturm_created__;
    Sign_at_functor sign_at_p(this);
    s_lower = apply(sign_at_p, ivl.lower_bound());
    s_upper = apply(sign_at_p, ivl.upper_bound());

    //      s_lower = sign_at( *this, Interval::LOWER );
    //      s_upper = sign_at( *this, Interval::UPPER );
  }
  //-------
  Sturm_root_rep(const Self& other)
    : idx(other.idx), ivl(), p_(other.p_), sseq(other.sseq),
      s_lower(other.s_lower),   s_upper(other.s_upper),
      multiplicity_(other.multiplicity_), tr_(other.tr_), refined_(other.refined_){
    if ( other.idx >= 1 ) {
      ivl = other.ivl;
    }
  }

  //=======================
  // ASSIGNMENT OPERATOR
  //=======================
  const Self& operator=(const Self& other) {
    if ( &other == this ) { return *this; }
    idx = other.idx;
    p_ = other.p_;
    sseq = other.sseq;
    s_lower = other.s_lower;
    s_upper = other.s_upper;
    multiplicity_ = other.multiplicity_;
    tr_ = other.tr_;

    if ( idx >= 1 ) {
      ivl = other.ivl;
    }
    return *this;
  }

  //=========================
  // COMPARISON OPERATORS
  //=========================
  bool operator!=(const Self& r) const
  {
    return compare(r) != CGAL::EQUAL;
  }

  bool operator==(const Self& r) const
  {
    return compare(r) == CGAL::EQUAL;
  }

  bool operator<=(const Self& r) const
  {
    return compare(r) != CGAL::LARGER;
  }

  bool operator>=(const Self& r) const
  {
    return compare(r) != CGAL::SMALLER;
  }

  bool operator>(const Self& r) const
  {
    return compare(r) == CGAL::LARGER;
  }

  bool operator<(const Self& r) const
  {
    return compare(r) == CGAL::SMALLER;
  }

  Self operator-() const
  {
    if ( idx == -1 ) {
      return Self(-2, 0, tr_);
    }
    else if ( idx == -2 ) {
      return Self(-1,0, tr_);
    }
    else if (idx == 1) {
      // HACK by Daniel
      return Self(-ivl, tr_);
    }
    else {

      Interval m_ivl = -ivl;
      Polynomial m_p = tr_.negate_variable_object()(p_);
      Standard_sequence m_sseq= tr_.standard_sequence_object(m_p);

      Self m_root(m_ivl, m_p, m_sseq, idx, tr_);

      m_root.multiplicity_ = multiplicity_;

      return m_root;
    }
  }

  //===========================
  // THE MULTIPLICITY METHOD
  //===========================
  unsigned int multiplicity() const
  {
    if ( idx < 0 ) { return 0; }
    if ( multiplicity_ == 0 ) {
      compute_multiplicity();
    }last_zero_=false;
    if (multiplicity_%2==0) return multiplicity_/2;
    else return multiplicity_;
  }

  //====================
  // THE PARITY METHOD
  //====================
  bool is_even_multiplicity() const
  {
    return false;
    /*if ( idx < 0 ) { return false; }

    if ( is_exact() ) {
      return (multiplicity() % 2 == 0);
    }
    else {
      Sign_at_functor sign_at_f(this);

      Sign sl = apply(sign_at_f, ivl.lower_bound());
      Sign su = apply(sign_at_f, ivl.upper_bound());
      return (sl == su);
      }*/
  }

  //==========================
  // ACCESS TO THE INTERVAL
  //==========================
  typename Interval::NT lower_bound() const { 
    CGAL_precondition(idx >= 0);
    return ivl.lower_bound(); }
  typename Interval::NT upper_bound() const { 
    CGAL_precondition(idx >= 0);
    return ivl.upper_bound(); }
  Interval    interval()    const { return ivl; }
  int         index()       const { return idx; }
  const Polynomial&     polynomial() const { return p_; }
  //    Polynomial  simple_polynomial() const { return p; }

public:
  //======================
  // CONVERTOR TO double
  //======================
  double compute_double(double acc = 1e-10) const
  {
    if (idx < 0) {
      double inf=std::numeric_limits<double>::has_infinity? std::numeric_limits<double>::infinity() : (std::numeric_limits<double>::max)();
      if ( idx == -1 ){
	return inf;
      } else return -inf;
    }
    /*if ( idx == -2 ) {
      return -DBL_MAX * DBL_MAX;
      }*/

    if ( is_exact() ) {
      return CGAL_POLYNOMIAL_TO_DOUBLE(lower_bound());
    }

    Exact_nt xacc(acc);


    while ( ivl.approximate_width() > acc ) {
      subdivide();
      /*if (!refined_) {
	++sturm_refined__;
	refined_=true;
	std::cout << "Refining " << *this << " to compute double " << std::endl;
	}*/
      if ( is_exact() ) { break; }
    }

    return CGAL_POLYNOMIAL_TO_DOUBLE( lower_bound() );
  }

  const std::pair<NT, NT>& isolating_interval() const {
    CGAL_precondition(idx >=0);
    return ivl.to_pair();
  }

  //========================
  // CONVERTOR TO interval
  //========================
  std::pair<double,double> compute_interval() const
  {
    if (*this == infinity()){ 
      return std::make_pair(std::numeric_limits<double>::infinity(),
			    std::numeric_limits<double>::infinity());
    }
    else if (*this == -infinity()) {
      return std::make_pair(-std::numeric_limits<double>::infinity(),
			    -std::numeric_limits<double>::infinity());
    } else {
      compute_double();
	    
      std::pair<double,double> i_low =
	CGAL_POLYNOMIAL_TO_INTERVAL(ivl.lower_bound());
      std::pair<double,double> i_high =
	CGAL_POLYNOMIAL_TO_INTERVAL(ivl.upper_bound());
	    
      return std::pair<double,double>(i_low.first, i_high.second);
    }
  }

  //================
  // STREAM WRITER
  //================
  template<class Stream>
  Stream& write(Stream& os) const
  {
    if ( idx == -2 ) {
      os << "-inf";
    }
    else if ( idx == -1 ) {
      os << "inf";
    }
    else {
      os << "{" << ivl << ", " << idx << "}";
    }
    if (idx != -2 && idx != -1) {
      Self copy = *this;
      os << " = " << copy.compute_double();
    }
    return os;
  }

protected:
  bool is_up_;
  int  
  int                    idx;
  mutable Interval       ivl;
  Polynomial             p_;
  Standard_sequence      sseq;
  mutable Sign           s_lower, s_upper;
  mutable unsigned int   multiplicity_;
  Traits                 tr_;
  mutable bool refined_;
};

template<class Stream, class S, class I>
Stream&
operator<<(Stream& os, const Sturm_root_rep<S,I>& r)
{
  return r.write(os);
}

/*template<class S, class I>
  std::pair<double,double>
  to_interval(const Sturm_root_rep<S,I>& r)
  {
  return r.compute_interval();
  }*/

CGAL_POLYNOMIAL_END_INTERNAL_NAMESPACE

CGAL_BEGIN_NAMESPACE

/*template<class S, class I>
  double
  to_double(const CGAL_POLYNOMIAL_NS::internal::Sturm_root_rep<S,I>& r)
  {
  return r.compute_double();
  }


  template<class S, class I>
  std::pair<double,double>
  to_interval(const CGAL_POLYNOMIAL_NS::internal::Sturm_root_rep<S,I>& r)
  {
  return r.compute_interval();
  }*/


template <class T, class I>
class Real_embeddable_traits< CGAL::POLYNOMIAL::internal::Sturm_root_rep<T,I> > 
  : public Real_embeddable_traits_base< CGAL::POLYNOMIAL::internal::Sturm_root_rep<T,I> > {
public:
  typedef CGAL::POLYNOMIAL::internal::Sturm_root_rep<T,I>  Type;
  class Abs 
    : public Unary_function< Type, Type > {
  public:
    Type operator()( const Type& x ) const {
      if (x < Type(0)) return -x;
      else return x;
    }
  };
    
  class Sign 
    : public Unary_function< Type, ::CGAL::Sign > {
  public:
    ::CGAL::Sign operator()( const Type& x ) const {
      return static_cast<CGAL::Sign>(x.compare(0));
    }        
  };
    
  class Compare 
    : public Binary_function< Type, Type,
			      Comparison_result > {
  public:
    Comparison_result operator()( const Type& x, 
				  const Type& y ) const {
      return x.compare(y);
    }
        
    CGAL_IMPLICIT_INTEROPERABLE_BINARY_OPERATOR_WITH_RT( Type,
							 Comparison_result )
        
      };
    
  class To_double 
    : public Unary_function< Type, double > {
  public:
    double operator()( const Type& x ) const {
      // this call is required to get reasonable values for the double
      // approximation
      return x.compute_double();
    }
  };
    
  class To_interval 
    : public Unary_function< Type, std::pair< double, double > > {
  public:
    std::pair<double, double> operator()( const Type& x ) const {

      return x.compute_interval();
    }          
  };
};


CGAL_END_NAMESPACE

namespace std
{
  template <class S, class I>
  class numeric_limits<CGAL_POLYNOMIAL_NS::internal::Sturm_root_rep<S,I> >:
    public numeric_limits<typename CGAL_POLYNOMIAL_NS::internal::Sturm_root_rep<S,I>::NT >
  {
  public:
    typedef numeric_limits<typename CGAL_POLYNOMIAL_NS::internal::Sturm_root_rep<S,I>::NT > P;
    typedef CGAL_POLYNOMIAL_NS::internal::Sturm_root_rep<S,I> T;
    static const bool is_specialized = true;
    static T min BOOST_PREVENT_MACRO_SUBSTITUTION () throw() {return T((P::min)());}
    static T max BOOST_PREVENT_MACRO_SUBSTITUTION () throw() {return T((P::max)());}
    /*static const int digits =0;
      static const int digits10 =0;
      static const bool is_signed = true;
      static const bool is_integer = false;
      static const bool is_exact = true;
      static const int radix =0;
      static T epsilon() throw(){return T(0);}
      static T round_error() throw(){return T(0);}
      static const int min_exponent=0;
      static const int min_exponent10=0;
      static const int max_exponent=0;
      static const int max_exponent10=0;*/
    static const bool has_infinity=true;
    /*static const bool has_quiet_NaN = false;
      static const bool has_signaling_NaN= false;
      static const float_denorm_style has_denorm= denorm_absent;
      static const bool has_denorm_loss = false;
    */
    static T infinity() throw() {return T::infinity();}
    /*static T quiet_NaN() throw(){return T(0);}
      static T denorm_min() throw() {return T(0);}
      static const bool is_iec559=false;
      static const bool is_bounded =false;
      static const bool is_modulo= false;
      static const bool traps = false;
      static const bool tinyness_before =false;
      static const float_round_style round_stype = round_toward_zero;*/
  };
};
#endif                                            // CGAL_POLYNOMIAL_STURM_ROOT_REP_H