Gaudi Framework, version v23r6
Home   Generated: Wed Jan 30 2013
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
Lomont.cpp
Go to the documentation of this file.
1 // $Id: Lomont.cpp,v 1.2 2009/05/23 17:08:16 ibelyaev Exp $
2 // ============================================================================
3 // Include files
4 // ============================================================================
5 // STD & STL
6 // ============================================================================
7 #include <iostream>
8 #include <limits>
9 #include <cassert>
10 #include <cmath>
11 #include <vector>
12 // ============================================================================
13 // GaudiKernel
14 // ============================================================================
15 #include "GaudiKernel/Lomont.h"
16 // ============================================================================
17 // Boost
18 // ============================================================================
19 #include "boost/static_assert.hpp"
20 #include "boost/integer_traits.hpp"
21 // ============================================================================
22 namespace
23 {
24  // prerequisites for "correct" Float
25  BOOST_STATIC_ASSERT( std::numeric_limits<float> ::is_specialized &&
26  boost::integer_traits<int> ::is_specialized &&
27  boost::integer_traits<unsigned int> ::is_specialized &&
28  sizeof(float)==sizeof(int) &&
29  sizeof(float)==sizeof(unsigned int) &&
30  32 == boost::integer_traits<unsigned int>::digits ) ;
31  // ==========================================================================
32  // define proepr double
33  // ==========================================================================
34 #ifdef _WIN32
35  struct _Longs
36  {
37  typedef __int64 Long ;
38  typedef unsigned __int64 ULong ;
39  } ;
40 #else
41  template <bool I>
42  struct __Longs ;
43  template <>
44  struct __Longs<true>
45  {
46  typedef long Long ;
47  typedef unsigned long ULong ;
48  } ;
49  template <>
50  struct __Longs<false>
51  {
52  typedef long long Long ;
53  typedef unsigned long long ULong ;
54  } ;
55  struct _Longs : public __Longs<sizeof(double)==sizeof(long)>{};
56 #endif
57 
59  typedef _Longs::Long Long ;
60  typedef _Longs::ULong ULong ;
61  // ==========================================================================
63  BOOST_STATIC_ASSERT( std::numeric_limits<double> ::is_specialized &&
64  std::numeric_limits<Long> ::is_specialized &&
65  std::numeric_limits<ULong> ::is_specialized &&
66  boost::integer_traits<ULong> ::is_specialized &&
67  boost::integer_traits<Long> ::is_specialized &&
68  sizeof(double)==sizeof(Long) &&
69  sizeof(double)==sizeof(ULong) &&
70  64 == std::numeric_limits<ULong>::digits ) ;
71  // ==========================================================================
72 
73  // ==========================================================================
79  struct Cast_F
80  {
81  public :
82  // ========================================================================
83  // prerequisites:
84  BOOST_STATIC_ASSERT( std::numeric_limits<float> ::is_specialized &&
85  boost::integer_traits<int> ::is_specialized &&
86  boost::integer_traits<unsigned int> ::is_specialized &&
87  sizeof(float)==sizeof(int) &&
88  32 == boost::integer_traits<unsigned int>::digits ) ;
89  // ========================================================================
90  public:
91  // ========================================================================
93  float i2f ( const int i ) { m_f.i = i ; return m_f.f ; } // int -> float
95  int f2i ( const float f ) { m_f.f = f ; return m_f.i ; } // float -> int
96  // ========================================================================
97  private:
98  // ========================================================================
100  union Float_U // Helper union to avoid reinterpret cast
101  {
102  float f ; // float value
103  int i ; // int value
104  } ;
105  // ========================================================================
106  private:
107  // ========================================================================
109  Float_U m_f ; // the helper union
110  // ========================================================================
111  } ;
112  // ==========================================================================
118  struct Cast_D
119  {
120  // ========================================================================
122  BOOST_STATIC_ASSERT( std::numeric_limits<double> ::is_specialized &&
123  std::numeric_limits<Long> ::is_specialized &&
124  std::numeric_limits<ULong> ::is_specialized &&
125  boost::integer_traits<ULong> ::is_specialized &&
126  boost::integer_traits<Long> ::is_specialized &&
127  sizeof(double)==sizeof(Long) &&
128  sizeof(double)==sizeof(ULong) &&
129  64 == std::numeric_limits<ULong>::digits ) ;
130  // ========================================================================
131  public:
132  // ========================================================================
134  double l2d ( const Long l ) { m_d.l = l ; return m_d.d ; } // long -> double
136  Long d2l ( const double d ) { m_d.d = d ; return m_d.l ; } // double -> long
137  // ========================================================================
138  private:
139  // ========================================================================
141  union Double_U // Helper union to avoid reinterpret cast
142  {
143  double d ; // double value
144  Long l ; // long value
145  } ;
146  // ========================================================================
147  private:
148  // ========================================================================
150  Double_U m_d ; // the helper union
151  // ========================================================================
152  } ;
153  // ==========================================================================
154 } // end of anonymous namespace
155 // ============================================================================
156 /* equality comparion of float numbers using as the metric the maximal
157  * number of Units in the Last Place (ULP).
158  * It is a slightly modified version of very efficient implementation
159  * of the initial Bruce Dawson's algorithm by Chris Lomont.
160  *
161  * @see www.lomont.org
162  * @see http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
163  *
164  * Lomont claims the algorithm is factor 2-10 more efficient
165  * with respect to Knuth's algorithm fomr comparions of floating number
166  * using the relative precision.
167  *
168  * The effective relative difference depends on the choice of
169  * <c>maxULPS</c>:
170  * - For the case of maxULPs=1, (of cource it is totally unphysical case!!!)
171  * the effective relative precision r = |a-b|/(|a|+|b|)is
172  * between 3.5e-8 and 5.5e-8 for |a|,|b|>1.e-37, and
173  * then it quickly goes to ~1
174  * - For the case of maxULPS=10
175  * the effective relative precision is
176  * between 3e-8 and 6e-7 for |a|,|b|>1.e-37, and
177  * then it quickly goes to ~1
178  * - For the case of maxULPS=100
179  * the effective relative precision is
180  * around ~6e-6 for |a|,|b|>1.e-37, and
181  * then it quickly goes to ~1
182  * - For the case of maxULPS=1000
183  * the effective relative precision is
184  * around ~6e-5 for |a|,|b|>1.e-37, and
185  * then it quickly goes to ~1
186  *
187  * @param af the first number
188  * @param bf the second number
189  * @param maxULPS the maximal metric deciation in the terms of
190  * maximal number of units in the last place
191  * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl
192  */
193 // ============================================================================
194 bool Gaudi::Math::lomont_compare_float
195 ( const float af ,
196  const float bf ,
197  const unsigned short maxULPs )
198 {
199  // ==========================================================================
200  // prerequisites:
201  BOOST_STATIC_ASSERT( std::numeric_limits<float> ::is_specialized &&
202  boost::integer_traits<int> ::is_specialized &&
203  boost::integer_traits<unsigned int> ::is_specialized &&
204  sizeof(float)==sizeof(int) &&
205  sizeof(float)==sizeof(unsigned int) &&
206  32 == boost::integer_traits<unsigned int>::digits ) ;
207  // ==========================================================================
208 
209  Cast_F caster ;
210 
211  //int ai = *reinterpret_cast<const int*>( &af ) ;
212  //int bi = *reinterpret_cast<const int*>( &bf ) ;
213  int ai = caster.f2i ( af ) ;
214  int bi = caster.f2i ( bf ) ;
215 
216  int test = (((unsigned int)(ai^bi))>>31)-1;
217 
218  // assert ( (0==test) || ( boost::integer_traits<unsigned int>::const_max == test ) ) ;
219 
220  int diff = ((( boost::integer_traits<int>::const_min - ai ) & (~test)) | ( ai& test )) - bi ;
221 
222  int maxDiff_ = maxULPs ;
223 
224  int v1 = maxDiff_ + diff ;
225  int v2 = maxDiff_ - diff ;
226 
227  return 0<=(v1|v2) ;
228 }
229 // ============================================================================
230 /* get the floating number that representation
231  * is different with respect to the argument for
232  * the certain number of "Units in the Last Position".
233  * For ulps=1, it is just next float number, for ulps=-1 is is the
234  * previous one.
235  *
236  * This routine is very convinient to test the parameter maxULPS for
237  * the routine Gaudi::Math::lomont_compare
238  *
239  * @see Gaudi::Math::lomont_compare
240  * @param af the reference number
241  * @param ulps the bias
242  * @return the biased float number (on distance "ulps")
243  * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl
244  * @date 2008-11-08
245  */
246 // ============================================================================
247 float Gaudi::Math::next_float ( const float af , const short ulps )
248 {
250  BOOST_STATIC_ASSERT( std::numeric_limits<float> ::is_specialized &&
251  std::numeric_limits<int> ::is_specialized &&
252  sizeof(float)==sizeof(int) ) ;
253  // ==========================================================================
254  Cast_F caster ;
255 
256  // int ai = *reinterpret_cast<const int*>( &af ) ;
257  int ai = caster.f2i ( af ) ;
258  ai += ulps ;
259  // return *reinterpret_cast<float*>(&ai) ;
260  return caster.i2f ( ai ) ;
261 }
262 // ============================================================================
263 /* equality comparison of float numbers using as the metric the maximal
264  * number of Units in the Last Place (ULP).
265  * It is a slightly modified version of very efficient implementation
266  * of the initial Bruce Dawson's algorithm by Chris Lomont.
267  *
268  * @see www.lomont.org
269  * @see http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
270  *
271  * C.Lomont claims the algorithm is factor 2-10 more efficient
272  * with respect to Knuth's algorithm fomr comparions of floating number
273  * using the relative precision.
274  *
275  * The effective relative difference depends on the choice of
276  * <c>maxULPS</c>:
277  * - For the case of maxULPs=1, (of cource it is totally unphysical case!!!)
278  * the effective relative precision r = |a-b|/(|a|+|b|)is
279  * ~6e-16 for |a|,|b|>1.e-304, and
280  * then it quickly goes to ~1
281  *
282  * @param af the first number
283  * @param bf the second number
284  * @param maxULPS the maximal metric deciation in the terms of
285  * maximal number of units in the last place
286  * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl
287  * @date 2008-11-08
288  */
289 // ============================================================================
290 bool Gaudi::Math::lomont_compare_double
291 ( const double af ,
292  const double bf ,
293  const unsigned int maxULPs )
294 {
295  // ==========================================================================
297  BOOST_STATIC_ASSERT( std::numeric_limits<double> ::is_specialized &&
298  std::numeric_limits<Long> ::is_specialized &&
299  std::numeric_limits<ULong> ::is_specialized &&
300  boost::integer_traits<ULong> ::is_specialized &&
301  boost::integer_traits<Long> ::is_specialized &&
302  sizeof(double)==sizeof(Long) &&
303  sizeof(double)==sizeof(ULong) &&
304  64 == std::numeric_limits<ULong>::digits ) ;
305  // ==========================================================================
306  Cast_D caster ;
307 
308  //Long ai = *reinterpret_cast<const Long*>( &af ) ;
309  //Long bi = *reinterpret_cast<const Long*>( &bf ) ;
310  Long ai = caster.d2l ( af ) ;
311  Long bi = caster.d2l ( bf ) ;
312 
313  Long test = (((ULong)(ai^bi))>>63)-1;
314 
315  // assert ( (0==test) || ( boost::integer_traits<ULong>::const_max == test ) ) ;
316 
317  Long diff = ((( boost::integer_traits<Long>::const_min - ai ) & (~test)) | ( ai& test )) - bi ;
318 
319  Long maxDiff_ = maxULPs ;
320 
321  Long v1 = maxDiff_ + diff ;
322  Long v2 = maxDiff_ - diff ;
323 
324  return 0<=(v1|v2) ;
325 }
326 // ============================================================================
327 /* Get the floating number that representation
328  * is different with respect to the argument for
329  * the certain number of "Units in the Last Position".
330  * For ulps=1, it is just next float number, for ulps=-1 is is the
331  * previous one.
332  *
333  * This routine is very convinient to test the parameter maxULPS for
334  * the routine LHCb::Math::lomont_compare_float
335  *
336  * @see Gaudi::Math::lomont_compare
337  * @param af the reference number
338  * @param ulps the bias
339  * @return the biased float number (on distance "ulps")
340  * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl
341  * @date 2008-11-08
342  */
343 // ============================================================================
344 double Gaudi::Math::next_double ( const double ad , const short ulps )
345 {
346  // ==========================================================================
348  BOOST_STATIC_ASSERT( std::numeric_limits<double> ::is_specialized &&
349  std::numeric_limits<Long> ::is_specialized &&
350  sizeof(double)==sizeof(Long) ) ;
351  // ==========================================================================
352  Cast_D caster ;
353 
354  //Long al = *reinterpret_cast<const Long*>( &ad ) ;
355  Long al = caster.d2l ( ad ) ;
356  al += ulps ;
357  //return *reinterpret_cast<double*>(&al) ;
358  return caster.l2d ( al ) ;
359 }
360 // ============================================================================
361 
362 
363 // ============================================================================
364 // The END
365 // ============================================================================
366 
Generated at Wed Jan 30 2013 17:13:40 for Gaudi Framework, version v23r6 by Doxygen version 1.8.2 written by Dimitri van Heesch, © 1997-2004