compare_quadratic.h

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

  if ( CGAL::is_positive(Jp) )  { return SMALLER; }

  FT P4 = value_of_P4(J, Jp, G);

  Sign s_P4 = CGAL::sign(P4);
  if ( s_P4 == POSITIVE )  { return LARGER; }
  if ( s_P4 == NEGATIVE )  { return SMALLER; }
  return EQUAL;
}

template < class FT >
inline
Comparison_result
ke_compare_l1_r2(const FT& a1, const FT& b1, const FT& c1,
		 const FT& a2, const FT& b2, const FT& c2)
{
  FT J = value_of_J(a1, b1, a2, b2);
  Sign s_J = CGAL::sign(J);

#if 0
  if ( s_J == ZERO ) {
    FT D1 = value_of_D(a1, b1, c1);
    if ( CGAL::is_positive(D1) ) { return SMALLER; }

    FT D2 = value_of_D(a2, b2, c2);
    if ( CGAL::is_positive(D2) ) { return SMALLER; }

    return EQUAL;
  }
#endif

  if ( s_J == POSITIVE ) { return SMALLER; }

  FT a1c2 = a1 * c2;
  FT a2c1 = a2 * c1;
  FT K = value_of_K(a1c2, a2c1, b1 * b2);
  Sign s_K = CGAL::sign(K);

  if ( s_K == NEGATIVE ) { return SMALLER; }
#if 0
  if ( s_K == ZERO ) {
    FT D1 = value_of_D(a1, b1, c1);
    if ( CGAL::is_zero(D1) ) { return EQUAL; }

    FT D2 = value_of_D(a2, b2, c2);
    if ( CGAL::is_zero(D2) ) { return EQUAL; }

    return SMALLER;
  }
#endif

  FT Jp = value_of_Jp(b1, c1, b2, c2);
    
  if ( CGAL::is_positive(Jp) ) { return LARGER; }

  FT P4 = value_of_P4(J, Jp, a1c2 - a2c1);
    
  Sign s_P4 = CGAL::sign(P4);
  if ( s_P4 == POSITIVE ) { return SMALLER; }
  if ( s_P4 == NEGATIVE ) { return LARGER; }
  return EQUAL;
}

template < class FT >
inline
Comparison_result
ke_compare_r1_l2(const FT& a1, const FT& b1, const FT& c1,
		 const FT& a2, const FT& b2, const FT& c2)
{
  FT J = value_of_J(a1, b1, a2, b2);
  Sign s_J = CGAL::sign(J);

#if 0
  if ( s_J == ZERO ) {
    FT D1 = value_of_D(a1, b1, c1);
    if ( CGAL::is_positive(D1) ) { return LARGER; }

    FT D2 = value_of_D(a2, b2, c2);
    if ( CGAL::is_positive(D2) ) { return LARGER; }

    return EQUAL;
  }
#endif

  if ( s_J == NEGATIVE ) { return LARGER; }

  FT a1c2 = a1 * c2;
  FT a2c1 = a2 * c1;
  FT K = value_of_K(a1c2, a2c1, b1 * b2);
  Sign s_K = CGAL::sign(K);

  if ( s_K == NEGATIVE ) { return LARGER; }

#if 0
  if ( s_K == ZERO ) {
    FT D1 = value_of_D(a1, b1, c1);
    if ( CGAL::is_zero(D1) ) { return EQUAL; }

    FT D2 = value_of_D(a2, b2, c2);
    if ( CGAL::is_zero(D2) ) { return EQUAL; }

    return LARGER;
  }
#endif

  FT Jp = value_of_Jp(b1, c1, b2, c2);
    
  if ( CGAL::is_negative(Jp) ) { return SMALLER; }
    
  FT P4 = value_of_P4(J, Jp, a1c2 - a2c1);

  Sign s_P4 = CGAL::sign(P4);
  if ( s_P4 == POSITIVE ) { return LARGER; }
  if ( s_P4 == NEGATIVE ) { return SMALLER; }
  return EQUAL;
}

template < class FT >
inline
Comparison_result
ke_compare_r1_r2(const FT& a1, const FT& b1, const FT& c1,
		 const FT& a2, const FT& b2, const FT& c2)
{
#ifdef COMPARATOR_PROFILER
  comparator_profiler::counter_rr++;
#endif

  FT J = value_of_J(a1, b1, a2, b2);
  Sign s_J = CGAL::sign(J);

  FT a1c2 = a1 * c2;
  FT a2c1 = a2 * c1;
  FT K = value_of_K(a1c2, a2c1, b1 * b2);
  Sign s_K = CGAL::sign(K);
    
  if ( s_J == POSITIVE ) {
    if ( s_K == POSITIVE )  { return SMALLER; }
    else if ( s_K == NEGATIVE ) {

#ifdef COMPARATOR_PROFILER
      comparator_profiler::counter_rr_p3inf++;
#endif

      FT G = a1c2 - a2c1;
      FT P3inf = value_of_P3inf(a1, b1, J, G);

      if ( !(CGAL::is_negative(P3inf)) ) { return SMALLER; }

#ifdef COMPARATOR_PROFILER
      comparator_profiler::counter_rr_p4++;
#endif

      FT Jp = value_of_Jp(b1, c1, b2, c2);

      if ( CGAL::is_negative(Jp) )  { return LARGER; }

      FT P4 = value_of_P4(J, Jp, G);

      Sign s_P4 = CGAL::sign(P4);
      if ( s_P4 == POSITIVE )  { return SMALLER; }
      if ( s_P4 == NEGATIVE )  { return LARGER; }
      return EQUAL;
    } else { // K = 0
      FT D2 = value_of_D(a2, b2, c2);
      if ( CGAL::is_zero(D2) ) { return EQUAL; }

      return SMALLER;
    }
  } else if ( s_J == NEGATIVE ) { // J < 0
    if ( s_K == POSITIVE )  { return LARGER; }
    else if ( s_K == NEGATIVE ) {

#ifdef COMPARATOR_PROFILER
      comparator_profiler::counter_rr_p3inf++;
#endif

      FT G = a1c2 - a2c1;
      FT P3inf = value_of_P3inf(a1, b1, J, G);

      if ( !(CGAL::is_positive(P3inf)) ) { return LARGER; }

#ifdef COMPARATOR_PROFILER
      comparator_profiler::counter_rr_p4++;
#endif

      FT Jp = value_of_Jp(b1, c1, b2, c2);

      if ( CGAL::is_positive(Jp) )   { return SMALLER; }

      FT P4 = value_of_P4(J, Jp, G);

      Sign s_P4 = CGAL::sign(P4);
      if ( s_P4 == POSITIVE )  { return LARGER; }
      if ( s_P4 == NEGATIVE )  { return SMALLER; }
      return EQUAL;
    } else { // K = 0
      FT D1 = value_of_D(a1, b1, c1);
      if ( CGAL::is_zero(D1) ) { return EQUAL; }

      return LARGER;
    }
  }

  // J = 0
  Sign s_G = CGAL::sign( value_of_G(a1, c1, a2, c2) );
  if ( s_G == NEGATIVE ) { return SMALLER; }
  if ( s_G == POSITIVE ) { return LARGER; }
  return EQUAL;
}


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

//--------------------------------------------------------------------
//--------------------------------------------------------------------
// the following functions do the filtering for the r1-r2 tree without 
// using C++ exceptions

template < class CT, class ET >
inline
Comparison_result
sqrt_compare_r1_r2_filtered(const CT& a1, const CT& b1, const CT& c1,
			    const CT& a2, const CT& b2, const CT& c2)
{
  typedef Interval_nt<false> IT;

  FPU_CW_t backup = FPU_get_cw();
  FPU_set_cw(CGAL_FE_UPWARD);

  IT a1i(a1), b1i(b1), c1i(c1);
  IT a2i(a2), b2i(b2), c2i(c2);

  IT D1 = value_of_D(a1i, b1i, c1i);
  IT D2 = value_of_D(a2i, b2i, c2i);

  IT r1 = ( b1i + CGAL::sqrt(D1) ) / a1i;
  IT r2 = ( b2i + CGAL::sqrt(D2) ) / a2i;

  FPU_set_cw(backup);

  if ( r1.sup() < r2.inf() ) { return SMALLER; }
  if ( r1.inf() > r2.sup() ) { return LARGER; }

  return sqrt_compare_r1_r2(ET(a1), ET(b1), ET(c1),
			    ET(a2), ET(b2), ET(c2));
}

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

template < class CT, class ET >
inline
Comparison_result
dfmt_compare_r1_r2_filtered(const CT& a1, const CT& b1, const CT& c1,
			    const CT& a2, const CT& b2, const CT& c2)
{
  typedef Interval_nt<false> IT;

  FPU_CW_t backup = FPU_get_cw();
  FPU_set_cw(CGAL_FE_UPWARD);

  IT a1i(a1), b1i(b1), c1i(c1);
  IT a2i(a2), b2i(b2), c2i(c2);

  IT J = value_of_J(a1i, b1i, a2i, b2i);
  IT K = value_of_K(a1i, b1i, c1i, a2i, b2i, c2i);

  if ( J.inf() > 0 ){
    if ( K.inf() > 0 ) {
      FPU_set_cw(backup);
      return SMALLER;
    }
    if ( K.sup() > 0 ) {
      FPU_set_cw(backup);
      return dfmt_compare_r1_r2(ET(a1), ET(b1), ET(c1),
				ET(a2), ET(b2), ET(c2));
    }

    IT D1 = value_of_D(a1i, b1i, c1i);
    IT D2 = value_of_D(a2i, b2i, c2i);

    IT D = value_of_D(a1i, D1, a2i, D2);

    if ( D.sup() < 0 ) {
      FPU_set_cw(backup);
      return SMALLER;
    }
    if ( D.inf() < 0 ) {
      FPU_set_cw(backup);
      return dfmt_compare_r1_r2(ET(a1), ET(b1), ET(c1),
				ET(a2), ET(b2), ET(c2));
    }

    IT Jp = value_of_Jp(b1i, c1i, b2i, c2i);

    if ( Jp.sup() < 0 ) {
      FPU_set_cw(backup);
      return LARGER;
    }
    if ( Jp.inf() < 0 ) {
      FPU_set_cw(backup);
      return dfmt_compare_r1_r2(ET(a1), ET(b1), ET(c1),
				ET(a2), ET(b2), ET(c2));
    }
    

    IT R0 = value_of_R0(D1, D2, K);
    FPU_set_cw(backup);

    if ( R0.sup() < 0 ) { return SMALLER; }
    if ( R0.inf() > 0 ) { return LARGER; }
    return dfmt_compare_r1_r2(ET(a1), ET(b1), ET(c1),
			      ET(a2), ET(b2), ET(c2));
  } else if ( J.sup() > 0 ) {
    FPU_set_cw(backup);
    return dfmt_compare_r1_r2(ET(a1), ET(b1), ET(c1),
			      ET(a2), ET(b2), ET(c2));
  } else { // J < 0
    if ( K.inf() > 0 ) {
      FPU_set_cw(backup);
      return LARGER;
    }
    if ( K.sup() > 0 ) {
      FPU_set_cw(backup);
      return dfmt_compare_r1_r2(ET(a1), ET(b1), ET(c1),
				ET(a2), ET(b2), ET(c2));
    }

    IT D1 = value_of_D(a1i, b1i, c1i);
    IT D2 = value_of_D(a2i, b2i, c2i);

    IT D = value_of_D(a1i, D1, a2i, D2);

    if ( D.inf() > 0 ) {
      FPU_set_cw(backup);
      return LARGER;
    }
    if ( D.sup() > 0 ) {
      FPU_set_cw(backup);
      return dfmt_compare_r1_r2(ET(a1), ET(b1), ET(c1),
				ET(a2), ET(b2), ET(c2));
    }


    IT Jp = value_of_Jp(b1i, c1i, b2i, c2i);

    if ( Jp.inf() > 0 ) {
       FPU_set_cw(backup);
      return SMALLER;
    }
    if ( Jp.sup() > 0 ) {
      FPU_set_cw(backup);
      return dfmt_compare_r1_r2(ET(a1), ET(b1), ET(c1),
				ET(a2), ET(b2), ET(c2));
    }

    IT R0 = value_of_P4(D1, D2, K);
    FPU_set_cw(backup);

    if ( R0.sup() < 0 ) { return LARGER; }
    if ( R0.inf() > 0 ) { return SMALLER; }

    return dfmt_compare_r1_r2(ET(a1), ET(b1), ET(c1),
			      ET(a2), ET(b2), ET(c2));
  }
}

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

template < class CT, class ET >
inline
Comparison_result
ke_compare_r1_r2_filtered(const CT& a1, const CT& b1, const CT& c1,
			  const CT& a2, const CT& b2, const CT& c2)
{
  typedef Interval_nt<false> IT;

  FPU_CW_t backup = FPU_get_cw();
  FPU_set_cw(CGAL_FE_UPWARD);

  IT a1i(a1), b1i(b1), c1i(c1);
  IT a2i(a2), b2i(b2), c2i(c2);

  IT J = value_of_J(a1i, b1i, a2i, b2i);

  IT a1c2 = a1i * c2i;
  IT a2c1 = a2i * c1i;

  IT K = value_of_K(a1c2, a2c1, b1i * b2i);
    
  if ( J.inf() > 0 ) {
    if ( K.inf() > 0 )  { 
      FPU_set_cw(backup);
      return SMALLER;
    } else if ( K.sup() < 0 ) {
      IT G = a1c2 - a2c1;
      IT P3inf = value_of_P3inf(a1i, b1i, J, G);

      if ( P3inf.inf() >= 0 ) {
	FPU_set_cw(backup);
	return SMALLER;
      }
      if ( P3inf.sup() > 0 ) {
	FPU_set_cw(backup);
	return ke_compare_r1_r2(ET(a1), ET(b1), ET(c1),
				ET(a2), ET(b2), ET(c2));
      }

      IT Jp = value_of_Jp(b1i, c1i, b2i, c2i);

      if ( Jp.sup() < 0 ) {
	FPU_set_cw(backup);
	return LARGER;
      }
      if ( Jp.inf() < 0 ) {
	FPU_set_cw(backup);
	return ke_compare_r1_r2(ET(a1), ET(b1), ET(c1),
				ET(a2), ET(b2), ET(c2));
      }

      IT P4 = value_of_P4(J, Jp, G);
      FPU_set_cw(backup);

      if ( P4.inf() > 0 ) { return SMALLER; }
      if ( P4.sup() < 0 ) { return LARGER; }

      return ke_compare_r1_r2(ET(a1), ET(b1), ET(c1),
			      ET(a2), ET(b2), ET(c2));
    } else {
      FPU_set_cw(backup);
      return ke_compare_r1_r2(ET(a1), ET(b1), ET(c1),
			      ET(a2), ET(b2), ET(c2));
    }
  } else if ( J.sup() < 0 ) { // J < 0
    if ( K.inf() > 0 ) {
      FPU_set_cw(backup);
      return LARGER;
    } else if ( K.sup() < 0 ) {
      IT G = a1c2 - a2c1;
      IT P3inf = value_of_P3inf(a1i, b1i, J, G);

      if ( P3inf.sup() <= 0 ) {
	FPU_set_cw(backup);
	return LARGER;
      }
      if ( P3inf.inf() < 0 ) {
	FPU_set_cw(backup);
	return ke_compare_r1_r2(ET(a1), ET(b1), ET(c1),
				ET(a2), ET(b2), ET(c2));
      }

      IT Jp = value_of_Jp(b1i, c1i, b2i, c2i);

      if ( Jp.inf() > 0 ) {
	FPU_set_cw(backup);
	return SMALLER;
      }
      if ( Jp.sup() > 0 ) {
	FPU_set_cw(backup);
	return ke_compare_r1_r2(ET(a1), ET(b1), ET(c1),
				ET(a2), ET(b2), ET(c2));
      }

      IT P4 = value_of_P4(J, Jp, G);
      FPU_set_cw(backup);

      if ( P4.inf() > 0 ) { return LARGER; }
      if ( P4.sup() < 0 ) { return SMALLER; }

      return ke_compare_r1_r2(ET(a1), ET(b1), ET(c1),
			      ET(a2), ET(b2), ET(c2));
    } else {
      FPU_set_cw(backup);
      return ke_compare_r1_r2(ET(a1), ET(b1), ET(c1),
			      ET(a2), ET(b2), ET(c2));
    }
  } else { // J = ?
    FPU_set_cw(backup);
    return ke_compare_r1_r2(ET(a1), ET(b1), ET(c1),
			    ET(a2), ET(b2), ET(c2));
  }
}





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

CGAL_APOLLONIUS_GRAPH_2_END_NAMESPACE

CGAL_END_NAMESPACE

#endif // CGAL_APOLLONIUS_GRAPH_2_COMPARE_QUADRATIC_H