|
Gaudi Framework, version v21r6 |
| Home | Generated: 11 Nov 2009 |
|
|
Gaudi Framework, version v21r6 |
| Home | Generated: 11 Nov 2009 |
00001 // $Id: Lomont.cpp,v 1.2 2009/05/23 17:08:16 ibelyaev Exp $ 00002 // ============================================================================ 00003 // Include files 00004 // ============================================================================ 00005 // STD & STL 00006 // ============================================================================ 00007 #include <iostream> 00008 #include <limits> 00009 #include <cassert> 00010 #include <cmath> 00011 #include <vector> 00012 // ============================================================================ 00013 // GaudiKernel 00014 // ============================================================================ 00015 #include "GaudiKernel/Lomont.h" 00016 // ============================================================================ 00017 // Boost 00018 // ============================================================================ 00019 #include "boost/static_assert.hpp" 00020 #include "boost/integer_traits.hpp" 00021 // ============================================================================ 00022 namespace 00023 { 00024 // prerequisites for "correct" Float 00025 BOOST_STATIC_ASSERT( std::numeric_limits<float> ::is_specialized && 00026 boost::integer_traits<int> ::is_specialized && 00027 boost::integer_traits<unsigned int> ::is_specialized && 00028 sizeof(float)==sizeof(int) && 00029 sizeof(float)==sizeof(unsigned int) && 00030 32 == boost::integer_traits<unsigned int>::digits ) ; 00031 // ========================================================================== 00032 // define proepr double 00033 // ========================================================================== 00034 #ifdef _WIN32 00035 struct _Longs 00036 { 00037 typedef __int64 Long ; 00038 typedef unsigned __int64 ULong ; 00039 } ; 00040 #else 00041 template <bool I> 00042 struct __Longs ; 00043 template <> 00044 struct __Longs<true> 00045 { 00046 typedef long Long ; 00047 typedef unsigned long ULong ; 00048 } ; 00049 template <> 00050 struct __Longs<false> 00051 { 00052 typedef long long Long ; 00053 typedef unsigned long long ULong ; 00054 } ; 00055 struct _Longs : public __Longs<sizeof(double)==sizeof(long)>{}; 00056 #endif 00057 00059 typedef _Longs::Long Long ; 00060 typedef _Longs::ULong ULong ; 00061 // ========================================================================== 00063 BOOST_STATIC_ASSERT( std::numeric_limits<double> ::is_specialized && 00064 std::numeric_limits<Long> ::is_specialized && 00065 std::numeric_limits<ULong> ::is_specialized && 00066 boost::integer_traits<ULong> ::is_specialized && 00067 boost::integer_traits<Long> ::is_specialized && 00068 sizeof(double)==sizeof(Long) && 00069 sizeof(double)==sizeof(ULong) && 00070 64 == std::numeric_limits<ULong>::digits ) ; 00071 // ========================================================================== 00072 00073 // ========================================================================== 00079 struct Cast_F 00080 { 00081 public : 00082 // ======================================================================== 00083 // prerequisites: 00084 BOOST_STATIC_ASSERT( std::numeric_limits<float> ::is_specialized && 00085 boost::integer_traits<int> ::is_specialized && 00086 boost::integer_traits<unsigned int> ::is_specialized && 00087 sizeof(float)==sizeof(int) && 00088 32 == boost::integer_traits<unsigned int>::digits ) ; 00089 // ======================================================================== 00090 public: 00091 // ======================================================================== 00093 float i2f ( const int i ) { m_f.i = i ; return m_f.f ; } // int -> float 00095 int f2i ( const float f ) { m_f.f = f ; return m_f.i ; } // float -> int 00096 // ======================================================================== 00097 private: 00098 // ======================================================================== 00100 union Float_U // Helper union to avoid reinterpret cast 00101 { 00102 float f ; // float value 00103 int i ; // int value 00104 } ; 00105 // ======================================================================== 00106 private: 00107 // ======================================================================== 00109 Float_U m_f ; // the helper union 00110 // ======================================================================== 00111 } ; 00112 // ========================================================================== 00118 struct Cast_D 00119 { 00120 // ======================================================================== 00122 BOOST_STATIC_ASSERT( std::numeric_limits<double> ::is_specialized && 00123 std::numeric_limits<Long> ::is_specialized && 00124 std::numeric_limits<ULong> ::is_specialized && 00125 boost::integer_traits<ULong> ::is_specialized && 00126 boost::integer_traits<Long> ::is_specialized && 00127 sizeof(double)==sizeof(Long) && 00128 sizeof(double)==sizeof(ULong) && 00129 64 == std::numeric_limits<ULong>::digits ) ; 00130 // ======================================================================== 00131 public: 00132 // ======================================================================== 00134 double l2d ( const Long l ) { m_d.l = l ; return m_d.d ; } // long -> double 00136 Long d2l ( const double d ) { m_d.d = d ; return m_d.l ; } // double -> long 00137 // ======================================================================== 00138 private: 00139 // ======================================================================== 00141 union Double_U // Helper union to avoid reinterpret cast 00142 { 00143 double d ; // double value 00144 Long l ; // long value 00145 } ; 00146 // ======================================================================== 00147 private: 00148 // ======================================================================== 00150 Double_U m_d ; // the helper union 00151 // ======================================================================== 00152 } ; 00153 // ========================================================================== 00154 } // end of anonymous namespace 00155 // ============================================================================ 00156 /* equality comparion of float numbers using as the metric the maximal 00157 * number of Units in the Last Place (ULP). 00158 * It is a slightly modified version of very efficient implementation 00159 * of the initial Bruce Dawson's algorithm by Chris Lomont. 00160 * 00161 * @see www.lomont.org 00162 * @see http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm 00163 * 00164 * Lomont claims the algorithm is factor 2-10 more efficient 00165 * with respect to Knuth's algorithm fomr comparions of floating number 00166 * using the relative precision. 00167 * 00168 * The effective relative difference depends on the choice of 00169 * <c>maxULPS</c>: 00170 * - For the case of maxULPs=1, (of cource it is totally unphysical case!!!) 00171 * the effective relative precision r = |a-b|/(|a|+|b|)is 00172 * between 3.5e-8 and 5.5e-8 for |a|,|b|>1.e-37, and 00173 * then it quickly goes to ~1 00174 * - For the case of maxULPS=10 00175 * the effective relative precision is 00176 * between 3e-8 and 6e-7 for |a|,|b|>1.e-37, and 00177 * then it quickly goes to ~1 00178 * - For the case of maxULPS=100 00179 * the effective relative precision is 00180 * around ~6e-6 for |a|,|b|>1.e-37, and 00181 * then it quickly goes to ~1 00182 * - For the case of maxULPS=1000 00183 * the effective relative precision is 00184 * around ~6e-5 for |a|,|b|>1.e-37, and 00185 * then it quickly goes to ~1 00186 * 00187 * @param af the first number 00188 * @param bf the second number 00189 * @param maxULPS the maximal metric deciation in the terms of 00190 * maximal number of units in the last place 00191 * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl 00192 */ 00193 // ============================================================================ 00194 bool Gaudi::Math::lomont_compare_float 00195 ( const float af , 00196 const float bf , 00197 const unsigned short maxULPs ) 00198 { 00199 // ========================================================================== 00200 // prerequisites: 00201 BOOST_STATIC_ASSERT( std::numeric_limits<float> ::is_specialized && 00202 boost::integer_traits<int> ::is_specialized && 00203 boost::integer_traits<unsigned int> ::is_specialized && 00204 sizeof(float)==sizeof(int) && 00205 sizeof(float)==sizeof(unsigned int) && 00206 32 == boost::integer_traits<unsigned int>::digits ) ; 00207 // ========================================================================== 00208 00209 Cast_F caster ; 00210 00211 //int ai = *reinterpret_cast<const int*>( &af ) ; 00212 //int bi = *reinterpret_cast<const int*>( &bf ) ; 00213 int ai = caster.f2i ( af ) ; 00214 int bi = caster.f2i ( bf ) ; 00215 00216 int test = (((unsigned int)(ai^bi))>>31)-1; 00217 00218 // assert ( (0==test) || ( boost::integer_traits<unsigned int>::const_max == test ) ) ; 00219 00220 int diff = ((( boost::integer_traits<int>::const_min - ai ) & (~test)) | ( ai& test )) - bi ; 00221 00222 int maxDiff_ = maxULPs ; 00223 00224 int v1 = maxDiff_ + diff ; 00225 int v2 = maxDiff_ - diff ; 00226 00227 return 0<=(v1|v2) ; 00228 } 00229 // ============================================================================ 00230 /* get the floating number that representation 00231 * is different with respect to the argument for 00232 * the certain number of "Units in the Last Position". 00233 * For ulps=1, it is just next float number, for ulps=-1 is is the 00234 * previous one. 00235 * 00236 * This routine is very convinient to test the parameter maxULPS for 00237 * the routine Gaudi::Math::lomont_compare 00238 * 00239 * @see Gaudi::Math::lomont_compare 00240 * @param af the reference number 00241 * @param ulps the bias 00242 * @return the biased float number (on distance "ulps") 00243 * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl 00244 * @date 2008-11-08 00245 */ 00246 // ============================================================================ 00247 float Gaudi::Math::next_float ( const float af , const short ulps ) 00248 { 00250 BOOST_STATIC_ASSERT( std::numeric_limits<float> ::is_specialized && 00251 std::numeric_limits<int> ::is_specialized && 00252 sizeof(float)==sizeof(int) ) ; 00253 // ========================================================================== 00254 Cast_F caster ; 00255 00256 // int ai = *reinterpret_cast<const int*>( &af ) ; 00257 int ai = caster.f2i ( af ) ; 00258 ai += ulps ; 00259 // return *reinterpret_cast<float*>(&ai) ; 00260 return caster.i2f ( ai ) ; 00261 } 00262 // ============================================================================ 00263 /* equality comparison of float numbers using as the metric the maximal 00264 * number of Units in the Last Place (ULP). 00265 * It is a slightly modified version of very efficient implementation 00266 * of the initial Bruce Dawson's algorithm by Chris Lomont. 00267 * 00268 * @see www.lomont.org 00269 * @see http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm 00270 * 00271 * C.Lomont claims the algorithm is factor 2-10 more efficient 00272 * with respect to Knuth's algorithm fomr comparions of floating number 00273 * using the relative precision. 00274 * 00275 * The effective relative difference depends on the choice of 00276 * <c>maxULPS</c>: 00277 * - For the case of maxULPs=1, (of cource it is totally unphysical case!!!) 00278 * the effective relative precision r = |a-b|/(|a|+|b|)is 00279 * ~6e-16 for |a|,|b|>1.e-304, and 00280 * then it quickly goes to ~1 00281 * 00282 * @param af the first number 00283 * @param bf the second number 00284 * @param maxULPS the maximal metric deciation in the terms of 00285 * maximal number of units in the last place 00286 * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl 00287 * @date 2008-11-08 00288 */ 00289 // ============================================================================ 00290 bool Gaudi::Math::lomont_compare_double 00291 ( const double af , 00292 const double bf , 00293 const unsigned int maxULPs ) 00294 { 00295 // ========================================================================== 00297 BOOST_STATIC_ASSERT( std::numeric_limits<double> ::is_specialized && 00298 std::numeric_limits<Long> ::is_specialized && 00299 std::numeric_limits<ULong> ::is_specialized && 00300 boost::integer_traits<ULong> ::is_specialized && 00301 boost::integer_traits<Long> ::is_specialized && 00302 sizeof(double)==sizeof(Long) && 00303 sizeof(double)==sizeof(ULong) && 00304 64 == std::numeric_limits<ULong>::digits ) ; 00305 // ========================================================================== 00306 Cast_D caster ; 00307 00308 //Long ai = *reinterpret_cast<const Long*>( &af ) ; 00309 //Long bi = *reinterpret_cast<const Long*>( &bf ) ; 00310 Long ai = caster.d2l ( af ) ; 00311 Long bi = caster.d2l ( bf ) ; 00312 00313 Long test = (((ULong)(ai^bi))>>63)-1; 00314 00315 // assert ( (0==test) || ( boost::integer_traits<ULong>::const_max == test ) ) ; 00316 00317 Long diff = ((( boost::integer_traits<Long>::const_min - ai ) & (~test)) | ( ai& test )) - bi ; 00318 00319 Long maxDiff_ = maxULPs ; 00320 00321 Long v1 = maxDiff_ + diff ; 00322 Long v2 = maxDiff_ - diff ; 00323 00324 return 0<=(v1|v2) ; 00325 } 00326 // ============================================================================ 00327 /* Get the floating number that representation 00328 * is different with respect to the argument for 00329 * the certain number of "Units in the Last Position". 00330 * For ulps=1, it is just next float number, for ulps=-1 is is the 00331 * previous one. 00332 * 00333 * This routine is very convinient to test the parameter maxULPS for 00334 * the routine LHCb::Math::lomont_compare_float 00335 * 00336 * @see Gaudi::Math::lomont_compare 00337 * @param af the reference number 00338 * @param ulps the bias 00339 * @return the biased float number (on distance "ulps") 00340 * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl 00341 * @date 2008-11-08 00342 */ 00343 // ============================================================================ 00344 double Gaudi::Math::next_double ( const double ad , const short ulps ) 00345 { 00346 // ========================================================================== 00348 BOOST_STATIC_ASSERT( std::numeric_limits<double> ::is_specialized && 00349 std::numeric_limits<Long> ::is_specialized && 00350 sizeof(double)==sizeof(Long) ) ; 00351 // ========================================================================== 00352 Cast_D caster ; 00353 00354 //Long al = *reinterpret_cast<const Long*>( &ad ) ; 00355 Long al = caster.d2l ( ad ) ; 00356 al += ulps ; 00357 //return *reinterpret_cast<double*>(&al) ; 00358 return caster.l2d ( al ) ; 00359 } 00360 // ============================================================================ 00361 00362 00363 // ============================================================================ 00364 // The END 00365 // ============================================================================ 00366