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