Elements  6.0.1
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 <cmath> // for round
76 #include <limits> // for numeric_limits
77 #include <type_traits> // for is_floating_point
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 template <std::size_t size>
100 public:
101  // This prevents the user from using TypeWithSize<N> with incorrect
102  // values of N.
103  using UInt = void;
104 };
105 
106 // The specialisation for size 4.
107 template <>
109 public:
110  // unsigned int has size 4 in both gcc and MSVC.
111  //
112  // As base/basictypes.h doesn't compile on Windows, we cannot use
113  // uint32, uint64, and etc here.
114  using Int = int;
115  using UInt = unsigned int;
116 };
117 
118 // The specialisation for size 8.
119 template <>
121 public:
122  using Int = long long; // NOLINT
123  using UInt = unsigned long long; // NOLINT
124 };
125 
126 template <typename RawType>
128  return FLT_DEFAULT_MAX_ULPS;
129 }
130 
131 template <>
133  return FLT_DEFAULT_MAX_ULPS;
134 }
135 
136 template <>
138  return DBL_DEFAULT_MAX_ULPS;
139 }
140 
141 // This template class represents an IEEE floating-point number
142 // (either single-precision or double-precision, depending on the
143 // template parameters).
144 //
145 // The purpose of this class is to do more sophisticated number
146 // comparison. (Due to round-off error, etc, it's very unlikely that
147 // two floating-points will be equal exactly. Hence a naive
148 // comparison by the == operation often doesn't work.)
149 //
150 // Format of IEEE floating-point:
151 //
152 // The most-significant bit being the leftmost, an IEEE
153 // floating-point looks like
154 //
155 // sign_bit exponent_bits fraction_bits
156 //
157 // Here, sign_bit is a single bit that designates the sign of the
158 // number.
159 //
160 // For float, there are 8 exponent bits and 23 fraction bits.
161 //
162 // For double, there are 11 exponent bits and 52 fraction bits.
163 //
164 // More details can be found at
165 // http://en.wikipedia.org/wiki/IEEE_floating-point_standard.
166 //
167 // Template parameter:
168 //
169 // RawType: the raw floating-point type (either float or double)
170 template <typename RawType>
172 public:
173  // Defines the unsigned integer type that has the same size as the
174  // floating point number.
176 
177  // Constants.
178 
179  // # of bits in a number.
180  static const std::size_t s_bitcount = 8 * sizeof(RawType);
181 
182  // # of fraction bits in a number.
183  static const std::size_t s_fraction_bitcount = std::numeric_limits<RawType>::digits - 1;
184 
185  // # of exponent bits in a number.
186  static const std::size_t s_exponent_bitcount = s_bitcount - 1 - s_fraction_bitcount;
187 
188  // The mask for the sign bit.
189  static const Bits s_sign_bitmask = static_cast<Bits>(1) << (s_bitcount - 1);
190 
191  // The mask for the fraction bits.
192  static const Bits s_fraction_bitmask = ~static_cast<Bits>(0) >> (s_exponent_bitcount + 1);
193 
194  // The mask for the exponent bits.
195  static const Bits s_exponent_bitmask = ~(s_sign_bitmask | s_fraction_bitmask);
196 
197  // How many ULP's (Units in the Last Place) we want to tolerate when
198  // comparing two numbers. The larger the value, the more error we
199  // allow. A 0 value means that two numbers must be exactly the same
200  // to be considered equal.
201  //
202  // The maximum error of a single floating-point operation is 0.5
203  // units in the last place. On Intel CPU's, all floating-point
204  // calculations are done with 80-bit precision, while double has 64
205  // bits. Therefore, 4 should be enough for ordinary use.
206  //
207  // See the following article for more details on ULP:
208  // http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm.
209  static const std::size_t m_max_ulps = defaultMaxUlps<RawType>();
210 
211  // Constructs a FloatingPoint from a raw floating-point number.
212  //
213  // On an Intel CPU, passing a non-normalised NAN (Not a Number)
214  // around may change its bits, although the new value is guaranteed
215  // to be also a NAN. Therefore, don't expect this constructor to
216  // preserve the bits in x when x is a NAN.
217  explicit FloatingPoint(const RawType& x) {
218  m_u.m_value = x;
219  }
220 
221  // Static methods
222 
223  // Reinterprets a bit pattern as a floating-point number.
224  //
225  // This function is needed to test the AlmostEquals() method.
226  static RawType ReinterpretBits(const Bits& bits) {
227  FloatingPoint fp(0);
228  fp.m_u.m_bits = bits;
229  return fp.m_u.m_value;
230  }
231 
232  // Returns the floating-point number that represent positive infinity.
233  static RawType Infinity() {
234  return ReinterpretBits(s_exponent_bitmask);
235  }
236 
237  // Non-static methods
238 
239  // Returns the bits that represents this number.
240  const Bits& bits() const {
241  return m_u.m_bits;
242  }
243 
244  // Returns the exponent bits of this number.
245  Bits exponentBits() const {
246  return s_exponent_bitmask & m_u.m_bits;
247  }
248 
249  // Returns the fraction bits of this number.
250  Bits fractionBits() const {
251  return s_fraction_bitmask & m_u.m_bits;
252  }
253 
254  // Returns the sign bit of this number.
255  Bits signBit() const {
256  return s_sign_bitmask & m_u.m_bits;
257  }
258 
259  // Returns true iff this is NAN (not a number).
260  bool isNan() const {
261  // It's a NAN if the exponent bits are all ones and the fraction
262  // bits are not entirely zeros.
263  return (exponentBits() == s_exponent_bitmask) && (fractionBits() != 0);
264  }
265 
266  // Returns true iff this number is at most kMaxUlps ULP's away from
267  // rhs. In particular, this function:
268  //
269  // - returns false if either number is (or both are) NAN.
270  // - treats really large numbers as almost equal to infinity.
271  // - thinks +0.0 and -0.0 are 0 DLP's apart.
272  bool AlmostEquals(const FloatingPoint& rhs) const {
273  // The IEEE standard says that any comparison operation involving
274  // a NAN must return false.
275  if (isNan() || rhs.isNan()) {
276  return false;
277  }
278  return distanceBetweenSignAndMagnitudeNumbers(m_u.m_bits, rhs.m_u.m_bits) <= m_max_ulps;
279  }
280 
281  // Converts an integer from the sign-and-magnitude representation to
282  // the biased representation. More precisely, let N be 2 to the
283  // power of (kBitCount - 1), an integer x is represented by the
284  // unsigned number x + N.
285  //
286  // For instance,
287  //
288  // -N + 1 (the most negative number representable using
289  // sign-and-magnitude) is represented by 1;
290  // 0 is represented by N; and
291  // N - 1 (the biggest number representable using
292  // sign-and-magnitude) is represented by 2N - 1.
293  //
294  // Read http://en.wikipedia.org/wiki/Signed_number_representations
295  // for more details on signed number representations.
296  static Bits signAndMagnitudeToBiased(const Bits& sam) {
297  if (s_sign_bitmask & sam) {
298  // sam represents a negative number.
299  return ~sam + 1;
300  } else {
301  // sam represents a positive number.
302  return s_sign_bitmask | sam;
303  }
304  }
305 
306  // Given two numbers in the sign-and-magnitude representation,
307  // returns the distance between them as an unsigned number.
308  static Bits distanceBetweenSignAndMagnitudeNumbers(const Bits& sam1, const Bits& sam2) {
309  const Bits biased1 = signAndMagnitudeToBiased(sam1);
310  const Bits biased2 = signAndMagnitudeToBiased(sam2);
311  return (biased1 >= biased2) ? (biased1 - biased2) : (biased2 - biased1);
312  }
313 
314 private:
315  // The data type used to store the actual floating-point number.
317  RawType m_value; // The raw floating-point number.
318  Bits m_bits; // The bits that represent the number.
319  };
320 
322 };
323 
324 // Usable AlmostEqual function
325 
326 template <typename FloatType>
327 bool almostEqual2sComplement(ELEMENTS_UNUSED const FloatType& a, ELEMENTS_UNUSED const FloatType& b,
328  ELEMENTS_UNUSED const std::size_t& max_ulps = 0) {
329  return false;
330 }
331 
332 template <typename RawType>
333 bool isNan(const RawType& x) {
334 
335  using Bits = typename TypeWithSize<sizeof(RawType)>::UInt;
336  Bits x_bits = *reinterpret_cast<const Bits*>(&x);
337 
338  Bits x_exp_bits = FloatingPoint<RawType>::s_exponent_bitmask & x_bits;
339  Bits x_frac_bits = FloatingPoint<RawType>::s_fraction_bitmask & x_bits;
340 
341  return (x_exp_bits == FloatingPoint<RawType>::s_exponent_bitmask) && (x_frac_bits != 0);
342 }
343 
344 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
345 bool isEqual(const RawType& left, const RawType& right) {
346 
347  bool is_equal{false};
348 
349  if (not(isNan<RawType>(left) or isNan<RawType>(right))) {
350  using Bits = typename TypeWithSize<sizeof(RawType)>::UInt;
351  Bits l_bits = *reinterpret_cast<const Bits*>(&left);
352  Bits r_bits = *reinterpret_cast<const Bits*>(&right);
353  is_equal = (FloatingPoint<RawType>::distanceBetweenSignAndMagnitudeNumbers(l_bits, r_bits) <= max_ulps);
354  }
355 
356  return is_equal;
357 }
358 
359 template <std::size_t max_ulps>
360 inline bool isEqual(const float& left, const float& right) {
361  return (isEqual<float, max_ulps>(left, right));
362 }
363 
364 template <std::size_t max_ulps>
365 inline bool isEqual(const double& left, const double& right) {
366  return (isEqual<double, max_ulps>(left, right));
367 }
368 
369 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
370 inline bool isNotEqual(const RawType& left, const RawType& right) {
371  return (not isEqual<RawType, max_ulps>(left, right));
372 }
373 
374 template <std::size_t max_ulps>
375 inline bool isNotEqual(const float& left, const float& right) {
376  return (isNotEqual<float, max_ulps>(left, right));
377 }
378 
379 template <std::size_t max_ulps>
380 inline bool isNotEqual(const double& left, const double& right) {
381  return (isNotEqual<double, max_ulps>(left, right));
382 }
383 
384 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
385 bool isLess(const RawType& left, const RawType& right) {
386  bool is_less{false};
387 
388  if (left < right && (not isEqual<RawType, max_ulps>(left, right))) {
389  is_less = true;
390  }
391 
392  return is_less;
393 }
394 
395 template <std::size_t max_ulps>
396 inline bool isLess(const float& left, const float& right) {
397  return (isLess<float, max_ulps>(left, right));
398 }
399 
400 template <std::size_t max_ulps>
401 inline bool isLess(const double& left, const double& right) {
402  return (isLess<double, max_ulps>(left, right));
403 }
404 
405 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
406 bool isGreater(const RawType& left, const RawType& right) {
407  bool is_greater{false};
408 
409  if (left > right && (not isEqual<RawType, max_ulps>(left, right))) {
410  is_greater = true;
411  }
412 
413  return is_greater;
414 }
415 
416 template <std::size_t max_ulps>
417 inline bool isGreater(const float& left, const float& right) {
418  return (isGreater<float, max_ulps>(left, right));
419 }
420 
421 template <std::size_t max_ulps>
422 inline bool isGreater(const double& left, const double& right) {
423  return (isGreater<double, max_ulps>(left, right));
424 }
425 
426 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
427 bool isLessOrEqual(const RawType& left, const RawType& right) {
428  bool is_loe{false};
429 
430  if (not isGreater<RawType, max_ulps>(left, right)) {
431  is_loe = true;
432  }
433 
434  return is_loe;
435 }
436 
437 template <std::size_t max_ulps>
438 inline bool isLessOrEqual(const float& left, const float& right) {
439  return (isLessOrEqual<float, max_ulps>(left, right));
440 }
441 
442 template <std::size_t max_ulps>
443 inline bool isLessOrEqual(const double& left, const double& right) {
444  return (isLessOrEqual<double, max_ulps>(left, right));
445 }
446 
447 template <typename RawType, std::size_t max_ulps = defaultMaxUlps<RawType>()>
448 bool isGreaterOrEqual(const RawType& left, const RawType& right) {
449  bool is_goe{false};
450 
451  if (not isLess<RawType, max_ulps>(left, right)) {
452  is_goe = true;
453  }
454 
455  return is_goe;
456 }
457 
458 template <std::size_t max_ulps>
459 inline bool isGreaterOrEqual(const float& left, const float& right) {
460  return (isGreaterOrEqual<float, max_ulps>(left, right));
461 }
462 
463 template <std::size_t max_ulps>
464 inline bool isGreaterOrEqual(const double& left, const double& right) {
465  return (isGreaterOrEqual<double, max_ulps>(left, right));
466 }
467 
484 ELEMENTS_API bool almostEqual2sComplement(const float& left, const float& right,
485  const int& max_ulps = FLT_DEFAULT_MAX_ULPS);
502 ELEMENTS_API bool almostEqual2sComplement(const double& left, const double& right,
503  const int& max_ulps = DBL_DEFAULT_MAX_ULPS);
504 
518 template <typename RawType>
519 ELEMENTS_API bool realBitWiseEqual(const RawType& left, const RawType& right) {
520 #pragma GCC diagnostic push
521 #pragma GCC diagnostic ignored "-Wfloat-equal"
522  return (left == right);
523 #pragma GCC diagnostic pop
524 }
525 
526 } // namespace Elements
527 
528 #endif // ELEMENTSKERNEL_ELEMENTSKERNEL_REAL_H_
529 
Bits exponentBits() const
Definition: Real.h:245
bool isGreaterOrEqual(const RawType &left, const RawType &right)
Definition: Real.h:448
static Bits signAndMagnitudeToBiased(const Bits &sam)
Definition: Real.h:296
Bits fractionBits() const
Definition: Real.h:250
Macro to silence unused variables warnings from the compiler.
static RawType ReinterpretBits(const Bits &bits)
Definition: Real.h:226
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:327
typename TypeWithSize< sizeof(RawType)>::UInt Bits
Definition: Real.h:175
bool isGreater(const RawType &left, const RawType &right)
Definition: Real.h:406
Bits signBit() const
Definition: Real.h:255
unsigned long long UInt
Definition: Real.h:123
defines the macros to be used for explicit export of the symbols
bool isLessOrEqual(const RawType &left, const RawType &right)
Definition: Real.h:427
bool AlmostEquals(const FloatingPoint &rhs) const
Definition: Real.h:272
constexpr std::size_t defaultMaxUlps< double >()
Definition: Real.h:137
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:519
#define ELEMENTS_API
Dummy definitions for the backward compatibility mode.
Definition: Export.h:74
const Bits & bits() const
Definition: Real.h:240
static Bits distanceBetweenSignAndMagnitudeNumbers(const Bits &sam1, const Bits &sam2)
Definition: Real.h:308
bool isNan() const
Definition: Real.h:260
constexpr std::size_t defaultMaxUlps< float >()
Definition: Real.h:132
static RawType Infinity()
Definition: Real.h:233
bool isEqual(const RawType &left, const RawType &right)
Definition: Real.h:345
bool isNan(const RawType &x)
Definition: Real.h:333
FloatingPointUnion m_u
Definition: Real.h:321
bool isNotEqual(const RawType &left, const RawType &right)
Definition: Real.h:370
#define ELEMENTS_UNUSED
Definition: Unused.h:39
constexpr std::size_t defaultMaxUlps()
Definition: Real.h:127
bool isLess(const RawType &left, const RawType &right)
Definition: Real.h:385
unsigned int UInt
Definition: Real.h:115
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:217