From d139b4ff7cd8823923e8616e8fef954e09e46322 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Martin=20Ho=C5=99e=C5=88ovsk=C3=BD?= Date: Sun, 19 Nov 2023 21:00:53 +0100 Subject: [PATCH] Add implementation of helpers for uniform float distribution Specifically we add * `gamma(a, b)`, which returns the magnitude of largest 1-ULP step in range [a, b]. * `count_equidistant_float(a, b, distance)`, which returns the number of equi-distant floats in range [a, b]. --- src/CMakeLists.txt | 1 + src/catch2/benchmark/detail/catch_stats.cpp | 19 +--- src/catch2/catch_all.hpp | 1 + .../internal/catch_floating_point_helpers.cpp | 11 +++ .../internal/catch_floating_point_helpers.hpp | 5 + .../catch_random_floating_point_helpers.hpp | 94 +++++++++++++++++++ src/catch2/meson.build | 1 + .../FloatingPoint.tests.cpp | 59 ++++++++++++ 8 files changed, 175 insertions(+), 16 deletions(-) create mode 100644 src/catch2/internal/catch_random_floating_point_helpers.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a70eef03..01f05b4c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -109,6 +109,7 @@ set(IMPL_HEADERS ${SOURCES_DIR}/internal/catch_polyfills.hpp ${SOURCES_DIR}/internal/catch_preprocessor.hpp ${SOURCES_DIR}/internal/catch_preprocessor_remove_parens.hpp + ${SOURCES_DIR}/internal/catch_random_floating_point_helpers.hpp ${SOURCES_DIR}/internal/catch_random_number_generator.hpp ${SOURCES_DIR}/internal/catch_random_seed_generation.hpp ${SOURCES_DIR}/internal/catch_reporter_registry.hpp diff --git a/src/catch2/benchmark/detail/catch_stats.cpp b/src/catch2/benchmark/detail/catch_stats.cpp index ac255e36..6030e55c 100644 --- a/src/catch2/benchmark/detail/catch_stats.cpp +++ b/src/catch2/benchmark/detail/catch_stats.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include @@ -184,20 +185,6 @@ namespace Catch { return std::sqrt( variance ); } -#if defined( __GNUC__ ) || defined( __clang__ ) -# pragma GCC diagnostic push -# pragma GCC diagnostic ignored "-Wfloat-equal" -#endif - // Used when we know we want == comparison of two doubles - // to centralize warning suppression - static bool directCompare( double lhs, double rhs ) { - return lhs == rhs; - } -#if defined( __GNUC__ ) || defined( __clang__ ) -# pragma GCC diagnostic pop -#endif - - static sample jackknife( double ( *estimator )( double const*, double const* ), double* first, @@ -234,7 +221,7 @@ namespace Catch { double g = idx - j; std::nth_element(first, first + j, last); auto xj = first[j]; - if ( directCompare( g, 0 ) ) { + if ( Catch::Detail::directCompare( g, 0 ) ) { return xj; } @@ -338,7 +325,7 @@ namespace Catch { [point]( double x ) { return x < point; } ) / static_cast( n ); // degenerate case with uniform samples - if ( directCompare( prob_n, 0. ) ) { + if ( Catch::Detail::directCompare( prob_n, 0. ) ) { return { point, point, point, confidence_level }; } diff --git a/src/catch2/catch_all.hpp b/src/catch2/catch_all.hpp index a7969fef..3621c505 100644 --- a/src/catch2/catch_all.hpp +++ b/src/catch2/catch_all.hpp @@ -91,6 +91,7 @@ #include #include #include +#include #include #include #include diff --git a/src/catch2/internal/catch_floating_point_helpers.cpp b/src/catch2/internal/catch_floating_point_helpers.cpp index e30ee434..9631ed6d 100644 --- a/src/catch2/internal/catch_floating_point_helpers.cpp +++ b/src/catch2/internal/catch_floating_point_helpers.cpp @@ -27,6 +27,17 @@ namespace Catch { return i; } +#if defined( __GNUC__ ) || defined( __clang__ ) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wfloat-equal" +#endif + bool directCompare( float lhs, float rhs ) { return lhs == rhs; } + bool directCompare( double lhs, double rhs ) { return lhs == rhs; } +#if defined( __GNUC__ ) || defined( __clang__ ) +# pragma GCC diagnostic pop +#endif + + } // end namespace Detail } // end namespace Catch diff --git a/src/catch2/internal/catch_floating_point_helpers.hpp b/src/catch2/internal/catch_floating_point_helpers.hpp index ca883c61..b2143726 100644 --- a/src/catch2/internal/catch_floating_point_helpers.hpp +++ b/src/catch2/internal/catch_floating_point_helpers.hpp @@ -22,6 +22,11 @@ namespace Catch { uint32_t convertToBits(float f); uint64_t convertToBits(double d); + // Used when we know we want == comparison of two doubles + // to centralize warning suppression + bool directCompare( float lhs, float rhs ); + bool directCompare( double lhs, double rhs ); + } // end namespace Detail diff --git a/src/catch2/internal/catch_random_floating_point_helpers.hpp b/src/catch2/internal/catch_random_floating_point_helpers.hpp new file mode 100644 index 00000000..d75b4218 --- /dev/null +++ b/src/catch2/internal/catch_random_floating_point_helpers.hpp @@ -0,0 +1,94 @@ + +// Copyright Catch2 Authors +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE.txt or copy at +// https://www.boost.org/LICENSE_1_0.txt) + +// SPDX-License-Identifier: BSL-1.0 + +#ifndef CATCH_RANDOM_FLOATING_POINT_HELPERS_HPP_INCLUDED +#define CATCH_RANDOM_FLOATING_POINT_HELPERS_HPP_INCLUDED + +#include + +#include +#include +#include +#include +#include + +namespace Catch { + + namespace Detail { + /** + * Returns the largest magnitude of 1-ULP distance inside the [a, b] range. + * + * Assumes `a < b`. + */ + template + FloatType gamma(FloatType a, FloatType b) { + static_assert( std::is_floating_point::value, + "gamma returns the largest ULP magnitude within " + "floating point range [a, b]. This only makes sense " + "for floating point types" ); + assert( a < b ); + + const auto gamma_up = Catch::nextafter( a, std::numeric_limits::infinity() ) - a; + const auto gamma_down = b - Catch::nextafter( b, -std::numeric_limits::infinity() ); + + return gamma_up < gamma_down ? gamma_down : gamma_up; + } + + template + struct DistanceTypePicker; + template <> + struct DistanceTypePicker { + using type = std::uint32_t; + }; + template <> + struct DistanceTypePicker { + using type = std::uint64_t; + }; + + template + using DistanceType = typename DistanceTypePicker::type; + +#if defined( __GNUC__ ) || defined( __clang__ ) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wfloat-equal" +#endif + /** + * Computes the number of equi-distant floats in [a, b] + * + * Since not every range can be split into equidistant floats + * exactly, we actually compute ceil(b/distance - a/distance), + * because in those cases we want to overcount. + * + * Uses modified Dekker's FastTwoSum algorithm to handle rounding. + */ + template + DistanceType + count_equidistant_floats( FloatType a, FloatType b, FloatType distance ) { + assert( a < b ); + // We get distance as gamma for our uniform float distribution, + // so this will round perfectly. + const auto ag = a / distance; + const auto bg = b / distance; + + const auto s = bg - ag; + const auto err = ( std::fabs( a ) <= std::fabs( b ) ) + ? -ag - ( s - bg ) + : bg - ( s + ag ); + const auto ceil_s = static_cast>( std::ceil( s ) ); + + return ( ceil_s != s ) ? ceil_s : ceil_s + ( err > 0 ); + } +#if defined( __GNUC__ ) || defined( __clang__ ) +# pragma GCC diagnostic pop +#endif + + } + +} // end namespace Catch + +#endif // CATCH_RANDOM_FLOATING_POINT_HELPERS_HPP_INCLUDED diff --git a/src/catch2/meson.build b/src/catch2/meson.build index 09a3db1b..fee7a622 100644 --- a/src/catch2/meson.build +++ b/src/catch2/meson.build @@ -115,6 +115,7 @@ internal_headers = [ 'internal/catch_preprocessor.hpp', 'internal/catch_preprocessor_internal_stringify.hpp', 'internal/catch_preprocessor_remove_parens.hpp', + 'internal/catch_random_floating_point_helpers.hpp', 'internal/catch_random_number_generator.hpp', 'internal/catch_random_seed_generation.hpp', 'internal/catch_reporter_registry.hpp', diff --git a/tests/SelfTest/IntrospectiveTests/FloatingPoint.tests.cpp b/tests/SelfTest/IntrospectiveTests/FloatingPoint.tests.cpp index 08a579c9..eaf22a44 100644 --- a/tests/SelfTest/IntrospectiveTests/FloatingPoint.tests.cpp +++ b/tests/SelfTest/IntrospectiveTests/FloatingPoint.tests.cpp @@ -9,7 +9,9 @@ #include #include #include +#include +#include TEST_CASE("convertToBits", "[floating-point][conversion]") { using Catch::Detail::convertToBits; @@ -72,3 +74,60 @@ TEST_CASE("UlpDistance", "[floating-point][ulp][approvals]") { CHECK( ulpDistance( 1.f, 2.f ) == 0x80'00'00 ); CHECK( ulpDistance( -2.f, 2.f ) == 0x80'00'00'00 ); } + + + +TEMPLATE_TEST_CASE("gamma", "[approvals][floating-point][ulp][gamma]", float, double) { + using Catch::Detail::gamma; + using Catch::Detail::directCompare; + + // We need to butcher the equal tests with the directCompare helper, + // because the Wfloat-equal triggers in decomposer rather than here, + // so we cannot locally disable it. Goddamn GCC. + CHECK( directCompare( gamma( TestType( -1. ), TestType( 1. ) ), + gamma( TestType( 0.2332 ), TestType( 1.0 ) ) ) ); + CHECK( directCompare( gamma( TestType( -2. ), TestType( 0 ) ), + gamma( TestType( 1. ), TestType( 1.5 ) ) ) ); + CHECK( gamma( TestType( 0. ), TestType( 1.0 ) ) < + gamma( TestType( 1.0 ), TestType( 1.5 ) ) ); + CHECK( gamma( TestType( 0 ), TestType( 1. ) ) < + std::numeric_limits::epsilon() ); + CHECK( gamma( TestType( -1. ), TestType( -0. ) ) < + std::numeric_limits::epsilon() ); + CHECK( directCompare( gamma( TestType( 1. ), TestType( 2. ) ), + std::numeric_limits::epsilon() ) ); + CHECK( directCompare( gamma( TestType( -2. ), TestType( -1. ) ), + std::numeric_limits::epsilon() ) ); +} + +TEMPLATE_TEST_CASE("count_equidistant_floats", + "[approvals][floating-point][distance]", + float, + double) { + using Catch::Detail::count_equidistant_floats; + auto count_steps = []( TestType a, TestType b ) { + return count_equidistant_floats( a, b, Catch::Detail::gamma( a, b ) ); + }; + + CHECK( count_steps( TestType( -1. ), TestType( 1. ) ) == + 2 * count_steps( TestType( 0. ), TestType( 1. ) ) ); +} + +TEST_CASE( "count_equidistant_floats", + "[approvals][floating-point][distance]" ) { + using Catch::Detail::count_equidistant_floats; + auto count_floats_with_scaled_ulp = []( auto a, auto b ) { + return count_equidistant_floats( a, b, Catch::Detail::gamma( a, b ) ); + }; + + CHECK( count_floats_with_scaled_ulp( 1., 1.5 ) == 1ull << 51 ); + CHECK( count_floats_with_scaled_ulp( 1.25, 1.5 ) == 1ull << 50 ); + CHECK( count_floats_with_scaled_ulp( 1.f, 1.5f ) == 1 << 22 ); + + STATIC_REQUIRE( std::is_same::value ); + STATIC_REQUIRE( std::is_same::value ); +}