compare_quadratic.h

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

// Copyright (c) 2003,2006  INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL 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/Apollonius_graph_2/include/CGAL/Apollonius_graph_2/compare_quadratic.h $
// $Id: compare_quadratic.h 32776 2006-07-30 13:02:57Z spion $
// 
//
// Author(s)     : Menelaos Karavelas <[e-mail]>



#ifndef CGAL_APOLLONIUS_GRAPH_2_COMPARE_QUADRATIC_H
#define CGAL_APOLLONIUS_GRAPH_2_COMPARE_QUADRATIC_H

#include <CGAL/Apollonius_graph_2/basic.h>
#include <CGAL/functions_on_signs.h>

CGAL_BEGIN_NAMESPACE

CGAL_APOLLONIUS_GRAPH_2_BEGIN_NAMESPACE

#ifdef COMPARATOR_PROFILER
#include <CGAL/Apollonius_graph_2/comparator_profiler.h>
#endif


//--------------------------------------------------------------------
// help functions for the compulation of various quantites

template < class FT >
inline FT
value_of_D(const FT& a, const FT& b, const FT& c)
{
  return CGAL::square(b) - a * c;
}

template < class FT >
inline FT
value_of_J(const FT& a1, const FT& b1, const FT& a2, const FT& b2)
{
  return (a1 * b2 - a2 * b1);
}

template < class FT >
inline FT
value_of_Jp(const FT& b1, const FT& c1, const FT& b2, const FT& c2)
{
  return (b1 * c2 - b2 * c1);
}

template < class FT >
inline FT
value_of_G(const FT& a1, const FT& c1, const FT& a2, const FT& c2)
{
  return (a1 * c2 - a2 * c1);
}

template < class FT >
inline FT
value_of_K(const FT& a1c2, const FT& a2c1, const FT& b1b2)
{
  return a1c2 + a2c1 - FT(2) * b1b2;
}

template < class FT >
inline FT
value_of_K(const FT& a1, const FT& b1, const FT& c1,
	   const FT& a2, const FT& b2, const FT& c2)
{
  return c1 * a2 + a1 * c2 - FT(2) * b1 * b2;
}

template < class FT >
inline FT
value_of_R0(const FT&  J, const FT& Jp,
	    const FT& a1, const FT& c1,
	    const FT& a2, const FT& c2)
{
  return CGAL::square(a1 * c2 - c1 * a2) - FT(4) * J * Jp;
}

template < class FT >
inline FT
value_of_R0(const FT& D1, const FT& D2, const FT& K)
{
  return CGAL::square(K) - FT(4) * D1 * D2;
}

template < class FT >
inline FT
value_of_P4(const FT& D1, const FT& D2, const FT& K)
{
  return FT(4) * D1 * D2 - CGAL::square(K);
}

template < class FT >
inline FT
value_of_D(const FT& a1, const FT& D1, const FT& a2, const FT& D2)
{
  return D1 * CGAL::square(a2) - D2 * CGAL::square(a1);
}

template < class FT >
inline FT
value_of_P3inf(const FT& a1, const FT& b1, const FT& J, const FT& G)
{
  return FT(2) * b1 * J - a1 * G;
}

template < class FT >
inline FT
value_of_P3inf_bis(const FT& a1, const FT& K, const FT& a2, const FT& D1)
{
  return -(a1 * K + FT(2) * a2 * D1);
}

template < class FT >
inline FT
value_of_P3pinf(const FT& c1, const FT& J, const FT& a1, const FT& Jp)
{
  return c1 * J - a1 * Jp;
}

template < class FT >
inline FT
value_of_L(const FT& a1, const FT& c2, const FT& b1, const FT& b2)
{
  return (a1 * c2 - b1 * b2);
}

template < class FT >
inline FT
value_of_Lp(const FT& a2, const FT& c1, const FT& b1, const FT& b2)
{
  return (a2 * c1 - b1 * b2);
}

template < class FT >
inline FT
value_of_A(const FT& b1, const FT& J, const FT& a1, const FT& L)
{
  return b1 * J + a1 * L;
}

template < class FT >
inline FT
value_of_Q3(const FT& b2, const FT& a2, const FT& J, const FT& L)
{
  return (b2 * J + a2 * L);
}

template < class FT >
inline FT
value_of_Q3(const FT& a2, const FT& b2, const FT& J, const FT& G,
	    const FT& K)
{
  return FT(2) * b2 * J - a2 * (G - K);
} 

template < class FT >
inline FT
value_of_Q3p(const FT& a1, const FT& b1, const FT& J, const FT& G,
	     const FT& K)
{
  return a1 * (G + K) - FT(2) * b1 * J;
} 


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

// the trivial method that uses square roots

template < class FT >
inline
Comparison_result
sqrt_compare_l1_l2(const FT& a1, const FT& b1, const FT& c1,
		   const FT& a2, const FT& b2, const FT& c2)
{
  FT D1 = value_of_D(a1, b1, c1);
  FT D2 = value_of_D(a2, b2, c2);

  FT l1 = ( b1 - CGAL::sqrt(D1) ) / a1;
  FT l2 = ( b2 - CGAL::sqrt(D2) ) / a2;

  return CGAL::compare(l1, l2);
}

template < class FT >
inline
Comparison_result
sqrt_compare_l1_r2(const FT& a1, const FT& b1, const FT& c1,
		   const FT& a2, const FT& b2, const FT& c2)
{
  FT D1 = value_of_D(a1, b1, c1);
  FT D2 = value_of_D(a2, b2, c2);

  FT l1 = ( b1 - CGAL::sqrt(D1) ) / a1;
  FT r2 = ( b2 + CGAL::sqrt(D2) ) / a2;

  return CGAL::compare(l1, r2);
}


template < class FT >
inline
Comparison_result
sqrt_compare_r1_l2(const FT& a1, const FT& b1, const FT& c1,
		   const FT& a2, const FT& b2, const FT& c2)
{
  FT D1 = value_of_D(a1, b1, c1);
  FT D2 = value_of_D(a2, b2, c2);

  FT r1 = ( b1 + CGAL::sqrt(D1) ) / a1;
  FT l2 = ( b2 - CGAL::sqrt(D2) ) / a2;

  return CGAL::compare(r1, l2);
}

template < class FT >
inline
Comparison_result
sqrt_compare_r1_r2(const FT& a1, const FT& b1, const FT& c1,
		   const FT& a2, const FT& b2, const FT& c2)
{
  FT D1 = value_of_D(a1, b1, c1);
  FT D2 = value_of_D(a2, b2, c2);

  FT r1 = ( b1 + CGAL::sqrt(D1) ) / a1;
  FT r2 = ( b2 + CGAL::sqrt(D2) ) / a2;

  return CGAL::compare(r1, r2);
}



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


// the DFMT evaluation trees

template < class FT >
inline
Comparison_result
dfmt_compare_l1_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);
  FT K = value_of_K(a1, b1, c1, a2, b2, c2);
    
  if ( CGAL::is_positive(J) ) {
    if ( CGAL::is_positive(K) )  return SMALLER;  // l1 < l2
	
    FT D1 = value_of_D(a1, b1, c1);
    FT D2 = value_of_D(a2, b2, c2);

    FT D = value_of_D(a1, D1, a2, D2);

    if ( CGAL::is_positive(D) )  return SMALLER;  // l1 < l2
	
    FT Jp = value_of_Jp(b1, c1, b2, c2);

    if ( CGAL::is_negative(Jp) )  return LARGER;   // l1 > l2
	
    FT R0 = value_of_R0(D1, D2, K);
	
    Sign s_R0 = CGAL::sign(R0);
    if ( s_R0 == NEGATIVE )  return SMALLER;  // l1 < l2
    if ( s_R0 == POSITIVE )  return LARGER;   // l1 > l2
    return EQUAL;
  } else { // J<0
    if ( CGAL::is_positive(K) )  return LARGER;   // l1 > l2
	
    FT D1 = value_of_D(a1, b1, c1);
    FT D2 = value_of_D(a2, b2, c2);

    FT D = value_of_D(a1, D1, a2, D2);

    if ( CGAL::is_negative(D) )  return LARGER;   // l1 > l2

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

    if ( CGAL::is_positive(Jp) )  return SMALLER;  // l1 < l2

    FT R0 = value_of_R0(D1, D2, K);

    Sign s_R0 = CGAL::sign(R0);
    if ( s_R0 == NEGATIVE )  return LARGER;   // l1 > l2
    if ( s_R0 == POSITIVE )  return SMALLER;  // l1 < l2
    return EQUAL;
  }
}

template < class FT >
inline
Comparison_result
dfmt_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);

  if ( CGAL::is_positive(J) ) return SMALLER;   // l1 < r2

  FT K = value_of_K(a1, b1, c1, a2, b2, c2);

  if ( CGAL::is_negative(K) ) return SMALLER;   // l1 < r2

  FT Jp = value_of_Jp(b1, c1, b2, c2);
    
  if ( CGAL::is_positive(Jp) ) return LARGER;  // l1 > r2
    
  FT R0 = value_of_R0(J, Jp, a1, c1, a2, c2);
    
  Sign s_R0 = CGAL::sign(R0);
  if ( s_R0 == NEGATIVE ) return SMALLER;   // l1 < r2
  if ( s_R0 == POSITIVE ) return LARGER;    // l1 > r2
  return EQUAL;
}



template < class FT >
inline
Comparison_result
dfmt_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);

  if ( CGAL::is_negative(J) ) return LARGER;   // r1 > l2

  FT K = value_of_K(a1, b1, c1, a2, b2, c2);

  if ( CGAL::is_negative(K) ) return LARGER;   // r1 > l2

  FT Jp = value_of_Jp(b1, c1, b2, c2);
    
  if ( CGAL::is_negative(Jp) ) return SMALLER;  // r1 < l2
    
  FT R0 = value_of_R0(J, Jp, a1, c1, a2, c2);
    
  Sign s_R0 = CGAL::sign(R0);
  if ( s_R0 == NEGATIVE ) return LARGER;   // r1 > l2
  if ( s_R0 == POSITIVE ) return SMALLER;  // r1 < l2
  return EQUAL;
}

template < class FT >
inline
Comparison_result
dfmt_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);
  FT K = value_of_K(a1, b1, c1, a2, b2, c2);
    
  if ( CGAL::is_positive(J) ){
    if ( CGAL::is_positive(K) )  return SMALLER;   // r1 < r2   1,2

#ifdef COMPARATOR_PROFILER
    comparator_profiler::counter_rr_e++;
#endif
	
    FT D1 = value_of_D(a1, b1, c1);
    FT D2 = value_of_D(a2, b2, c2);

    FT D = value_of_D(a1, D1, a2, D2);

    if ( CGAL::is_negative(D) )  return SMALLER;   // r1 < r2   2,3b

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

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

    if ( CGAL::is_negative(Jp) )  return LARGER;    // r1 > r2   3a

    FT R0 = value_of_R0(D1, D2, K);

    Sign s_R0 = CGAL::sign(R0);
    if ( s_R0 == NEGATIVE )  return SMALLER;   // r1 < r2   2
    if ( s_R0 == POSITIVE )  return LARGER;    // r1 > r2   3a
    return EQUAL;
  } else { // J<0
    if ( CGAL::is_positive(K) )  return LARGER;   // r1 > r2    4,5

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

    FT D1 = value_of_D(a1, b1, c1);
    FT D2 = value_of_D(a2, b2, c2);

    FT D = value_of_D(a1, D1, a2, D2);
	
    if ( CGAL::is_positive(D) )  return LARGER;   // r1 > r2    3a,4
	
#ifdef COMPARATOR_PROFILER
    comparator_profiler::counter_rr_r0++;
#endif

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

    if ( CGAL::is_positive(Jp) )  return SMALLER;   // r1 < r2     3b

    FT R0 = value_of_P4(D1, D2, K);

    Sign s_R0 = CGAL::sign(R0);
    if ( s_R0 == NEGATIVE )  return LARGER;   // r1 > r2     4
    if ( s_R0 == POSITIVE )  return SMALLER;  // r1 < r2     3b
    return EQUAL;
  }
}

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

// the KE evaluation trees

template < class FT >
inline
Comparison_result
ke_compare_l1_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 ( s_J == ZERO ) {
    Sign s_G = CGAL::sign( value_of_G(a1, c1, a2, c2) );
    if ( s_G == POSITIVE ) { return SMALLER; }
    if ( s_G == NEGATIVE ) { return LARGER; }
    return EQUAL;
  }

  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; }

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

      return SMALLER;
    }

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

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

    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;
  }

  // J < 0
  if ( s_K == POSITIVE )  { return LARGER; }

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

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

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

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