| The Gaudi Framework v27r1

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

Lomont.cpp
Go to the documentation of this file.
1 // ============================================================================
2 // Include files
3 // ============================================================================
4 // STD & STL
5 // ============================================================================
6 #include <iostream>
7 #include <limits>
8 #include <cassert>
9 #include <cmath>
10 #include <vector>
11 // ============================================================================
12 // GaudiKernel
13 // ============================================================================
14 #include "GaudiKernel/Lomont.h"
15 // ============================================================================
16 // Boost
17 // ============================================================================
18 #include "boost/integer_traits.hpp"
19 // ============================================================================
20 namespace
21 {
22 // prerequisites for "correct" Float
23 static_assert( std::numeric_limits<float> ::is_specialized &&
24 boost::integer_traits<int> ::is_specialized &&
25 boost::integer_traits<unsigned int> ::is_specialized &&
26 sizeof(float)==sizeof(int) &&
27 sizeof(float)==sizeof(unsigned int) &&
28 32 == boost::integer_traits<unsigned int>::digits , "FAILED ASSUMPTIONS") ;
29 // ==========================================================================
30 // define proepr double
31 // ==========================================================================
32 #ifdef _WIN32
33 struct _Longs
34 {
35 typedef __int64 Long ;
36 typedef unsigned __int64 ULong ;
37 } ;
38 #else
39 template <bool I>
40 struct __Longs ;
41 template <>
42 struct __Longs<true>
43 {
44 typedef long Long ;
45 typedef unsigned long ULong ;
46 } ;
47 template <>
48 struct __Longs<false>
49 {
50 typedef long long Long ;
51 typedef unsigned long long ULong ;
52 } ;
53 struct _Longs : public __Longs<sizeof(double)==sizeof(long)>{};
54 #endif
55
57 typedef _Longs::Long Long ;
58 typedef _Longs::ULong ULong ;
59 // ==========================================================================
61 static_assert( std::numeric_limits<double> ::is_specialized &&
62 std::numeric_limits<Long> ::is_specialized &&
63 std::numeric_limits<ULong> ::is_specialized &&
64 boost::integer_traits<ULong> ::is_specialized &&
65 boost::integer_traits<Long> ::is_specialized &&
66 sizeof(double)==sizeof(Long) &&
67 sizeof(double)==sizeof(ULong) &&
68 64 == std::numeric_limits<ULong>::digits , "FAILED ASSUMPTIONS") ;
69 // ==========================================================================
70
71 // ==========================================================================
77 struct Cast_F
78 {
79 public :
80 // ========================================================================
81 // prerequisites:
82 static_assert( std::numeric_limits<float> ::is_specialized &&
83 boost::integer_traits<int> ::is_specialized &&
84 boost::integer_traits<unsigned int> ::is_specialized &&
85 sizeof(float)==sizeof(int) &&
86 32 == boost::integer_traits<unsigned int>::digits , "FAILED ASSUMPTIONS") ;
87 // ========================================================================
88 public:
89 // ========================================================================
91 float i2f ( const int i ) { m_f.i = i ; return m_f.f ; } // int -> float
93 int f2i ( const float f ) { m_f.f = f ; return m_f.i ; } // float -> int
94 // ========================================================================
95 private:
96 // ========================================================================
98 union Float_U // Helper union to avoid reinterpret cast
99 {
100 float f ; // float value
101 int i ; // int value
102 } ;
103 // ========================================================================
104 private:
105 // ========================================================================
107 Float_U m_f ; // the helper union
108 // ========================================================================
109 } ;
110 // ==========================================================================
116 struct Cast_D
117 {
118 // ========================================================================
120 static_assert( std::numeric_limits<double> ::is_specialized &&
121 std::numeric_limits<Long> ::is_specialized &&
122 std::numeric_limits<ULong> ::is_specialized &&
123 boost::integer_traits<ULong> ::is_specialized &&
124 boost::integer_traits<Long> ::is_specialized &&
125 sizeof(double)==sizeof(Long) &&
126 sizeof(double)==sizeof(ULong) &&
127 64 == std::numeric_limits<ULong>::digits , "FAILED ASSUMPTIONS") ;
128 // ========================================================================
129 public:
130 // ========================================================================
132 double l2d ( const Long l ) { m_d.l = l ; return m_d.d ; } // long -> double
134 Long d2l ( const double d ) { m_d.d = d ; return m_d.l ; } // double -> long
135 // ========================================================================
136 private:
137 // ========================================================================
139 union Double_U // Helper union to avoid reinterpret cast
140 {
141 double d ; // double value
142 Long l ; // long value
143 } ;
144 // ========================================================================
145 private:
146 // ========================================================================
148 Double_U m_d ; // the helper union
149 // ========================================================================
150 } ;
151 // ==========================================================================
152 } // end of anonymous namespace
153 // ============================================================================
154 /* equality comparion of float numbers using as the metric the maximal
155 * number of Units in the Last Place (ULP).
156 * It is a slightly modified version of very efficient implementation
157 * of the initial Bruce Dawson's algorithm by Chris Lomont.
158 *
159 * @see www.lomont.org
160 * @see http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
161 *
162 * Lomont claims the algorithm is factor 2-10 more efficient
163 * with respect to Knuth's algorithm fomr comparions of floating number
164 * using the relative precision.
165 *
166 * The effective relative difference depends on the choice of
167 * <c>maxULPS</c>:
168 * - For the case of maxULPs=1, (of cource it is totally unphysical case!!!)
169 * the effective relative precision r = |a-b|/(|a|+|b|)is
170 * between 3.5e-8 and 5.5e-8 for |a|,|b|>1.e-37, and
171 * then it quickly goes to ~1
172 * - For the case of maxULPS=10
173 * the effective relative precision is
174 * between 3e-8 and 6e-7 for |a|,|b|>1.e-37, and
175 * then it quickly goes to ~1
176 * - For the case of maxULPS=100
177 * the effective relative precision is
178 * around ~6e-6 for |a|,|b|>1.e-37, and
179 * then it quickly goes to ~1
180 * - For the case of maxULPS=1000
181 * the effective relative precision is
182 * around ~6e-5 for |a|,|b|>1.e-37, and
183 * then it quickly goes to ~1
184 *
185 * @param af the first number
186 * @param bf the second number
187 * @param maxULPS the maximal metric deciation in the terms of
188 * maximal number of units in the last place
189 * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl
190 */
191 // ============================================================================
192 bool Gaudi::Math::lomont_compare_float
193 ( const float af ,
194 const float bf ,
195 const unsigned short maxULPs )
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 }
227 // ============================================================================
228 /* get the floating number that representation
229 * is different with respect to the argument for
230 * the certain number of "Units in the Last Position".
231 * For ulps=1, it is just next float number, for ulps=-1 is is the
232 * previous one.
233 *
234 * This routine is very convinient to test the parameter maxULPS for
235 * the routine Gaudi::Math::lomont_compare
236 *
237 * @see Gaudi::Math::lomont_compare
238 * @param af the reference number
239 * @param ulps the bias
240 * @return the biased float number (on distance "ulps")
241 * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl
242 * @date 2008-11-08
243 */
244 // ============================================================================
245 float Gaudi::Math::next_float ( const float af , const short ulps )
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 }
260 // ============================================================================
261 /* equality comparison of float numbers using as the metric the maximal
262 * number of Units in the Last Place (ULP).
263 * It is a slightly modified version of very efficient implementation
264 * of the initial Bruce Dawson's algorithm by Chris Lomont.
265 *
266 * @see www.lomont.org
267 * @see http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
268 *
269 * C.Lomont claims the algorithm is factor 2-10 more efficient
270 * with respect to Knuth's algorithm fomr comparions of floating number
271 * using the relative precision.
272 *
273 * The effective relative difference depends on the choice of
274 * <c>maxULPS</c>:
275 * - For the case of maxULPs=1, (of cource it is totally unphysical case!!!)
276 * the effective relative precision r = |a-b|/(|a|+|b|)is
277 * ~6e-16 for |a|,|b|>1.e-304, and
278 * then it quickly goes to ~1
279 *
280 * @param af the first number
281 * @param bf the second number
282 * @param maxULPS the maximal metric deciation in the terms of
283 * maximal number of units in the last place
284 * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl
285 * @date 2008-11-08
286 */
287 // ============================================================================
288 bool Gaudi::Math::lomont_compare_double
289 ( const double af ,
290 const double bf ,
291 const unsigned int maxULPs )
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 }
324 // ============================================================================
325 /* Get the floating number that representation
326 * is different with respect to the argument for
327 * the certain number of "Units in the Last Position".
328 * For ulps=1, it is just next float number, for ulps=-1 is is the
329 * previous one.
330 *
331 * This routine is very convinient to test the parameter maxULPS for
332 * the routine LHCb::Math::lomont_compare_float
333 *
334 * @see Gaudi::Math::lomont_compare
335 * @param af the reference number
336 * @param ulps the bias
337 * @return the biased float number (on distance "ulps")
338 * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl
339 * @date 2008-11-08
340 */
341 // ============================================================================
342 double Gaudi::Math::next_double ( const double ad , const short ulps )
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 }
358 // ============================================================================
359
360
361 // ============================================================================
362 // The END
363 // ============================================================================
364
Gaudi::Math::lomont_compare_float
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 Plac...
Definition: Lomont.cpp:193
Gaudi::Math::next_float
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...
Definition: Lomont.cpp:245
compareRootHistos.test
string test
Definition: compareRootHistos.py:15
Gaudi::Math::lomont_compare_double
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 Pla...
Definition: Lomont.cpp:289
gaudirun.l
dictionary l
Definition: gaudirun.py:421
Lomont.h
ana.i
list i
Definition: ana.py:128
Gaudi::Math::next_double
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...
Definition: Lomont.cpp:342
std::numeric_limits