simple_interval_root.h

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

  bool contains_root(Sign ls, Sign us) const
  {
    return (is_up() && ls==NEGATIVE && us==POSITIVE) || (is_down() && ls==POSITIVE && us==NEGATIVE);
  }

  void refine() const
  {
    if (ii_.first == ii_.second) return;
    //CGAL_Polynomial_precondition(!is_rational());
    //CGAL_Polynomial_precondition(is_normal());
    NT mp = (ii_.first + ii_.second)*NT(0.5);
    Sign sn= sign_at(mp);
    if (sn == ZERO) {
      ii_= std::make_pair(mp,mp);
    } else if (contains_root(lower_sign(), sn)) {
      ii_.second= mp;
    } else {
      ii_.first= mp;
    }
  }

  void refine_using(const std::pair<NT, NT> &o) const
  {
    if (ii_.first== ii_.second) return;
    std::vector<NT > plist;
    plist.push_back(ii_.first);
    if (o.first < ii_.second && o.first > ii_.first) {
      plist.push_back(o.first);
    }
    if (o.first != o.second && o.second < ii_.second && o.second > ii_.first) {
      plist.push_back(o.second);
    }
    plist.push_back(ii_.second);

    /*for (unsigned int i=0; i< plist.size(); ++i) {
      std::cout << plist[i] << "   ";
    }
    std::cout << std::endl;*/

    if (plist.size()==2) return;
    
    CGAL::Sign ps= lower_sign();
    //std::cout << "ps is " << ps << std::endl;
    CGAL_assertion(ps != CGAL::ZERO);
    for (unsigned int i=1; i< plist.size()-1; ++i) {
      CGAL::Sign sn= sign_at(plist[i]);
      //std::cout << "sn is " << ps << std::endl;
      if (sn==0) {
	ii_= std::make_pair(plist[i], plist[i]);
	audit();
	return;
      } else if (sn != ps) {
	ii_= std::make_pair(plist[i-1], plist[i]);
	audit();
	return;
      }
    }

    ii_= std::make_pair(plist[plist.size()-2], plist[plist.size()-1]);
    CGAL_postcondition(sign_at(plist[plist.size()-2]) == ps);
    CGAL_postcondition(sign_at(plist[plist.size()-1]) == CGAL::Sign(-ps));
    audit();
  }

  Comparison_result compare_finite(const This &o) const
  {
    audit();
    o.audit();

    refine_using(o.ii_);
    o.refine_using(ii_);

    do {
      audit(); o.audit();

      CGAL::Comparison_result cmp;
      if (try_compare(o, cmp)) return cmp;
    
      refine();
      o.refine();
    } while (ii_.second-ii_.first  > NT(.0000001));

    //std::cout << "Using sturm to compare " << *this << " and " << o << std::endl;

    typename Traits::Compare_isolated_roots_in_interval pred=kernel_.compare_isolated_roots_in_interval_object(function_,o.function_);
    Comparison_result co= pred(ii_.first, ii_.second);
    //std::cout << "The result is " << co << std::endl;
    return co;
  }

  //! Check that everything is correct
  void audit() const
  {
    CGAL_Polynomial_assertion(!is_null());
#ifndef NDEBUG
    bool problem=false;
    if (type_&INF) {

    } else if (ii_.first == ii_.second) {
      
    } else {
      if (is_up()) {
	problem = problem || (sign_at(ii_.first) != NEGATIVE);
	problem = problem || (sign_at(ii_.second) != POSITIVE);
      } else {
	problem = problem || (sign_at(ii_.first) != POSITIVE);
	problem = problem || (sign_at(ii_.second) != NEGATIVE);
      }
    }
    if (problem) {
      std::cerr << "Problem with interval.\n";
      std::cerr << "Type is " << type_ << std::endl;
      std::cerr << ii_.first << "..." << ii_.second << ", " << sign_at(ii_.first) 
		<< sign_at(ii_.second) << std::endl;
      CGAL_Polynomial_exactness_assertion(0);
    }
#endif
  }



  CGAL::Sign sign_at(const NT &nt) const
  {
    return kernel_.sign_at_object()(function_, nt);
  }

  //! The representation of negative infinity
  /*static This negative_infinity(){
    This ret(INF);
    //ret.set_type(INF);
    //ret.compute_approximation();
    return ret;
    }*/

  bool try_compare(const This &o, CGAL::Comparison_result &cmp) const {
    if (ii_.first == ii_.second){
      if (o.ii_.first == o.ii_.second) {
	cmp= CGAL::compare(ii_.first, o.ii_.second);
	return true;
      } else {
	if (o.ii_.first < ii_.first && o.ii_.second > ii_.first) {
	  return false;
	} else {
	  cmp= CGAL::compare(ii_.first, o.ii_.first);
	  if (cmp == CGAL::EQUAL){
	    cmp=  CGAL::compare(ii_.first, o.ii_.second);
	  }
	  //CGAL_assertion(CGAL::compare(ii_.second, o.ii_.second) == cmp);
	  return true;
	}
      }
    } else if (o.ii_.first == o.ii_.second) {
      if (ii_.first < o.ii_.first && ii_.second > o.ii_.first) {
	  return false;
	} else {
	  cmp= CGAL::compare(ii_.first, o.ii_.first);
	  if (cmp == CGAL::EQUAL) {
	    cmp=  CGAL::compare(ii_.second, o.ii_.first);
	  }
	  //CGAL_assertion(CGAL::compare(ii_.second, o.ii_.second) == cmp);
	  return true;
	}
    } else {
      if (ii_.first >= o.ii_.second) {
	cmp= CGAL::LARGER;
	return true;
      } else if (ii_.second <= o.ii_.first) {
	cmp= CGAL::SMALLER;
	return true;
      } else return false;
    }
  }

  bool is_up() const
  {
    return type_&UP;
  }
  bool is_down() const
  {
    return !(type_&UP);
  }

  static double double_inf_rep() {
    if (std::numeric_limits<double>::has_infinity) {
      return (std::numeric_limits<double>::infinity());
    } else return ((std::numeric_limits<double>::max)());
  }

  Sign lower_sign() const
  {
    CGAL_Polynomial_precondition(!is_infinite());
    if (is_up()) return NEGATIVE;
    else return POSITIVE;
  }

  //! comptute a value that can be inspected in the compiler
  void compute_approximation() {
#ifndef NDEBUG
    approximation_=immutable_double_approximation(0.0001);
#endif
  }

  double immutable_double_approximation(double accuracy=0.00001) const
  {
    CGAL_Polynomial_expensive_precondition(!is_null());
    This temp = *this;
    return temp.internal_compute_interval(accuracy).first;
  }

  //! return true if the this is uninitialized
  bool is_null() const
  {
    return (type_&INVALID);
  }

  mutable std::pair<NT, NT> ii_;
  //Function f_;
  Type type_;
  Polynomial function_;
  Traits kernel_;
  mutable std::pair<double,double> interval_;
#ifndef NDEBUG
  double approximation_;
#endif
};

/*
  template <class F>
  double to_double(const Simple_interval_root<F> &f){
  return f.to_double();
  }

  template <class F>
  std::pair<double,double> to_interval(const Simple_interval_root<F> &f){
  return f.to_ii_;
  }*/

template <class F>
std::ostream &operator<<(std::ostream &out, const Simple_interval_root<F> &f)
{
  f.write(out);
  return out;
}

/*
template <class F>
bool root_is_even_multiplicity(const Simple_interval_root<F> &f)
{
  return f.is_even_multiplicity();
}


template <class F>
typename F::NT to_rational(const Simple_interval_root<F> &f)
{
  return f.to_rational();
}


template <class F>
bool is_rational(const Simple_interval_root<F> &f)
{
  return f.is_rational();
  }*/


CGAL_POLYNOMIAL_END_INTERNAL_NAMESPACE

CGAL_BEGIN_NAMESPACE


template <class T>
class Real_embeddable_traits< CGAL::POLYNOMIAL::internal::Simple_interval_root<T> > 
  : public Real_embeddable_traits_base< CGAL::POLYNOMIAL::internal::Simple_interval_root<T> > {
public:
  typedef CGAL::POLYNOMIAL::internal::Simple_interval_root<T>  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(.00000001);
    }
  };
    
  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(.00001);
    }          
  };
};



CGAL_END_NAMESPACE

namespace std
{
  template <class Tr>
  class numeric_limits<CGAL_POLYNOMIAL_NS::internal::Simple_interval_root<Tr> >
  {
  public:
    typedef CGAL_POLYNOMIAL_NS::internal::Simple_interval_root<Tr> T;
    static const bool is_specialized = true;
    static T min BOOST_PREVENT_MACRO_SUBSTITUTION () throw () {return -T::infinity();}
    static T max BOOST_PREVENT_MACRO_SUBSTITUTION () throw () {return T::infinity();}
    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 = true;
    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();}
    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