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].
This commit is contained in:
Martin Hořeňovský 2023-11-19 21:00:53 +01:00
parent bfd9f0f5a6
commit d139b4ff7c
No known key found for this signature in database
GPG Key ID: DE48307B8B0D381A
8 changed files with 175 additions and 16 deletions

View File

@ -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

View File

@ -10,6 +10,7 @@
#include <catch2/benchmark/detail/catch_stats.hpp>
#include <catch2/internal/catch_compiler_capabilities.hpp>
#include <catch2/internal/catch_floating_point_helpers.hpp>
#include <algorithm>
#include <cassert>
@ -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<double>( n );
// degenerate case with uniform samples
if ( directCompare( prob_n, 0. ) ) {
if ( Catch::Detail::directCompare( prob_n, 0. ) ) {
return { point, point, point, confidence_level };
}

View File

@ -91,6 +91,7 @@
#include <catch2/internal/catch_preprocessor.hpp>
#include <catch2/internal/catch_preprocessor_internal_stringify.hpp>
#include <catch2/internal/catch_preprocessor_remove_parens.hpp>
#include <catch2/internal/catch_random_floating_point_helpers.hpp>
#include <catch2/internal/catch_random_number_generator.hpp>
#include <catch2/internal/catch_random_seed_generation.hpp>
#include <catch2/internal/catch_reporter_registry.hpp>

View File

@ -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

View File

@ -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

View File

@ -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 <catch2/internal/catch_polyfills.hpp>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <limits>
#include <type_traits>
namespace Catch {
namespace Detail {
/**
* Returns the largest magnitude of 1-ULP distance inside the [a, b] range.
*
* Assumes `a < b`.
*/
template <typename FloatType>
FloatType gamma(FloatType a, FloatType b) {
static_assert( std::is_floating_point<FloatType>::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<FloatType>::infinity() ) - a;
const auto gamma_down = b - Catch::nextafter( b, -std::numeric_limits<FloatType>::infinity() );
return gamma_up < gamma_down ? gamma_down : gamma_up;
}
template <typename FloatingPoint>
struct DistanceTypePicker;
template <>
struct DistanceTypePicker<float> {
using type = std::uint32_t;
};
template <>
struct DistanceTypePicker<double> {
using type = std::uint64_t;
};
template <typename T>
using DistanceType = typename DistanceTypePicker<T>::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 <typename FloatType>
DistanceType<FloatType>
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<DistanceType<FloatType>>( 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

View File

@ -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',

View File

@ -9,7 +9,9 @@
#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_template_test_macros.hpp>
#include <catch2/internal/catch_floating_point_helpers.hpp>
#include <catch2/internal/catch_random_floating_point_helpers.hpp>
#include <limits>
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<TestType>::epsilon() );
CHECK( gamma( TestType( -1. ), TestType( -0. ) ) <
std::numeric_limits<TestType>::epsilon() );
CHECK( directCompare( gamma( TestType( 1. ), TestType( 2. ) ),
std::numeric_limits<TestType>::epsilon() ) );
CHECK( directCompare( gamma( TestType( -2. ), TestType( -1. ) ),
std::numeric_limits<TestType>::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<std::uint64_t,
decltype( count_floats_with_scaled_ulp(
0., 1. ) )>::value );
STATIC_REQUIRE( std::is_same<std::uint32_t,
decltype( count_floats_with_scaled_ulp(
0.f, 1.f ) )>::value );
}