Elements  5.8
A C++ base framework for the Euclid Software.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
Real.h
Go to the documentation of this file.
1 
40 // Copyright 2005, Google Inc.
41 // All rights reserved.
42 //
43 // Redistribution and use in source and binary forms, with or without
44 // modification, are permitted provided that the following conditions are
45 // met:
46 //
47 // * Redistributions of source code must retain the above copyright
48 // notice, this list of conditions and the following disclaimer.
49 // * Redistributions in binary form must reproduce the above
50 // copyright notice, this list of conditions and the following disclaimer
51 // in the documentation and/or other materials provided with the
52 // distribution.
53 // * Neither the name of Google Inc. nor the names of its
54 // contributors may be used to endorse or promote products derived from
55 // this software without specific prior written permission.
56 //
57 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
58 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
59 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
60 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
61 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
62 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
63 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
64 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
65 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
66 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
67 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
68 //
69 // Authors: wan@google.com (Zhanyong Wan), eefacm@gmail.com (Sean Mcafee)
70 //
71 // The Google C++ Testing Framework (Google Test)
72 #ifndef ELEMENTSKERNEL_ELEMENTSKERNEL_REAL_H_
73 #define ELEMENTSKERNEL_ELEMENTSKERNEL_REAL_H_
74 
75 #include <limits> // for numeric_limits
76 #include <type_traits> // for is_floating_point
77 #include <cmath> // for round
78 
79 #include "ElementsKernel/Export.h" // ELEMENTS_API
80 #include "ElementsKernel/Unused.h" // ELEMENTS_UNUSED
81 
83 
84 namespace Elements {
85 
90 
91 // For testing purposes only. Rather use the isEqual functions for real
92 // life comparison
94 ELEMENTS_API extern const double FLT_DEFAULT_TEST_TOLERANCE;
96 ELEMENTS_API extern const double DBL_DEFAULT_TEST_TOLERANCE;
97 
98 
99 template<std::size_t size>
101 public:
102  // This prevents the user from using TypeWithSize<N> with incorrect
103  // values of N.
104  using UInt = void;
105 };
106 
107 // The specialisation for size 4.
108 template<>
110 public:
111  // unsigned int has size 4 in both gcc and MSVC.
112  //
113  // As base/basictypes.h doesn't compile on Windows, we cannot use
114  // uint32, uint64, and etc here.
115  using Int = int;
116  using UInt = unsigned int;
117 };
118 
119 // The specialisation for size 8.
120 template<>
122 public:
123  using Int = long long; // NOLINT
124  using UInt = unsigned long long; // NOLINT
125 };
126 
127 template <typename RawType>
129  return FLT_DEFAULT_MAX_ULPS;
130 }
131 
132 template <>
134  return FLT_DEFAULT_MAX_ULPS;
135 }
136 
137 template <>
139  return DBL_DEFAULT_MAX_ULPS;
140 }
141 
142 
143 // This template class represents an IEEE floating-point number
144 // (either single-precision or double-precision, depending on the
145 // template parameters).
146 //
147 // The purpose of this class is to do more sophisticated number
148 // comparison. (Due to round-off error, etc, it's very unlikely that
149 // two floating-points will be equal exactly. Hence a naive
150 // comparison by the == operation often doesn't work.)
151 //
152 // Format of IEEE floating-point:
153 //
154 // The most-significant bit being the leftmost, an IEEE
155 // floating-point looks like
156 //
157 // sign_bit exponent_bits fraction_bits
158 //
159 // Here, sign_bit is a single bit that designates the sign of the
160 // number.
161 //
162 // For float, there are 8 exponent bits and 23 fraction bits.
163 //
164 // For double, there are 11 exponent bits and 52 fraction bits.
165 //
166 // More details can be found at
167 // http://en.wikipedia.org/wiki/IEEE_floating-point_standard.
168 //
169 // Template parameter:
170 //
171 // RawType: the raw floating-point type (either float or double)
172 template<typename RawType>
174 public:
175  // Defines the unsigned integer type that has the same size as the
176  // floating point number.
178 
179  // Constants.
180 
181  // # of bits in a number.
182  static const std::size_t s_bitcount = 8 * sizeof(RawType);
183 
184  // # of fraction bits in a number.
185  static const std::size_t s_fraction_bitcount = std::numeric_limits<RawType>::digits - 1;
186 
187  // # of exponent bits in a number.
188  static const std::size_t s_exponent_bitcount = s_bitcount - 1 - s_fraction_bitcount;
189 
190  // The mask for the sign bit.
191  static const Bits s_sign_bitmask = static_cast<Bits>(1) << (s_bitcount - 1);
192 
193  // The mask for the fraction bits.
194  static const Bits s_fraction_bitmask = ~static_cast<Bits>(0) >> (s_exponent_bitcount + 1);
195 
196  // The mask for the exponent bits.
197  static const Bits s_exponent_bitmask = ~(s_sign_bitmask | s_fraction_bitmask);
198 
199  // How many ULP's (Units in the Last Place) we want to tolerate when
200  // comparing two numbers. The larger the value, the more error we
201  // allow. A 0 value means that two numbers must be exactly the same
202  // to be considered equal.
203  //
204  // The maximum error of a single floating-point operation is 0.5
205  // units in the last place. On Intel CPU's, all floating-point
206  // calculations are done with 80-bit precision, while double has 64
207  // bits. Therefore, 4 should be enough for ordinary use.
208  //
209  // See the following article for more details on ULP:
210  // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm.
211  static const std::size_t m_max_ulps = defaultMaxUlps<RawType>();
212 
213  // Constructs a FloatingPoint from a raw floating-point number.
214  //
215  // On an Intel CPU, passing a non-normalised NAN (Not a Number)
216  // around may change its bits, although the new value is guaranteed
217  // to be also a NAN. Therefore, don't expect this constructor to
218  // preserve the bits in x when x is a NAN.
219  explicit FloatingPoint(const RawType& x) {
220  m_u.m_value = x;
221  }
222 
223  // Static methods
224 
225  // Reinterprets a bit pattern as a floating-point number.
226  //
227  // This function is needed to test the AlmostEquals() method.
228  static RawType ReinterpretBits(const Bits &bits) {
229  FloatingPoint fp(0);
230  fp.m_u.m_bits = bits;
231  return fp.m_u.m_value;
232  }
233 
234  // Returns the floating-point number that represent positive infinity.
235  static RawType Infinity() {
236  return ReinterpretBits(s_exponent_bitmask);
237  }
238 
239  // Non-static methods
240 
241  // Returns the bits that represents this number.
242  const Bits &bits() const {
243  return m_u.m_bits;
244  }
245 
246  // Returns the exponent bits of this number.
247  Bits exponentBits() const {
248  return s_exponent_bitmask & m_u.m_bits;
249  }
250 
251  // Returns the fraction bits of this number.
252  Bits fractionBits() const {
253  return s_fraction_bitmask & m_u.m_bits;
254  }
255 
256  // Returns the sign bit of this number.
257  Bits signBit() const {
258  return s_sign_bitmask & m_u.m_bits;
259  }
260 
261  // Returns true iff this is NAN (not a number).
262  bool isNan() const {
263  // It's a NAN if the exponent bits are all ones and the fraction
264  // bits are not entirely zeros.
265  return (exponentBits() == s_exponent_bitmask) && (fractionBits() != 0);
266  }
267 
268  // Returns true iff this number is at most kMaxUlps ULP's away from
269  // rhs. In particular, this function:
270  //
271  // - returns false if either number is (or both are) NAN.
272  // - treats really large numbers as almost equal to infinity.
273  // - thinks +0.0 and -0.0 are 0 DLP's apart.
274  bool AlmostEquals(const FloatingPoint& rhs) const {
275  // The IEEE standard says that any comparison operation involving
276  // a NAN must return false.
277  if (isNan() || rhs.isNan()) {
278  return false;
279  }
280  return distanceBetweenSignAndMagnitudeNumbers(m_u.m_bits, rhs.m_u.m_bits) <= m_max_ulps;
281  }
282 
283  // Converts an integer from the sign-and-magnitude representation to
284  // the biased representation. More precisely, let N be 2 to the
285  // power of (kBitCount - 1), an integer x is represented by the
286  // unsigned number x + N.
287  //
288  // For instance,
289  //
290  // -N + 1 (the most negative number representable using
291  // sign-and-magnitude) is represented by 1;
292  // 0 is represented by N; and
293  // N - 1 (the biggest number representable using
294  // sign-and-magnitude) is represented by 2N - 1.
295  //
296  // Read http://en.wikipedia.org/wiki/Signed_number_representations
297  // for more details on signed number representations.
298  static Bits signAndMagnitudeToBiased(const Bits &sam) {
299  if (s_sign_bitmask & sam) {
300  // sam represents a negative number.
301  return ~sam + 1;
302  } else {
303  // sam represents a positive number.
304  return s_sign_bitmask | sam;
305  }
306  }
307 
308  // Given two numbers in the sign-and-magnitude representation,
309  // returns the distance between them as an unsigned number.
310  static Bits distanceBetweenSignAndMagnitudeNumbers(const Bits &sam1, const Bits &sam2) {
311  const Bits biased1 = signAndMagnitudeToBiased(sam1);
312  const Bits biased2 = signAndMagnitudeToBiased(sam2);
313  return (biased1 >= biased2) ? (biased1 - biased2) : (biased2 - biased1);
314  }
315 
316 private:
317  // The data type used to store the actual floating-point number.
319  RawType m_value; // The raw floating-point number.
320  Bits m_bits; // The bits that represent the number.
321  };
322 
324 };
325 
326 
327 // Usable AlmostEqual function
328 
329 template <typename FloatType>
330 bool almostEqual2sComplement(ELEMENTS_UNUSED const FloatType& a,
331  ELEMENTS_UNUSED const FloatType& b,
332  ELEMENTS_UNUSED const std::size_t& max_ulps = 0) {
333  return false;
334 }
335 
336 
337 template <typename RawType>
338 bool isNan(const RawType& x) {
339 
340  using Bits = typename TypeWithSize<sizeof(RawType)>::UInt;
341  Bits x_bits = *reinterpret_cast<const Bits *>(&x);
342 
343  Bits x_exp_bits = FloatingPoint<RawType>::s_exponent_bitmask & x_bits;
344  Bits x_frac_bits = FloatingPoint<RawType>::s_fraction_bitmask & x_bits;
345 
346  return (x_exp_bits == FloatingPoint<RawType>::s_exponent_bitmask) && (x_frac_bits != 0);
347 
348 }
349 
350 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
351 bool isEqual(const RawType& left, const RawType& right) {
352 
353  bool is_equal{false};
354 
355  if (not(isNan<RawType>(left) or isNan<RawType>(right))) {
356  using Bits = typename TypeWithSize<sizeof(RawType)>::UInt;
357  Bits l_bits = *reinterpret_cast<const Bits *>(&left);
358  Bits r_bits = *reinterpret_cast<const Bits *>(&right);
359  is_equal = (FloatingPoint<RawType>::distanceBetweenSignAndMagnitudeNumbers(l_bits, r_bits) <= max_ulps);
360  }
361 
362  return is_equal;
363 }
364 
365 template <std::size_t max_ulps>
366 inline bool isEqual(const float& left, const float& right) {
367  return (isEqual<float, max_ulps>(left, right));
368 }
369 
370 template <std::size_t max_ulps>
371 inline bool isEqual(const double& left, const double& right) {
372  return (isEqual<double, max_ulps>(left, right));
373 }
374 
375 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
376 inline bool isNotEqual(const RawType& left, const RawType& right) {
377  return (not isEqual<RawType, max_ulps>(left, right) );
378 }
379 
380 template <std::size_t max_ulps>
381 inline bool isNotEqual(const float& left, const float& right) {
382  return (isNotEqual<float, max_ulps>(left, right));
383 }
384 
385 template <std::size_t max_ulps>
386 inline bool isNotEqual(const double& left, const double& right) {
387  return (isNotEqual<double, max_ulps>(left, right));
388 }
389 
390 
391 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
392 bool isLess(const RawType& left, const RawType& right) {
393  bool is_less{false};
394 
395  if ( left < right && (not isEqual<RawType, max_ulps>(left, right)) ) {
396  is_less = true;
397  }
398 
399  return is_less;
400 }
401 
402 template<std::size_t max_ulps>
403 inline bool isLess(const float& left, const float& right) {
404  return (isLess<float, max_ulps>(left, right));
405 }
406 
407 template<std::size_t max_ulps>
408 inline bool isLess(const double& left, const double& right) {
409  return (isLess<double, max_ulps>(left, right));
410 }
411 
412 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
413 bool isGreater(const RawType& left, const RawType& right) {
414  bool is_greater{false};
415 
416  if ( left > right && (not isEqual<RawType, max_ulps>(left, right)) ) {
417  is_greater = true;
418  }
419 
420  return is_greater;
421 }
422 
423 template<std::size_t max_ulps>
424 inline bool isGreater(const float& left, const float& right) {
425  return (isGreater<float, max_ulps>(left, right));
426 }
427 
428 template<std::size_t max_ulps>
429 inline bool isGreater(const double& left, const double& right) {
430  return (isGreater<double, max_ulps>(left, right));
431 }
432 
433 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
434 bool isLessOrEqual(const RawType& left, const RawType& right) {
435  bool is_loe{false};
436 
437  if (not isGreater<RawType, max_ulps>(left, right)) {
438  is_loe = true;
439  }
440 
441  return is_loe;
442 }
443 
444 template<std::size_t max_ulps>
445 inline bool isLessOrEqual(const float& left, const float& right) {
446  return (isLessOrEqual<float, max_ulps>(left, right));
447 }
448 
449 template<std::size_t max_ulps>
450 inline bool isLessOrEqual(const double& left, const double& right) {
451  return (isLessOrEqual<double, max_ulps>(left, right));
452 }
453 
454 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
455 bool isGreaterOrEqual(const RawType& left, const RawType& right) {
456  bool is_goe{false};
457 
458  if (not isLess<RawType, max_ulps>(left, right)) {
459  is_goe = true;
460  }
461 
462  return is_goe;
463 }
464 
465 template<std::size_t max_ulps>
466 inline bool isGreaterOrEqual(const float& left, const float& right) {
467  return (isGreaterOrEqual<float, max_ulps>(left, right));
468 }
469 
470 template<std::size_t max_ulps>
471 inline bool isGreaterOrEqual(const double& left, const double& right) {
472  return (isGreaterOrEqual<double, max_ulps>(left, right));
473 }
474 
491 ELEMENTS_API bool almostEqual2sComplement(const float& left, const float& right,
492  const int& max_ulps = FLT_DEFAULT_MAX_ULPS);
509 ELEMENTS_API bool almostEqual2sComplement(const double& left, const double& right,
510  const int& max_ulps = DBL_DEFAULT_MAX_ULPS);
511 
525 template<typename RawType>
526 ELEMENTS_API bool realBitWiseEqual(const RawType& left, const RawType& right) {
527 #pragma GCC diagnostic push
528 #pragma GCC diagnostic ignored "-Wfloat-equal"
529  return (left == right);
530 #pragma GCC diagnostic pop
531 }
532 
533 } // namespace Elements
534 
535 
536 #endif // ELEMENTSKERNEL_ELEMENTSKERNEL_REAL_H_
537 
Bits exponentBits() const
Definition: Real.h:247
bool isGreaterOrEqual(const RawType &left, const RawType &right)
Definition: Real.h:455
static Bits signAndMagnitudeToBiased(const Bits &sam)
Definition: Real.h:298
Bits fractionBits() const
Definition: Real.h:252
Macro to silence unused variables warnings from the compiler.
static RawType ReinterpretBits(const Bits &bits)
Definition: Real.h:228
ELEMENTS_API const double DBL_DEFAULT_TEST_TOLERANCE
Double precision float default test tolerance.
Definition: Real.cpp:36
constexpr std::size_t FLT_DEFAULT_MAX_ULPS
Single precision float default maximum unit in the last place.
Definition: Real.h:87
bool almostEqual2sComplement(ELEMENTS_UNUSED const FloatType &a, ELEMENTS_UNUSED const FloatType &b, ELEMENTS_UNUSED const std::size_t &max_ulps=0)
Definition: Real.h:330
typename TypeWithSize< sizeof(RawType)>::UInt Bits
Definition: Real.h:177
bool isGreater(const RawType &left, const RawType &right)
Definition: Real.h:413
Bits signBit() const
Definition: Real.h:257
unsigned long long UInt
Definition: Real.h:124
defines the macros to be used for explicit export of the symbols
bool isLessOrEqual(const RawType &left, const RawType &right)
Definition: Real.h:434
bool AlmostEquals(const FloatingPoint &rhs) const
Definition: Real.h:274
constexpr std::size_t defaultMaxUlps< double >()
Definition: Real.h:138
ELEMENTS_API bool realBitWiseEqual(const RawType &left, const RawType &right)
This function compares 2 floating point numbers bitwise. These are the strict equivalent of the &quot;==&quot;...
Definition: Real.h:526
#define ELEMENTS_API
Dummy definitions for the backward compatibility mode.
Definition: Export.h:74
const Bits & bits() const
Definition: Real.h:242
static Bits distanceBetweenSignAndMagnitudeNumbers(const Bits &sam1, const Bits &sam2)
Definition: Real.h:310
bool isNan() const
Definition: Real.h:262
constexpr std::size_t defaultMaxUlps< float >()
Definition: Real.h:133
static RawType Infinity()
Definition: Real.h:235
bool isEqual(const RawType &left, const RawType &right)
Definition: Real.h:351
bool isNan(const RawType &x)
Definition: Real.h:338
FloatingPointUnion m_u
Definition: Real.h:323
bool isNotEqual(const RawType &left, const RawType &right)
Definition: Real.h:376
#define ELEMENTS_UNUSED
Definition: Unused.h:39
constexpr std::size_t defaultMaxUlps()
Definition: Real.h:128
bool isLess(const RawType &left, const RawType &right)
Definition: Real.h:392
unsigned int UInt
Definition: Real.h:116
constexpr std::size_t DBL_DEFAULT_MAX_ULPS
Double precision float default maximum unit in the last place.
Definition: Real.h:89
ELEMENTS_API const double FLT_DEFAULT_TEST_TOLERANCE
Single precision float default test tolerance.
Definition: Real.cpp:35
FloatingPoint(const RawType &x)
Definition: Real.h:219