FloatComparison.h

// Copyright © 2008-2023 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 <limits>
#include <cstdint>
#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