261 lines
8.0 KiB
C++
261 lines
8.0 KiB
C++
// Copyright © 2008-2021 Pioneer Developers. See AUTHORS.txt for details
|
|
// Licensed under the terms of the GPL v3. See licenses/GPL-3.txt
|
|
|
|
#ifndef _FLOATCOMPARISON_H
|
|
#define _FLOATCOMPARISON_H
|
|
|
|
#include <SDL_stdinc.h>
|
|
#include <limits>
|
|
#ifdef _MSC_VER
|
|
#include <float.h> // for _finite
|
|
#else
|
|
#include <cmath> // for std::isfinite
|
|
#endif
|
|
|
|
// Fuzzy floating point comparisons based on:
|
|
// http://realtimecollisiondetection.net/blog/?p=89
|
|
// (absolute & relative error tolerance)
|
|
//
|
|
// http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
|
|
// (ULP based tolerance)
|
|
//
|
|
// ULP-based tolerance implementation takes some architectural ideas from the
|
|
// implementation in the Google test framework, and
|
|
// http://stackoverflow.com/questions/17333/most-effective-way-for-float-and-double-comparison
|
|
|
|
// provides (for float & double):
|
|
|
|
// bool is_equal_exact(float a, float b);
|
|
// bool is_equal_ulps(float a, float b, int ulps = DefaultUlpTolerance);
|
|
// int32_t float_ulp_difference(float a, float b);
|
|
// bool is_equal_relative(float a, float b, float tolerance = DefaultRelTolerance());
|
|
// bool is_equal_absolute(float a, float b, float tolerance = DefaultAbsTolerance());
|
|
// bool is_equal_general(float a, float b, float tolerance = DefaultTolerance());
|
|
// bool is_equal_general(float a, float b, float relative_tolerance, float absolute_tolerance);
|
|
|
|
// bool is_zero_exact(float x);
|
|
// bool is_zero_general(float x, float tolerance = IEEEFloatTraits<float>::DefaultRelTolerance());
|
|
|
|
// bool is_nan(float x);
|
|
// bool is_finite(float x);
|
|
|
|
// ====================================================================
|
|
|
|
// in the following code, IEEEFloatTraits<T>::bool_type is used to limit
|
|
// the application of the functions by SFINAE
|
|
|
|
template <typename T>
|
|
struct IEEEFloatTraits;
|
|
|
|
// --- float function helpers
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::float_type float_abs(T x)
|
|
{
|
|
return (x < T(0)) ? (-x) : x;
|
|
}
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::float_type float_max(T x, T y)
|
|
{
|
|
return (y > x) ? y : x;
|
|
}
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::float_type float_max(T x, T y, T z)
|
|
{
|
|
return float_max(x, float_max(y, z));
|
|
}
|
|
|
|
// --- float property helpers
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::bool_type is_nan_bits(const typename IEEEFloatTraits<T>::uint_type &bits)
|
|
{
|
|
typedef typename IEEEFloatTraits<T>::uint_type uint_type;
|
|
const uint_type top = IEEEFloatTraits<T>::TopBit;
|
|
const uint_type ebits = IEEEFloatTraits<T>::ExponentBits;
|
|
|
|
// NaN has the exponent bits set, and at least one mantissa bit set
|
|
// (therefore, if you mask off the top bit, the result must be strictly greater than
|
|
// just the exponent bits set; if it's equal then it's just an infinity; if it's
|
|
// less, then it's a valid finite number)
|
|
return ((bits & ~top) > ebits);
|
|
}
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::bool_type is_finite_bits(const typename IEEEFloatTraits<T>::uint_type &bits)
|
|
{
|
|
typedef typename IEEEFloatTraits<T>::uint_type uint_type;
|
|
const uint_type ebits = IEEEFloatTraits<T>::ExponentBits;
|
|
return ((bits & ebits) != ebits);
|
|
}
|
|
|
|
// --- infinity
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::bool_type is_finite(T x)
|
|
{
|
|
#ifdef _MSC_VER
|
|
return _finite(x);
|
|
#else
|
|
return std::isfinite(x);
|
|
#endif
|
|
}
|
|
|
|
// --- exact comparisons, and checking for NaN
|
|
|
|
#ifdef __GNUC__
|
|
#pragma GCC diagnostic push
|
|
#pragma GCC diagnostic ignored "-Wfloat-equal"
|
|
#endif
|
|
inline bool is_equal_exact(float a, float b)
|
|
{
|
|
return (a == b);
|
|
}
|
|
inline bool is_equal_exact(double a, double b) { return (a == b); }
|
|
|
|
inline bool is_zero_exact(float x) { return (x == 0.0f); }
|
|
inline bool is_zero_exact(double x) { return (x == 0.0); }
|
|
|
|
inline bool is_nan(float x) { return (x != x); }
|
|
inline bool is_nan(double x) { return (x != x); }
|
|
#ifdef __GNUC__
|
|
#pragma GCC diagnostic pop
|
|
#endif
|
|
|
|
// --- relative & absolute error comparisons
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::bool_type is_equal_relative(T a, T b, T tol = IEEEFloatTraits<T>::DefaultRelTolerance())
|
|
{
|
|
return (float_abs(a - b) <= tol * float_max(float_abs(a), float_abs(b)));
|
|
}
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::bool_type is_equal_absolute(T a, T b, T tol = IEEEFloatTraits<T>::DefaultAbsTolerance())
|
|
{
|
|
return (float_abs(a - b) <= tol);
|
|
}
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::bool_type is_equal_general(T a, T b, T rel_tol, T abs_tol)
|
|
{
|
|
return (float_abs(a - b) <= float_max(abs_tol, rel_tol * float_max(float_abs(a), float_abs(b))));
|
|
}
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::bool_type is_equal_general(T a, T b, T tol = IEEEFloatTraits<T>::DefaultTolerance())
|
|
{
|
|
return (float_abs(a - b) <= tol * float_max(T(1), float_abs(a), float_abs(b)));
|
|
}
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::bool_type is_zero_general(T x, T tol = IEEEFloatTraits<T>::DefaultRelTolerance())
|
|
{
|
|
return (float_abs(x) <= tol);
|
|
}
|
|
|
|
// --- ulp-based comparisons
|
|
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::int_type float_ulp_difference(T a, T b)
|
|
{
|
|
typedef typename IEEEFloatTraits<T>::FloatOrInt union_type;
|
|
union_type afi, bfi;
|
|
afi.f = a;
|
|
bfi.f = b;
|
|
|
|
// transform from sign-magnitude to two's-complement
|
|
if (afi.i < 0) afi.ui = (IEEEFloatTraits<T>::TopBit - afi.ui);
|
|
if (bfi.i < 0) bfi.ui = (IEEEFloatTraits<T>::TopBit - bfi.ui);
|
|
|
|
return (bfi.i - afi.i);
|
|
}
|
|
|
|
// IEEEFloatTraits<T>::bool_type used for SFINAE
|
|
template <typename T>
|
|
inline typename IEEEFloatTraits<T>::bool_type is_equal_ulps(T a, T b, typename IEEEFloatTraits<T>::int_type max_ulps = IEEEFloatTraits<T>::DefaultUlpTolerance)
|
|
{
|
|
typedef typename IEEEFloatTraits<T>::FloatOrInt union_type;
|
|
typedef typename IEEEFloatTraits<T>::int_type int_type;
|
|
union_type afi, bfi;
|
|
afi.f = a;
|
|
bfi.f = b;
|
|
|
|
// Infinities aren't close to anything except themselves
|
|
if ((!is_finite_bits<T>(afi.ui) && is_finite_bits<T>(bfi.ui)) || (is_finite_bits<T>(afi.ui) && !is_finite_bits<T>(bfi.ui)))
|
|
return false;
|
|
|
|
// IEEE says NaNs are unequal to everything (even themselves)
|
|
if (is_nan_bits<T>(afi.ui) || is_nan_bits<T>(bfi.ui))
|
|
return false;
|
|
|
|
// transform from sign-magnitude to two's-complement
|
|
if (afi.i < 0) afi.ui = (IEEEFloatTraits<T>::TopBit - afi.ui);
|
|
if (bfi.i < 0) bfi.ui = (IEEEFloatTraits<T>::TopBit - bfi.ui);
|
|
|
|
int_type difference = (bfi.i - afi.i);
|
|
difference = (difference < int_type(0)) ? -difference : difference;
|
|
return (difference <= max_ulps);
|
|
}
|
|
|
|
// ====================================================================
|
|
|
|
template <typename T>
|
|
struct IEEEFloatTraits {};
|
|
|
|
template <>
|
|
struct IEEEFloatTraits<double> {
|
|
typedef double float_type;
|
|
typedef bool bool_type;
|
|
|
|
typedef int64_t int_type;
|
|
typedef uint64_t uint_type;
|
|
|
|
union FloatOrInt {
|
|
double f;
|
|
uint_type ui;
|
|
int_type i;
|
|
};
|
|
|
|
static const uint_type TopBit = static_cast<uint_type>(1) << (sizeof(double) * 8 - 1);
|
|
static const uint_type ExponentBits = (~static_cast<uint_type>(0) << std::numeric_limits<double>::digits) & ~TopBit;
|
|
static const uint_type MantissaBits = ~TopBit & ~ExponentBits;
|
|
|
|
static const int_type DefaultUlpTolerance = 16;
|
|
|
|
static double DefaultAbsTolerance() { return 1e-12; }
|
|
static double DefaultRelTolerance() { return 1e-6; }
|
|
static double DefaultTolerance() { return 1e-8; }
|
|
static double SmallestNormalisedValue() { return std::numeric_limits<double>::min(); }
|
|
};
|
|
|
|
template <>
|
|
struct IEEEFloatTraits<float> {
|
|
typedef float float_type;
|
|
typedef bool bool_type;
|
|
|
|
typedef int32_t int_type;
|
|
typedef uint32_t uint_type;
|
|
|
|
union FloatOrInt {
|
|
float f;
|
|
uint_type ui;
|
|
int_type i;
|
|
};
|
|
|
|
static const uint_type TopBit = uint_type(1) << (sizeof(float) * 8 - 1);
|
|
static const uint_type ExponentBits = (~uint_type(0) << std::numeric_limits<float>::digits) & ~TopBit;
|
|
static const uint_type MantissaBits = ~TopBit & ~ExponentBits;
|
|
|
|
static const int_type DefaultUlpTolerance = 4;
|
|
|
|
static float DefaultAbsTolerance() { return 1e-6f; }
|
|
static float DefaultRelTolerance() { return 1e-5f; }
|
|
static float DefaultTolerance() { return 1e-5f; }
|
|
static float SmallestNormalisedValue() { return std::numeric_limits<float>::min(); }
|
|
};
|
|
|
|
#endif
|