sturm_root_rep.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/Sturm_root_rep.h $
// $Id: Sturm_root_rep.h 35769 2007-01-21 23:40:54Z drussel $
// 
//
// Author(s)     : Daniel Russel <[e-mail]>

#ifndef CGAL_POLYNOMIAL_STURM_ROOT_REP_H
#define CGAL_POLYNOMIAL_STURM_ROOT_REP_H

#include <CGAL/Polynomial/basic.h>
#include <CGAL/Polynomial/internal/Rational/Sign_Sturm_sequence.h>
#include <limits>
#include <cfloat>

//#include <CGAL/Polynomial/internal/Bisection.h>


CGAL_POLYNOMIAL_BEGIN_INTERNAL_NAMESPACE

//==================
// the Root class
//==================
template<class Solver_t, class Interval_t>
class Sturm_root_rep
{
public:
  typedef Solver_t                            Solver;
  typedef Interval_t                          Interval;
  typedef typename Solver::Traits             Traits;

  //  typedef typename Solver::Exact_nt           Exact_nt;
  //  typedef typename Solver::Storage_function   Storage_function;
  //  typedef typename Solver::Exact_function     Exact_function;
  typedef typename Solver::Method_tag         Method_tag;

  //  typedef typename Solver::Function_handle    Function_handle;
  typedef typename Solver::Standard_sequence  Standard_sequence;

  typedef typename Solver::Polynomial         Polynomial;
  typedef typename Solver::NT                 NT;
  typedef typename Solver::Sign_at            Sign_at;
  typedef NT                                  Exact_nt;

  typedef Sturm_root_rep<Solver,Interval>  Self;

  class Sign_at_functor
  {
  public:
    typedef Sign result_type;

    Sign_at_functor(const Self* outer) : outer(outer) {}

    template<class T>
    result_type operator()(const T& x) const
    {

      Sign s1 = outer->sseq.sign_at(x, 0);
      Sign s2 = outer->sseq.sign_at_gcd(x);

      CGAL_assertion( s1 == CGAL::ZERO || s2 != CGAL::ZERO );
      return Sign(s1 * s2);
    }

  private:
    const Self* outer;
  };

  typedef typename Traits::Sign_Sturm_sequence   Sign_Sturm_sequence;

protected:
  typedef std::pair<Interval,Interval>               Interval_pair;

public:
  //==========================================
  // METHODS FOR COMPUTING THE MULTIPLICITY
  //==========================================

  bool is_rational() const
  {
    return is_exact();
  }

  NT to_rational() const
  {
    CGAL_Polynomial_precondition( is_rational() );
    return ivl.lower_bound();
  }

  bool is_exact() const
  {
    return ivl.is_singular();
  }

  // computes the multiplicity using the interval number type;
  // if computation was successful the method sets success to true;
  // if computation was not possible the method sets success to
  // false;

  void compute_multiplicity() const
  {
    CGAL_precondition( multiplicity_ == 0 );

    if ( p_.degree() == 1 ) {
      multiplicity_ = 1;
      return;
    }

    if ( is_exact() ) {
      compute_multiplicity_exact();
    }
    else {
      compute_multiplicity_interval();
    }
  }

  void compute_multiplicity_exact() const
  {
    Polynomial q = p_;

    while ( true ) {
      multiplicity_++;
      q = tr_.differentiate_object()(q);

      if ( q.is_zero() ) { break; }

      Sign_at  sign_at_q(q);
      Sign s = ivl.apply_to_endpoint(sign_at_q, Interval::LOWER);
      if ( s != CGAL::ZERO ) {
	break;
      }
    }
  }

  void compute_multiplicity_interval() const
  {
    bool is_even_multiplicity_ = is_even_multiplicity();

    typename Traits::Differentiate differentiate =
      tr_.differentiate_object();

    Polynomial q = p_;
    bool first_time = true;
    while ( true ) {
      if ( first_time ) {
	first_time = false;
	if ( !is_even_multiplicity_ ) {
	  multiplicity_++;
	  q = differentiate(q);
	}
	else {
	  multiplicity_ += 2;
	  q = differentiate( differentiate(q) );
	}
      }
      else {
	multiplicity_ += 2;
	q = differentiate( differentiate(q) );
      }

      if ( q.degree() <= 0 ) { break; }

      Sign_Sturm_sequence sign_sturm =
	tr_.sign_Sturm_sequence_object(p_, q);

      int sign = ivl.apply_to_interval(sign_sturm);

      if ( sign != 0 ) { break; }
    }
  }

  //===========================================
  // HELPER METHODS FOR UPDATING THE VALUES
  // ASSOCIATED WITH THE INTERVAL
  //===========================================

  void set_lower(const typename Interval::NT& l,
		 const Sign& s_l) const
  {
    ivl.set_lower(l);
    s_lower = s_l;
  }

  void set_upper(const typename Interval::NT& u,
		 const Sign& s_u) const
  {
    ivl.set_upper(u);
    s_upper = s_u;
  }

  void set_interval(const typename Interval::NT& x) const
  {
    ivl.set_lower(x);
    ivl.set_upper(x);
    s_lower = s_upper = CGAL::ZERO;
  }



  template <class This>
  Comparison_result
  compare_finite(const This &r, bool subdiv=true) const {
    // consider the cases that the root is known exactly;
    // this is equivalent to saying that the two endpoints for the
    // interval containing the root are the same;
    // moreover if the two interval endpoints are not the same we
    // make sure that the interval endpoints are not roots.
    // this will make life easier afterwards.

    if ( is_exact() ) {
      //	std::cout << "first is exact" << std::endl;
      if ( r.is_exact() ) {
	//	  std::cout << "second is exact" << std::endl;
	return CGAL::compare( lower_bound(), r.lower_bound() );
      }
      else {
	if ( upper_bound() <= r.lower_bound() ) {
	  return CGAL::SMALLER;
	}
	else if ( lower_bound() >= r.upper_bound() ) {
	  return CGAL::LARGER;
	} else if ( upper_bound() > r.lower_bound() &&
                    lower_bound() < r.upper_bound() ) {
	  Sign s_r_lb = r.sign_lower();
	  {
	    Sign s_r_ub = r.sign_upper();
	    CGAL_assertion( s_r_lb != CGAL::ZERO && s_r_ub != CGAL::ZERO );
	    if(0) s_r_ub= CGAL::ZERO;
	  }

	  Sign s_at_r = r.sign_at( *this, Interval::LOWER );

	  if ( s_at_r == CGAL::ZERO ) { return CGAL::EQUAL; }
	  return ( s_at_r == s_r_lb ) ? CGAL::SMALLER : CGAL::LARGER;
	}
      }
    }
    else if ( r.is_exact() ) {

      // we have already checked the case that this root is also
      // known exactly
      //	std::cout << "second is exact" << std::endl;
      if ( upper_bound() <= r.lower_bound() ) {
	//	  std::cout << "this.upper bound <= other.lower bound" << std::endl;
	return CGAL::SMALLER;
      }
      else if ( lower_bound() >= r.upper_bound() ) {
	//	  std::cout << "this.lower bound >= other.upper bound" << std::endl;
	return CGAL::LARGER;
      } else if (  upper_bound() > r.lower_bound() &&
		   lower_bound() < r.upper_bound() ) {
	//	   std::cout << "other in interval of this" << std::endl;
	Sign s_ub = sign_upper();
	{
	  Sign s_lb = sign_lower();
	  CGAL_assertion( s_lb != CGAL::ZERO && s_ub != CGAL::ZERO );
	  if(0) s_lb= CGAL::ZERO;
	}

	Sign s_at_this = sign_at( r, Interval::LOWER );

	if ( s_at_this == CGAL::ZERO ) { return CGAL::EQUAL; }
	return ( s_ub == s_at_this ) ? CGAL::SMALLER : CGAL::LARGER;
      }
    }
    else {
      // now the roots are in the interior of interval
      // check if the interiors of the intervals are disjoint
      //	std::cout << "both not exact" << std::endl;
      //	int prec = std::cout.precision();
      //	std::cout.precision(16);
      //	std::cout << "this:  " << lower_bound() << " " << upper_bound()
      //		  << std::endl;
      //	std::cout << "other: " << r.lower_bound() << " "
      //		  << r.upper_bound() << std::endl;
      //	std::cout.precision(prec);

      if ( upper_bound() <= r.lower_bound() ) {
	//	  std::cout << "this.upper bound <= other.lower bound" << std::endl;
	return CGAL::SMALLER;
      }
      else if ( lower_bound() >= r.upper_bound() ) {
	//	  std::cout << "this.lower bound >= other.upper bound" << std::endl;
	return CGAL::LARGER;
      } else {
	if (subdiv) {
	  int count=0;
	  while ( upper_bound() > r.lower_bound() &&
		  lower_bound() < r.upper_bound() ) {
	    subdivide();
	    ++count;
	    if (count==4) break;
	    //	  std::cout << "intersecting intervals" << std::endl;
	  }
	  return compare_finite(r, false);
	} else {
	  return compare_intersecting(r);
	}
      }
    }

    bool this_line_should_not_have_been_reached = false;
    CGAL_assertion( this_line_should_not_have_been_reached );
    if (0) this_line_should_not_have_been_reached= false;
    return CGAL::EQUAL;
  }
  

  //===========================================
  // HELPER METHODS FOR COMPARISON OPERATORS
  //===========================================
  template<class This>
  Comparison_result compare(const This& r) const
  {
    // check against positive and negative infinity
    if (idx == -2) {
      if ( r.idx == -2 ) { return CGAL::EQUAL; }
      return CGAL::SMALLER;
    }
    if (idx == -1) {
      if ( r.idx == -1 ) { return CGAL::EQUAL; }
      return CGAL::LARGER;
    }
    if (r.idx == -2) {
      // I have already checked if both are negative infinity
      return CGAL::LARGER;
    }
    if (r.idx == -1) {
      // I have already checked if both are positive infinity
      return CGAL::SMALLER;
    }
    return compare_finite(r);
  }

  
  //--------------------------------------------------------------


  template<class This>
  Comparison_result
  compare_intersecting(const This& r) const
  {
    // in this method we know that the intervals are not
    // trivial, are intersecting and that the endpoints are not
    // roots

    // Case 1: the intervals have common left endpoint
    if ( lower_bound() == r.lower_bound() ) {
      if ( upper_bound() > r.upper_bound() ) {
	Sign s = sign_at( r, Interval::UPPER );

	if ( s == CGAL::ZERO ) {
	  // the root of this is equal to the upper bound of r
	  //	    set_interval( r, Interval::UPPER );
	  set_interval( r.upper_bound() );
	  return CGAL::LARGER;
	}
	else {
	  if ( s == sign_lower() ) {
	    return CGAL::LARGER;
	  }
	  else {
	    //	      set_upper( r, Interval::UPPER, s );
	    set_upper( r.upper_bound(), s );
	    return compare_same_interval(r);
	  }
	}
      }
      else if ( r.upper_bound() > upper_bound() ) {
	const This* new_this = static_cast<const This*>(this);
	return opposite(r.compare_intersecting( *new_this ));
      }
      else if ( upper_bound() == r.upper_bound() ) {
	return compare_same_interval(r);
      }
    }

    // Case 2: the intervals have common right endpoint
    if ( upper_bound() == r.upper_bound() ) {
      if ( lower_bound() < r.lower_bound() ) {
	Sign s = sign_at( r, Interval::LOWER );
	if ( s == CGAL::ZERO ) {
	  // the root of this is equal to the lower bound of r
	  //	    set_interval( r, Interval::LOWER );
	  set_interval( r.lower_bound() );
	  return CGAL::SMALLER;
	}
	else {
	  if ( s == sign_upper() ) {
	    return CGAL::SMALLER;
	  }
	  else {
	    //	      set_lower( r, Interval::LOWER, s );
	    set_lower( r.lower_bound(), s );
	    return compare_same_interval(r);
	  }
	}
      }
      else if ( r.lower_bound() < lower_bound() ) {
	const This* new_this = static_cast<const This*>(this);
	return opposite(r.compare_intersecting( *new_this ));
      }
      else if ( lower_bound() == r.lower_bound() ) {
	return compare_same_interval(r);
      }
    }

    // Case 3: the upper bound of r is an interior point of the
    //         interval of this and the lower bound of r is before
    //         the lower bound of this
    if ( lower_bound() > r.lower_bound() &&
	 upper_bound() > r.upper_bound() ) {

      Sign s_this_at_r = r.sign_at( *this, Interval::LOWER );
      if ( s_this_at_r == CGAL::ZERO ) {
	//	  r.set_interval( *this, Interval::LOWER );
	r.set_interval( lower_bound() );
	return CGAL::LARGER;
      }

      Sign s_r_at_this = sign_at( r, Interval::UPPER );
      if ( s_r_at_this == CGAL::ZERO ) {
	//	  set_interval( r, Interval::UPPER );
	set_interval( r.upper_bound() );
	return CGAL::SMALLER;
      }

      if ( s_r_at_this != CGAL::ZERO && s_this_at_r != CGAL::ZERO ) {

	if ( s_this_at_r != r.sign_lower() ) {
	  //	    r.set_upper( *this, Interval::LOWER, s_this_at_r );
	  r.set_upper( lower_bound(), s_this_at_r );
	  return CGAL::LARGER;
	}

	if ( s_r_at_this != sign_upper() ) {
	  //	    set_lower( r, Interval::UPPER, s_r_at_this );
	  set_lower( r.upper_bound(), s_r_at_this );
	  return CGAL::LARGER;
	}

	//	  set_upper( r, Interval::UPPER, s_r_at_this );
	//	  r.set_lower( *this, Interval::LOWER, s_this_at_r );
	set_upper( r.upper_bound(), s_r_at_this );
	r.set_lower( lower_bound(), s_this_at_r );

	return compare_same_interval(r);
      }
    }

    // Case 4: the upper bound of this is an interior point of the
    //         interval of r and the lower bound of this is before
    //         the lower bound of r
    if ( lower_bound() < r.lower_bound() &&
	 upper_bound() < r.upper_bound() ) {
      const This* new_this = static_cast<const This*>(this);
      return opposite(r.compare_intersecting( *new_this ));
    }

    // Case 5: the interval of r is contained in the interval of
    //         this
    if ( lower_bound() < r.lower_bound() &&
	 upper_bound() > r.upper_bound() ) {

      Sign s_rl_at_this = sign_at( r, Interval::LOWER );
      if ( s_rl_at_this == CGAL::ZERO ) {
	//	  set_interval( r, Interval::LOWER );
	set_interval( r.lower_bound() );
	return CGAL::SMALLER;
      }

      Sign s_ru_at_this = sign_at( r, Interval::UPPER );
      if ( s_ru_at_this == CGAL::ZERO ) {
	//	  set_interval( r, Interval::UPPER );
	set_interval( r.upper_bound() );
	return CGAL::LARGER;
      }

      Sign s_l = sign_lower();
      if ( s_l != s_rl_at_this ) {
	//	  set_upper( r, Interval::LOWER, s_rl_at_this );
	set_upper( r.lower_bound(), s_rl_at_this );
	return CGAL::SMALLER;
      }

      Sign s_u = sign_upper();
      if ( s_u != s_ru_at_this ) {
	//	  set_lower( r, Interval::UPPER, s_ru_at_this );
	set_lower( r.upper_bound(), s_ru_at_this );
	return CGAL::LARGER;
      }

      if ( s_rl_at_this != CGAL::ZERO &&
	   s_ru_at_this != CGAL::ZERO   ) {
	//	  set_lower( r, Interval::LOWER, s_rl_at_this );
	//	  set_upper( r, Interval::UPPER, s_ru_at_this );
	set_lower( r.lower_bound(), s_rl_at_this );
	set_upper( r.upper_bound(), s_ru_at_this );

	return compare_same_interval(r);
      }
    }

    // Case 6: the interval of this is contained in
    //         the interval of r
    if ( lower_bound() > r.lower_bound() &&
	 upper_bound() < r.upper_bound() ) {
      const This* new_this = static_cast<const This*>(this);
      return opposite(r.compare_intersecting( *new_this ));
    }

    CGAL_postcondition_msg(0, "This line should not have been reached.\n");
    //bool this_line_should_not_have_been_reached = false;