simple_interval_root.h

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

// Copyright (c) 2005  Stanford University (USA).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; version 2.1 of the License.
// See the file LICENSE.LGPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Kinetic_data_structures/include/CGAL/Polynomial/internal/Simple_interval_root.h $
// $Id: Simple_interval_root.h 36369 2007-02-16 02:46:23Z drussel $
//
//
// Author(s)     : Daniel Russel <[e-mail]>

#ifndef CGAL_POLYNOMIAL_SIMPLE_INTERVAL_ROOT_H
#define CGAL_POLYNOMIAL_SIMPLE_INTERVAL_ROOT_H
#include <CGAL/Polynomial/basic.h>
#include <CGAL/Real_embeddable_traits.h>
#include <vector>
CGAL_POLYNOMIAL_BEGIN_INTERNAL_NAMESPACE;

//! A root represented as a bounding interval and a polynomial.
/*!
  Representing an interval which contains one root of a function.
  \todo cache sturm sequence
*/
template <class Traits>
class Simple_interval_root
{
  typedef Simple_interval_root<Traits> This;


  typedef enum Fields {INVALID= 4,UP=1, INF=2} Type_fields;
  typedef unsigned char Type;
  typedef typename Traits::Function Polynomial;
  typedef typename Traits::FT NT;
public:
  Simple_interval_root(){
    set_type(INVALID);
    CGAL_Polynomial_assertion(is_null());
  }

  Simple_interval_root(double nt): function_(0), interval_(nt,nt) {
    if (std::numeric_limits<double>::has_infinity && (nt == std::numeric_limits<double>::infinity()
						   || -nt == std::numeric_limits<double>::infinity())) {
      if (nt == std::numeric_limits<double>::infinity() ) {
	set_type(INF|UP);
      }
      else {
	set_type(INF);
      }
    } else {
      set_type(0);
      ii_= std::make_pair(NT(nt),NT(nt));
    }
    audit();
    compute_approximation();
  }

  Simple_interval_root(const Simple_interval_root<Traits> &o):ii_(o.ii_), type_(o.type_), 
							      function_(o.function_),
							      kernel_(o.kernel_), interval_(o.interval_)
#ifndef NDEBUG
							     , approximation_(o.approximation_)
#endif
  {
    
  }

  /*Simple_interval_root(Type t): type_(t) {
    audit();
    compute_approximation();
    };*/

  //! represent a rational root
  Simple_interval_root(const NT &nt): ii_(std::make_pair(nt,nt)), function_(0), interval_(0,-1) {
    set_type(0);
    audit();
    compute_approximation();
  }


  //! Represent a root by an interval and a polynomial
  Simple_interval_root(const std::pair<NT,NT> &ii,
		       const Polynomial &sa,
		       Sign slb, Sign,
		       Traits k): ii_(ii), function_(sa),
				  kernel_(k), interval_(0,-1){
    //CGAL_Polynomial_precondition(!ii_.is_singular());
    //Sign slb= sign_(ii_.lower_bound());
    if (slb == CGAL::NEGATIVE) {
      set_type(UP);
    } else {
      set_type(0);
    }
    compute_approximation();
    audit();
  }

  static This infinity() {
    static This ret(std::numeric_limits<double>::infinity());
    CGAL_Polynomial_postcondition(ret.is_infinite());
    //ret.compute_approximation();
    return ret;
  }

  bool operator<(const This &o) const
  {
    CGAL_Polynomial_expensive_precondition(!is_null() && !o.is_null());
    Comparison_result r= compare(o);
    audit(); o.audit();
    return r==SMALLER;
  }

  bool operator>(const This &o) const
  {
    CGAL_Polynomial_expensive_precondition(!is_null() && !o.is_null());
    Comparison_result r= compare(o);
    audit(); o.audit();
    return r==LARGER;
  }

  bool operator!=(const This &o) const
  {
    if (is_null()) return !o.is_null();
    else if (o.is_null()) return true;
    CGAL_Polynomial_expensive_precondition(!is_null() && !o.is_null());
    Comparison_result r= compare(o);
    audit(); o.audit();
    return r != EQUAL;
  }
  bool operator<=(const This &o) const
  {
    CGAL_Polynomial_expensive_precondition(!is_null() && !o.is_null());
    Comparison_result r= compare(o);
    audit(); o.audit();
    return r != LARGER;
  }
  bool operator>=(const This &o) const
  {
    CGAL_Polynomial_expensive_precondition(!is_null() && !o.is_null());
    Comparison_result r= compare(o);
    audit(); o.audit();
    return r != SMALLER;
  }

  bool operator==(const This &o) const
  {
    if (is_null()) return o.is_null();
    else if (o.is_null()) return false;
    CGAL_Polynomial_expensive_precondition(!is_null() && !o.is_null());
    Comparison_result r= compare(o);
    audit(); o.audit();
    return r==EQUAL;
  }

 
  const std::pair<NT, NT>& isolating_interval() const
  {
    CGAL_Polynomial_precondition(!is_infinite());
    CGAL_Polynomial_expensive_precondition(!is_null());
    return ii_;
  }

  //! Write lots of info about the interval
  void write(std::ostream &o) const
  {
#ifndef NDEBUG
    This t=*this;
    t.write_internal(o);
#else
    write_internal(o);
#endif
  }

  void print() const
  {
    write(std::cout);
  }

  //! Negate the interval.
  This operator-() const
  {
    CGAL_Polynomial_expensive_precondition(!is_null());
    if (is_pos_inf()) {
      This ret;
      ret.type_=INF;
      return ret;
    } else if (is_neg_inf()) return infinity();
    else {
      This copy= *this;
      copy.ii_= std::make_pair(-ii_.second, -ii_.first);
      typename Traits::Negate_variable nf= kernel_.negate_variable_object();
      copy.function_= nf(function_);
      if (type_&UP) {
	copy.type_= 0;
      } else {
	copy.type_=UP;
      }
      return copy;
    }
  }

  //! Return true if the root is +/- infinity.
  bool is_infinite() const
  {
    return type_&INF;
  }
 
  
  Comparison_result compare(const This &o) const
  {
    audit();
    o.audit();
    CGAL_Polynomial_precondition(-SMALLER == LARGER);
    if (is_infinite() || o.is_infinite()) {
      if (is_pos_inf() && o.is_pos_inf()) return EQUAL;
      if (is_neg_inf() && o.is_neg_inf()) return EQUAL;
      else if (is_pos_inf() || o.is_neg_inf()) return LARGER;
      else if (is_neg_inf() || o.is_pos_inf()) return SMALLER;
      else {
	CGAL_Polynomial_assertion(0); return EQUAL;
      }
    } else {
      CGAL::Comparison_result cmp;
      if (try_compare(o, cmp)) return cmp;
      else return compare_finite(o);
    }
  }
 std::pair<double, double> compute_interval(double accuracy=.0001) const
  {
    if (interval_.first > interval_.second) {
      //std::cout << "Computing interval for ";
      //write_raw(std::cout) << std::endl;
      std::pair<double,double> r= internal_compute_interval(accuracy);
      /*std::cout << "Got " << CGAL::to_interval(ii_.first).first << "..." 
	<< CGAL::to_interval(ii_.second).second << std::endl;*/
      interval_=std::make_pair(CGAL::to_interval(ii_.first).first, CGAL::to_interval(ii_.second).second);
    } 
    return interval_;
  }

  double compute_double(double accuracy) const
  {
    if (is_infinite()) {
      if (is_up()) {
	return double_inf_rep();
      }
      else {
	return -double_inf_rep();
      }
    }
    //This t= *this;
    std::pair<double, double> i= compute_interval(accuracy);
    return (i.first+i.second)/2.0;
  }


  NT rational_between(const This &o) const {
    CGAL_precondition(CGAL::compare(*this,o) == CGAL::SMALLER);
    CGAL_precondition(ii_.first != ii_.second || ii_.second != o.ii_.first || o.ii_.first != o.ii_.second);
    compare(o);
    CGAL_precondition(ii_.first != ii_.second || ii_.second != o.ii_.first || o.ii_.first != o.ii_.second);
    //write_internal(std::cout) << std::endl;
    //write_internal(std::cout) << std::endl;
    if (is_neg_inf()) return o.ii_.first -1;
        
    // hopefully disjoint
    NT ret = ii_.second;
    NT step= NT(.0000000596046447753906250000000);
    
    do {
      while (This(ret) <= *this) {
	ret+= step;
      }
      while (This(ret) >= o) {
	ret -= step;
      }
      step= step/NT(2.0);
      /*std::cout << ret << "  (" << step << ")" 
		<< o.ii_.first - ii_.second 
		<< " " << o.ii_.second- ii_.first << std::endl;*/
    } while (This(ret) >= o || This(ret) <= *this);
    return ret;
  }

protected:
  std::pair<double, double> internal_compute_interval(double accuracy) const
  {
 
    if (type_&INF) {
      if (is_up()) {
	return std::pair<double, double>(double_inf_rep(), double_inf_rep());
      }
      else {
	return std::pair<double, double>(-double_inf_rep(), -double_inf_rep());
      }
    }

    //double oaw;                           //= ii_.approximate_width();
    // int ct=0;
    while (ii_.second != ii_.first && ii_.second-ii_.first > NT(accuracy)) {
      refine();
      /*++ct;
      if (ct== 30) {
	std::cerr << "Error subdividing ";
	write_internal(std::cerr) << std::endl;
	break;
	}*/
    }
    return std::make_pair(CGAL::to_interval(ii_.first).first, CGAL::to_interval(ii_.second).second);
  }
  std::ostream& write_internal(std::ostream &o) const
  {
    if (is_pos_inf()) {
      o << "inf";
    }
    else if (is_neg_inf()) {
      o << "-inf";
    }
    else {
      if (ii_.first == ii_.second) {
	o << ii_.first;
      } else {
	o << function_ << " in [" << ii_.first << "," << ii_.second << "]";
	o << " = " << internal_compute_interval(.00001).first 
	  << "..." << internal_compute_interval(.00001).second;
      }
    }
    return o;
  }

  std::ostream & write_raw(std::ostream &o) const
  {
    if (is_pos_inf()) {
      o << "inf";
    }
    else if (is_neg_inf()) {
      o << "-inf";
    }
    else {
      if (ii_.first == ii_.second) {
	o << ii_.first;
      } else {
	o << function_ << " in [" << ii_.first << "," << ii_.second << "]";
      }
    }
    return o;
  }

  void set_type(Type t)
  {
    type_=t;
  }

  bool is_pos_inf() const
  {
    return (type_&INF) && (type_&UP);
  }

  bool is_neg_inf() const
  {
    return (type_&INF) && !(type_&UP);
  }