---

gul_float.h

Go to the documentation of this file.

#ifndef GUL_FLOAT_H
#define GUL_FLOAT_H

namespace gul {

/*---------------------------------------------------------------------------
  check if two floating point numbers are approximately equal. in this check
  the relative error is measured.
  (  remark: according to knuth
          |a-b| <= eps*|b|  =>  a  ~  b  ('approximately rel_equal')
     and  |a-b| >  eps*|b|  =>  a !~~ b  (not 'essentialy rel_equal', for eps < 1
                                          & normalized floating point numbers)
     this check mixes both relations, 'approximately rel_equal' and 
     'essentialy rel_equal', so it isn't 100% correct!)
-----------------------------------------------------------------------------*/
template< class T >
inline bool rel_equal( const T& a, const T& b, const T& eps )
{
  return rtr<T>::fabs(a-b) <= rtr<T>::fabs(b) * eps; 
}

/*--------------------------------------------------------------------------
  check if two vectors are approximately equal. in this check
  the relative error is measured.
  ( remark: the maximum norm {..} is used to estimate the euclidian norm |..|
    by the equation:  {v} <= |v| <= {v}*sqrt(dim(v)) <= eps
    if this check succeeds, it is still possible that |v| <= eps < {v}*sqrt(dim),
    so this check isn't 100% correct!)
---------------------------------------------------------------------------*/
template< class T >
inline bool rel_equal( const gul::hpoint<T>& a, const gul::hpoint<T>& b, 
                   const T& eps )
{
  T tol = eps / (T)2;

  if( !rel_equal(a.x, b.x, tol) ) return false;
  if( !rel_equal(a.y, b.y, tol) ) return false;
  if( !rel_equal(a.z, b.z, tol) ) return false;
  if( !rel_equal(a.w, b.w, tol) ) return false;

  return true;
}

template< class T >
inline bool rel_equal( const gul::hpoint2<T>& a, const gul::hpoint2<T>& b, 
                   const T& eps )
{
  T tol = eps / rtr<T>::root2_3();

  if( !rel_equal(a.x, b.x, tol) ) return false;
  if( !rel_equal(a.y, b.y, tol) ) return false;
  if( !rel_equal(a.w, b.w, tol) ) return false;

  return true;
}

template< class T >
inline bool rel_equal( const gul::hpoint1<T>& a, const gul::hpoint1<T>& b, 
                   const T& eps )
{
  T tol = eps / rtr<T>::root2_2();

  if( !rel_equal(a.x, b.x, tol) ) return false;
  if( !rel_equal(a.w, b.w, tol) ) return false;

  return true;
}

template< class T >
inline bool rel_equal( const gul::point<T>& a, const gul::point<T>& b, 
                   const T& eps )
{
  T tol = eps / rtr<T>::root2_3();

  if( !rel_equal(a.x, b.x, tol) ) return false;
  if( !rel_equal(a.y, b.y, tol) ) return false;
  if( !rel_equal(a.z, b.z, tol) ) return false;

  return true;
}

template< class T >
inline bool rel_equal( const gul::point2<T>& a, const gul::point2<T>& b, 
                   const T& eps )
{
  T tol = eps / rtr<T>::root2_2();

  if( !rel_equal(a.x, b.x, tol) ) return false;
  if( !rel_equal(a.y, b.y, tol) ) return false;

  return true;
}

template< class T >
inline bool rel_equal( const gul::point1<T>& a, const gul::point1<T>& b, 
                   const T& eps )
{
  T tol = eps;

  if( !rel_equal(a.x, b.x, tol) ) return false;

  return true;
}

/*-----------------------------------------------------------------------
  approximate normalization of a vector
------------------------------------------------------------------------*/
template< class T >
T max_norm( const point<T>& a ) 
{
  T m;

  m = rtr<T>::fabs(a.x);
  if( rtr<T>::fabs(a.y) > m ) m = rtr<T>::fabs(a.y);
  if( rtr<T>::fabs(a.z) > m ) m = rtr<T>::fabs(a.z); 
  return m;
}

template< class T >
T max_norm( const point2<T>& a ) 
{
  T m;

  m = rtr<T>::fabs(a.x);
  if( rtr<T>::fabs(a.y) > m ) m = rtr<T>::fabs(a.y); 
  return m;
}

/*-----------------------------------------------------------------------
  check if the cross-product of two vectors is approximately zero
------------------------------------------------------------------------*/
template< class T >
bool zero_sine( const point<T>& ina, const point<T>& inb, const T& eps )
{
  T la,lb;
  point<T> a,b;

  // normalize a,b (approximately)
  la = max_norm(ina);
  if( !la ) return true;
  a = ((T)1/la)*ina;

  lb = max_norm(inb);
  if( !lb ) return true;
  b = ((T)1/lb)*inb;

  if( rel_equal(a.y*b.z,b.y*a.z,eps) &&
      rel_equal(b.x*a.z,a.x*b.z,eps) &&
      rel_equal(a.x*b.y,b.x*a.y,eps) )
    return true;

  return false;
}

template< class T >
bool zero_sine( const point2<T>& ina, const point2<T>& inb, const T& eps )
{
  T la,lb;

  // normalize a,b (approximately)
  la = max_norm(a);
  if( !la ) return true;
  a = ((T)1/la)*ina;

  lb = max_norm(b);
  if( !lb ) return true;
  b = ((T)1/lb)*inb;

  if( rel_equal(a.x*b.y-b.x*a.y,eps) )
    return true;

  return false;
}

/*---------------------------------------------------------------------------
  check if two floating point numbers are approximately equal. in this check
  the absolute error is measured.
-----------------------------------------------------------------------------*/
template< class T >
inline bool abs_equal( const T& a, const T& b, const T& eps )
{
  return rtr<T>::fabs(a-b) <= eps; 
}

/*--------------------------------------------------------------------------
  check if two vectors are approximately equal. in this check
  the absolute error is measured.
  ( remark: the maximum norm {..} is used to estimate the euclidian norm |..|
    by the equation:  {v} <= |v| <= {v}*sqrt(dim(v)) <= eps
    if this check succeeds, it is still possible that |v| <= eps < {v}*sqrt(dim),
    so this check isn't 100% correct!)
---------------------------------------------------------------------------*/
template< class T >
inline bool abs_equal( const gul::hpoint<T>& a, const gul::hpoint<T>& b, 
                       const T& eps )
{
  T tol = eps / (T)2;

  if( !abs_equal(a.x, b.x, tol) ) return false;
  if( !abs_equal(a.y, b.y, tol) ) return false;
  if( !abs_equal(a.z, b.z, tol) ) return false;
  if( !abs_equal(a.w, b.w, tol) ) return false;

  return true;
}

template< class T >
inline bool abs_equal( const gul::hpoint2<T>& a, const gul::hpoint2<T>& b, 
                   const T& eps )
{
  T tol = eps / rtr<T>::root2_3();

  if( !abs_equal(a.x, b.x, tol) ) return false;
  if( !abs_equal(a.y, b.y, tol) ) return false;
  if( !abs_equal(a.w, b.w, tol) ) return false;

  return true;
}

template< class T >
inline bool abs_equal( const gul::hpoint1<T>& a, const gul::hpoint1<T>& b, 
                   const T& eps )
{
  T tol = eps / rtr<T>::root2_2();

  if( !abs_equal(a.x, b.x, tol) ) return false;
  if( !abs_equal(a.w, b.w, tol) ) return false;

  return true;
}

template< class T >
inline bool abs_equal( const gul::point<T>& a, const gul::point<T>& b, 
                   const T& eps )
{
  T tol = eps / rtr<T>::root2_3();

  if( !abs_equal(a.x, b.x, tol) ) return false;
  if( !abs_equal(a.y, b.y, tol) ) return false;
  if( !abs_equal(a.z, b.z, tol) ) return false;

  return true;
}

template< class T >
inline bool abs_equal( const gul::point2<T>& a, const gul::point2<T>& b, 
                   const T& eps )
{
  T tol = eps / rtr<T>::root2_2();

  if( !abs_equal(a.x, b.x, tol) ) return false;
  if( !abs_equal(a.y, b.y, tol) ) return false;

  return true;
}

template< class T >
inline bool abs_equal( const gul::point1<T>& a, const gul::point1<T>& b, 
                   const T& eps )
{
  T tol = eps;

  if( !abs_equal(a.x, b.x, tol) ) return false;

  return true;
}


}

#endif


#if 0
##################################################################################
  OLD STUFF
##################################################################################


/*
template< class T >
inline bool equal( const T& a, const T& b, const T& eps )
{
  T scaled = rtr<T>::fabs(a-b)/ 
             rtr<T>::scalbn( (T)1,
               Max( rtr<T>::ilogb(rtr<T>::fabs(a)),
                    rtr<T>::ilogb(rtr<T>::fabs(b))) );
  return  scaled <= eps;
}

template< class T >
inline bool equal( const T& a, const T& b, const T& eps, T& scaled )
{
  scaled = rtr<T>::fabs(a-b)/ 
           rtr<T>::scalbn( (T)1,
             Max( rtr<T>::ilogb(rtr<T>::fabs(a)),
                  rtr<T>::ilogb(rtr<T>::fabs(b))) );
  return  scaled <= eps;
}
*/

/*
bool strongly_equal( const T& a, const T& b, const T& eps )
{
  T diff = rtr<T>::fabs(a-b);
  return (diff/rtr<T>::fabs(a) < eps) && (diff/rtr<T>::fabs(b) < eps);
}

bool greater( const T& a, const T& b, const T& eps )
{
  T diff = a-b;
  return (diff/rtr<T>::fabs(a) > eps) && (diff/rtr<T>::fabs(b) > eps);
}

bool lower( const T& a, const T& b, const T& eps )
{
  T diff = b-a;
  return (diff/rtr<T>::fabs(a) > eps) && (diff/rtr<T>::fabs(b) > eps);
}
*/

##################################################################

                                 /* -- check, if vector is (almost) zero --- */
/*
template< class T >
inline bool is_zero( const point<T>& a, const rtl<T>& t )
{
  T x,y,z;

  x = rtr<T>::fabs(a.x);
  if( x > t.zero_tol ) return false;
  y = rtr<T>::fabs(a.y);
  if( y > t.zero_tol ) return false;    
  z = rtr<T>::fabs(a.z);
  if( z > t.zero_tol ) return false;
    
  if( (x <= t.zero_tol_v3) && (y <= t.zero_tol_v3) && 
      (z <= t.zero_tol_v3) )
    return true;

  if( rtr<T>::sqrt(x*x + y*y + z*z) > t.zero_tol ) return false;

  return true;
}
template< class T >
inline bool is_zero( const point2<T>& a, const rtl<T>& t )
{
  T x,y;

  x = rtr<T>::fabs(a.x);
  if( x > t.zero_tol ) return false;
  y = rtr<T>::fabs(a.y);
  if( y > t.zero_tol ) return false;    
  
  if( (x <= t.zero_tol_v2) && (y <= t.zero_tol_v2) ) 
    return true;

  if( rtr<T>::sqrt(x*x + y*y) > t.zero_tol ) return false;

  return true;
}
*/

     /* check, if the difference of two homogeneous vectors is (almost) zero */
/*
template< class T >
inline bool dist_is_zero( const hpoint<T>& a, const hpoint<T>& b,
                          const rtl<T>& t )
{
  T x,y,z,w;
  
  if( !equal(a.x, b.x, t.zero_tol, x) ) return false;
  if( !equal(a.y, b.y, t.zero_tol, y) ) return false;
  if( !equal(a.z, b.z, t.zero_tol, z) ) return false;
  if( !equal(a.w, b.w, t.zero_tol, w) ) return false;
*/
  /*
  x = rtr<T>::fabs(a.x-b.x);
  if( x > t.zero_tol ) return false;
  y = rtr<T>::fabs(a.y-b.y);
  if( y > t.zero_tol ) return false;    
  z = rtr<T>::fabs(a.z-b.z);
  if( z > t.zero_tol ) return false;
  w = rtr<T>::fabs(a.w-b.w);
  if( w > t.zero_tol ) return false;

  if( (x <= t.zero_tol_v4) && (y <= t.zero_tol_v4) &&
      (z <= t.zero_tol_v4) && (w <= t.zero_tol_v4) )
    return true;

  if( rtr<T>::sqrt( x*x + y*y + z*z + w*w ) > t.zero_tol ) return false;
  */
/*
  return true;      
}
*/
   /* check, if the difference of two homogeneous vectors has a greater length
      than tol. */
/*
template< class T >
inline bool dist_longer_than( const hpoint<T>& a, const hpoint<T>& b, 
                              const T tolerance )
{
  T x,y,z,w;

  if( !equal(a.x, b.x, tolerance, x) ) return true;
  if( !equal(a.y, b.y, tolerance, y) ) return true;
  if( !equal(a.z, b.z, tolerance, z) ) return true;
  if( !equal(a.w, b.w, tolerance, w) ) return true;
*/
  /*
  x = rtr<T>::fabs(a.x-b.x);
  if( x > tolerance ) return true;
  y = rtr<T>::fabs(a.y-b.y);
  if( y > tolerance ) return true;    
  z = rtr<T>::fabs(a.z-b.z);
  if( z > tolerance ) return true;
  w = rtr<T>::fabs(a.w-b.w);
  if( w > tolerance ) return true;

  tol_v4 = tolerance / (T)2;  // sqrt(4)=2

  if( (x <= tol_v4) && (y <= tol_v4) && 
      (z <= tol_v4) && (w <= tol_v4) ) 
    return false;
  
  if( rtr<T>::sqrt(x*x + y*y + z*z + w*w) > tolerance ) return true;
  */
/*
  return false;  
}
*/

/*
template< class T >
inline bool dist_longer_than( const hpoint2<T>& a, const hpoint2<T>& b, 
                              const T tolerance )
{
  T x,y,w;

  if( !equal(a.x, b.x, tolerance, x) ) return true;
  if( !equal(a.y, b.y, tolerance, y) ) return true;
  if( !equal(a.w, b.w, tolerance, w) ) return true;
*/
  /*
  x = rtr<T>::fabs(a.x-b.x);
  if( x > tolerance ) return true;
  y = rtr<T>::fabs(a.y-b.y);
  if( y > tolerance ) return true;    
  w = rtr<T>::fabs(a.w-b.w);
  if( w > tolerance ) return true;

  tol_v3 = tolerance / rtr<T>::root2_3();

  if( (x <= tol_v3) && (y <= tol_v3) && (w <= tol_v3) ) 
    return false;
  
  if( rtr<T>::sqrt(x*x + y*y + w*w) > tolerance ) return true;
  */
/*
  return false;  
}
*/
                  /* check, if a vector has a greater length than tolerance */
/*
template< class T >
inline bool longer_than( const point<T>& a, const T tolerance )
{
  T x,y,z,tol_v3;

  x = rtr<T>::fabs(a.x);
  if( x > tolerance ) return true;
  y = rtr<T>::fabs(a.y);
  if( y > tolerance ) return true;    
  z = rtr<T>::fabs(a.z);
  if( z > tolerance ) return true;

  tol_v3 = tolerance / rtr<T>::root2_3();

  if( (x <= tol_v3) && (y <= tol_v3) && (z <= tol_v3) ) return false;
  
  if( rtr<T>::sqrt(x*x + y*y + z*z) > tolerance ) return true;

  return false;  
}
template< class T >
inline bool longer_than( const point2<T>& a, const T tolerance )
{
  T x,y,tol_v2;

  x = rtr<T>::fabs(a.x);
  if( x > tolerance ) return true;
  y = rtr<T>::fabs(a.y);
  if( y > tolerance ) return true;    

  tol_v2 = tolerance / rtr<T>::root2_2();

  if( (x <= tol_v2) && (y <= tol_v2) ) return false;
  
  if( rtr<T>::sqrt(x*x + y*y) > tolerance ) return true;

  return false;  
}
*/
   /* check, if the difference of two vectors has a greater length than tol. */
/*
template< class T >
inline bool dist_longer_than( const point<T>& a, const point<T>& b, 
                              const T tolerance )
{
  T x,y,z;

  if( !equal(a.x, b.x, tolerance, x) ) return true;
  if( !equal(a.y, b.y, tolerance, y) ) return true;
  if( !equal(a.z, b.z, tolerance, z) ) return true;
*/
  /*
  x = rtr<T>::fabs(a.x-b.x);
  if( x > tolerance ) return true;
  y = rtr<T>::fabs(a.y-b.y);
  if( y > tolerance ) return true;    
  z = rtr<T>::fabs(a.z-b.z);
  if( z > tolerance ) return true;

  tol_v3 = tolerance / rtr<T>::root2_3();

  if( (x <= tol_v3) && (y <= tol_v3) && (z <= tol_v3) ) return false;
  
  if( rtr<T>::sqrt(x*x + y*y + z*z) > tolerance ) return true;
  */
/*
  return false;  
}
*/

/*
template< class T >
inline bool dist_longer_than( const point2<T>& a, const point2<T>& b, 
                              const T tolerance )
{
  T x,y;

  if( !equal(a.x, b.x, tolerance, x) ) return true;
  if( !equal(a.y, b.y, tolerance, y) ) return true;
*/
  /*
  x = rtr<T>::fabs(a.x-b.x);
  if( x > tolerance ) return true;
  y = rtr<T>::fabs(a.y-b.y);
  if( y > tolerance ) return true;    

  tol_v2 = tolerance / rtr<T>::root2_2();

  if( (x <= tol_v2) && (y <= tol_v2) ) return false;
  
  if( rtr<T>::sqrt(x*x + y*y) > tolerance ) return true;
  */
/*
  return false;  
}
*/

################################################################################
#endif

---