Classes | Functions
Gaudi::Math Namespace Reference

Classes

class  Lomont
 The equality comparison of double numbers using as the metric the maximal number of Units in the Last Place (ULP). More...
 
class  Lomont< double >
 the specialization for double numbers More...
 
class  Lomont< float >
 the specialization for float 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). More...
 
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). More...
 
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". More...
 
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". More...
 

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
afthe first number
bfthe second number
maxULPsthe maximal metric deviation in the terms of maximal number of units in the last place
Author
Vanya BELYAEV Ivan..nosp@m.Bely.nosp@m.aev@n.nosp@m.ikhe.nosp@m.f.nl
Date
2008-11-08

the final check

Definition at line 289 of file Lomont.cpp.

292 {
293  // ==========================================================================
295  static_assert( std::numeric_limits<double> ::is_specialized &&
296  std::numeric_limits<Long> ::is_specialized &&
297  std::numeric_limits<ULong> ::is_specialized &&
298  boost::integer_traits<ULong> ::is_specialized &&
299  boost::integer_traits<Long> ::is_specialized &&
300  sizeof(double)==sizeof(Long) &&
301  sizeof(double)==sizeof(ULong) &&
302  64 == std::numeric_limits<ULong>::digits , "FAILED ASSUMPTIONS") ;
303  // ==========================================================================
304  Cast_D caster ;
305 
306  //Long ai = *reinterpret_cast<const Long*>( &af ) ;
307  //Long bi = *reinterpret_cast<const Long*>( &bf ) ;
308  Long ai = caster.d2l ( af ) ;
309  Long bi = caster.d2l ( bf ) ;
310 
311  Long test = (((ULong)(ai^bi))>>63)-1;
312 
313  // assert ( (0==test) || ( boost::integer_traits<ULong>::const_max == test ) ) ;
314 
315  Long diff = ((( boost::integer_traits<Long>::const_min - ai ) & (~test)) | ( ai& test )) - bi ;
316 
317  Long maxDiff_ = maxULPs ;
318 
319  Long v1 = maxDiff_ + diff ;
320  Long v2 = maxDiff_ - diff ;
321 
322  return 0<=(v1|v2) ;
323 }
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
afthe first number
bfthe second number
maxULPsthe maximal metric deviation in the terms of maximal number of units in the last place
Author
Vanya BELYAEV Ivan..nosp@m.Bely.nosp@m.aev@n.nosp@m.ikhe.nosp@m.f.nl
Date
2008-11-08

Definition at line 193 of file Lomont.cpp.

196 {
197  // ==========================================================================
198  // prerequisites:
199  static_assert( std::numeric_limits<float> ::is_specialized &&
200  boost::integer_traits<int> ::is_specialized &&
201  boost::integer_traits<unsigned int> ::is_specialized &&
202  sizeof(float)==sizeof(int) &&
203  sizeof(float)==sizeof(unsigned int) &&
204  32 == boost::integer_traits<unsigned int>::digits , "FAILED ASSUMPTIONS") ;
205  // ==========================================================================
206 
207  Cast_F caster ;
208 
209  //int ai = *reinterpret_cast<const int*>( &af ) ;
210  //int bi = *reinterpret_cast<const int*>( &bf ) ;
211  int ai = caster.f2i ( af ) ;
212  int bi = caster.f2i ( bf ) ;
213 
214  int test = (((unsigned int)(ai^bi))>>31)-1;
215 
216  // assert ( (0==test) || ( boost::integer_traits<unsigned int>::const_max == test ) ) ;
217 
218  int diff = ((( boost::integer_traits<int>::const_min - ai ) & (~test)) | ( ai& test )) - bi ;
219 
220  int maxDiff_ = maxULPs ;
221 
222  int v1 = maxDiff_ + diff ;
223  int v2 = maxDiff_ - diff ;
224 
225  return 0<=(v1|v2) ;
226 }
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
afthe reference number
ulpsthe bias
Returns
the biased float number (on distance "ulps")
Author
Vanya BELYAEV Ivan..nosp@m.Bely.nosp@m.aev@n.nosp@m.ikhe.nosp@m.f.nl
Date
2008-11-08

the final check

Definition at line 342 of file Lomont.cpp.

343 {
344  // ==========================================================================
346  static_assert( std::numeric_limits<double> ::is_specialized &&
347  std::numeric_limits<Long> ::is_specialized &&
348  sizeof(double)==sizeof(Long) , "FAILED ASSUMPTIONS") ;
349  // ==========================================================================
350  Cast_D caster ;
351 
352  //Long al = *reinterpret_cast<const Long*>( &ad ) ;
353  Long al = caster.d2l ( ad ) ;
354  al += ulps ;
355  //return *reinterpret_cast<double*>(&al) ;
356  return caster.l2d ( al ) ;
357 }
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
afthe reference number
ulpsthe bias
Returns
the biased float number (on distance "ulps")
Author
Vanya BELYAEV Ivan..nosp@m.Bely.nosp@m.aev@n.nosp@m.ikhe.nosp@m.f.nl
Date
2008-11-08

the final check

Definition at line 245 of file Lomont.cpp.

246 {
248  static_assert( std::numeric_limits<float> ::is_specialized &&
249  std::numeric_limits<int> ::is_specialized &&
250  sizeof(float)==sizeof(int) , "FAILED ASSUMPTIONS") ;
251  // ==========================================================================
252  Cast_F caster ;
253 
254  // int ai = *reinterpret_cast<const int*>( &af ) ;
255  int ai = caster.f2i ( af ) ;
256  ai += ulps ;
257  // return *reinterpret_cast<float*>(&ai) ;
258  return caster.i2f ( ai ) ;
259 }