Gaudi::Math Namespace Reference

Classes
class	Lomont< float >
	the specialization for float numbers More...
class	Lomont< double >
	the specialization for double numbers More...
Functions
GAUDI_API bool	lomont_compare_float (const float af, const float bf, const unsigned short maxULPs)
	equality comparison of float numbers using as the metric the maximal number of Units in the Last Place (ULP).
GAUDI_API bool	lomont_compare_double (const double af, const double bf, const unsigned int maxULPs)
	equality comparison of double numbers using as the metric the maximal number of Units in the Last Place (ULP).
GAUDI_API float	next_float (const float af, const short ulps)
	Get the floating number that representation is different with respect to the argument for the certain number of "Units in the Last Position".
GAUDI_API double	next_double (const double af, const short ulps)
	Get the floating number that representation is different with respect to the argument for the certain number of "Units in the Last Position".

---

Function Documentation

bool Gaudi::Math::lomont_compare_double	(	const double	af,
		const double	bf,
		const unsigned int	maxULPs
	)

equality comparison of double numbers using as the metric the maximal number of Units in the Last Place (ULP).

It is a slightly modified version of very efficient implementation of the initial Bruce Dawson's algorithm by Chris Lomont.

See also:

www.lomont.org

http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm

C.Lomont claims the algorithm is factor 2-10 more efficient with respect to classical Knuth's algorithm for comparison of the floating number using the relative precision.

The effective relative difference depends on the choice of maxULPS:

• For the case of maxULPs=1, (of course it is totally unphysical case!!!) the effective relative precision r = |a-b|/(|a|+|b|)is ~6e-16 for |a|,|b|>1.e-304, and then it quickly goes to ~1

  const double a = ... ;
  const double b = ... ;

  const bool equal = Gaudi::Math::lomont_compare_double ( a , b ) ;

Parameters:

	af	the first number
	bf	the second number
	maxULPs	the maximal metric deviation in the terms of maximal number of units in the last place

Author:

Vanya BELYAEV Ivan.Belyaev@nikhef.nl

Date:

2008-11-08

the final check

Definition at line 291 of file Lomont.cpp.

{
  // ==========================================================================
  BOOST_STATIC_ASSERT( std::numeric_limits<double>  ::is_specialized &&
                       std::numeric_limits<Long>    ::is_specialized && 
                       std::numeric_limits<ULong>   ::is_specialized &&
                       boost::integer_traits<ULong> ::is_specialized &&
                       boost::integer_traits<Long>  ::is_specialized &&
                       sizeof(double)==sizeof(Long)                  && 
                       sizeof(double)==sizeof(ULong)                 &&
                       64 == std::numeric_limits<ULong>::digits      ) ;
  // ==========================================================================
  Cast_D caster ;
  
  //Long ai = *reinterpret_cast<const Long*>( &af ) ;
  //Long bi = *reinterpret_cast<const Long*>( &bf ) ;
  Long ai = caster.d2l ( af ) ;
  Long bi = caster.d2l ( bf ) ;
  
  Long test = (((ULong)(ai^bi))>>63)-1;
  
  // assert ( (0==test) || ( boost::integer_traits<ULong>::const_max == test ) ) ;
  
  Long diff = ((( boost::integer_traits<Long>::const_min - ai ) & (~test)) | ( ai& test )) - bi ;
  
  Long maxDiff_ = maxULPs ;
  
  Long v1 = maxDiff_ + diff ;
  Long v2 = maxDiff_ - diff ;
  
  return 0<=(v1|v2) ;
}

bool Gaudi::Math::lomont_compare_float	(	const float	af,
		const float	bf,
		const unsigned short	maxULPs
	)

equality comparison of float numbers using as the metric the maximal number of Units in the Last Place (ULP).

It is a slightly modified version of very efficient implementation of the initial Bruce Dawson's algorithm by Chris Lomont.

See also:

www.lomont.org

http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm

C.Lomont claims the algorithm is factor 2-10 more efficient with respect to classical Knuth's algorithm for comparison of the floating number using the relative precision.

The effective relative difference depends on the choice of maxULPS:

• For the case of maxULPs=1, (of course it is totally unphysical case!!!) the effective relative precision r = |a-b|/(|a|+|b|)is between 3.5e-8 and 5.5e-8 for |a|,|b|>1.e-37, and then it quickly goes to ~1
• For the case of maxULPS=10 the effective relative precision is between 3e-8 and 6e-7 for |a|,|b|>1.e-37, and then it quickly goes to ~1
• For the case of maxULPS=100 the effective relative precision is around ~6e-6 for |a|,|b|>1.e-37, and then it quickly goes to ~1
• For the case of maxULPS=1000 the effective relative precision is around ~6e-5 for |a|,|b|>1.e-37, and then it quickly goes to ~1

  const float a = ... ;
  const float b = ... ;

  const bool equal = Gaudi::Math::lomont_compare_float ( a , b ) ;

Parameters:

	af	the first number
	bf	the second number
	maxULPs	the maximal metric deviation in the terms of maximal number of units in the last place

Author:

Vanya BELYAEV Ivan.Belyaev@nikhef.nl

Date:

2008-11-08

Definition at line 195 of file Lomont.cpp.

{
  // ==========================================================================
  // prerequisites:
  BOOST_STATIC_ASSERT( std::numeric_limits<float>          ::is_specialized &&
                       boost::integer_traits<int>          ::is_specialized && 
                       boost::integer_traits<unsigned int> ::is_specialized &&
                       sizeof(float)==sizeof(int)                           && 
                       sizeof(float)==sizeof(unsigned int)                  &&
                       32 == boost::integer_traits<unsigned int>::digits    ) ;
  // ==========================================================================
  
  Cast_F caster ;
  
  //int ai = *reinterpret_cast<const int*>( &af ) ;
  //int bi = *reinterpret_cast<const int*>( &bf ) ;
  int ai = caster.f2i ( af ) ;
  int bi = caster.f2i ( bf ) ;
  
  int test = (((unsigned int)(ai^bi))>>31)-1;
  
  // assert ( (0==test) || ( boost::integer_traits<unsigned int>::const_max == test ) ) ;
  
  int diff = ((( boost::integer_traits<int>::const_min - ai ) & (~test)) | ( ai& test )) - bi ;
  
  int maxDiff_ = maxULPs ;
  
  int v1 = maxDiff_ + diff ;
  int v2 = maxDiff_ - diff ;
  
  return 0<=(v1|v2) ;
}

double Gaudi::Math::next_double	(	const double	af,
		const short	ulps
	)

Get the floating number that representation is different with respect to the argument for the certain number of "Units in the Last Position".

For ulps=1, it is just next float number, for ulps=-1 is is the previous one.

This routine is very convenient to test the parameter maxULPS for the routine Gaudi::Math::lomont_compare_double

See also:

Gaudi::Math::lomont_compare_double

Parameters:

	af	the reference number
	ulps	the bias

Returns:

the biased float number (on distance "ulps")

Author:

Vanya BELYAEV Ivan.Belyaev@nikhef.nl

Date:

2008-11-08

the final check

Definition at line 344 of file Lomont.cpp.

{
  // ==========================================================================
  BOOST_STATIC_ASSERT( std::numeric_limits<double> ::is_specialized &&
                       std::numeric_limits<Long>   ::is_specialized && 
                       sizeof(double)==sizeof(Long) ) ;
  // ==========================================================================
  Cast_D caster ;
  
  //Long al = *reinterpret_cast<const Long*>( &ad ) ;
  Long al = caster.d2l ( ad ) ;
  al += ulps ;
  //return *reinterpret_cast<double*>(&al) ;
  return caster.l2d ( al ) ;
}

float Gaudi::Math::next_float	(	const float	af,
		const short	ulps
	)

Get the floating number that representation is different with respect to the argument for the certain number of "Units in the Last Position".

For ulps=1, it is just next float number, for ulps=-1 is is the previous one.

This routine is very convenient to test the parameter maxULPS for the routine Gaudi::Math::lomont_compare_float

See also:

Gaudi:Math::lomont_compare_float

Parameters:

	af	the reference number
	ulps	the bias

Returns:

the biased float number (on distance "ulps")

Author:

Vanya BELYAEV Ivan.Belyaev@nikhef.nl

Date:

2008-11-08

the final check

Definition at line 247 of file Lomont.cpp.

{
  BOOST_STATIC_ASSERT( std::numeric_limits<float> ::is_specialized &&
                       std::numeric_limits<int>   ::is_specialized && 
                       sizeof(float)==sizeof(int) ) ;
  // ==========================================================================
  Cast_F caster ;
  
  // int ai = *reinterpret_cast<const int*>( &af ) ;
  int ai = caster.f2i ( af ) ;
  ai += ulps ;
  // return *reinterpret_cast<float*>(&ai) ;
  return   caster.i2f ( ai ) ;
}