| # Gaudi Framework, version v25r0 |

Home |  | Generated: Mon Feb 17 2014

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

- File List
- File Members

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