From e6400b6084c885cbf54c9e507483f92fdf6a32f1 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Wed, 9 Oct 2024 09:55:33 +0200 Subject: [PATCH] refactor: Remove quick math helpers (#3701) This has proven not to be useful as `std` math functions are not really a bottleneck in our current code. I propose to remove them here. If they become useful in the future we can just bring them back from here. --- Core/include/Acts/Propagator/EigenStepper.ipp | 9 +- Core/include/Acts/Utilities/QuickMath.hpp | 90 ------------------ Core/src/Propagator/SympyStepper.cpp | 9 +- Tests/Benchmarks/CMakeLists.txt | 1 - Tests/Benchmarks/QuickMathBenchmark.cpp | 92 ------------------- Tests/UnitTests/Core/Utilities/CMakeLists.txt | 1 - .../Core/Utilities/QuickMathTests.cpp | 64 ------------- 7 files changed, 2 insertions(+), 264 deletions(-) delete mode 100644 Core/include/Acts/Utilities/QuickMath.hpp delete mode 100644 Tests/Benchmarks/QuickMathBenchmark.cpp delete mode 100644 Tests/UnitTests/Core/Utilities/QuickMathTests.cpp diff --git a/Core/include/Acts/Propagator/EigenStepper.ipp b/Core/include/Acts/Propagator/EigenStepper.ipp index cd08613f722..84b305ae431 100644 --- a/Core/include/Acts/Propagator/EigenStepper.ipp +++ b/Core/include/Acts/Propagator/EigenStepper.ipp @@ -10,7 +10,6 @@ #include "Acts/Propagator/ConstrainedStep.hpp" #include "Acts/Propagator/EigenStepperError.hpp" #include "Acts/Propagator/detail/CovarianceEngine.hpp" -#include "Acts/Utilities/QuickMath.hpp" #include @@ -177,17 +176,11 @@ Acts::Result Acts::EigenStepper::step( // This is given by the order of the Runge-Kutta method constexpr double exponent = 0.25; - // Whether to use fast power function if available - constexpr bool tryUseFastPow{false}; - double x = state.options.stepping.stepTolerance / errorEstimate_; - if constexpr (exponent == 0.25 && !tryUseFastPow) { + if constexpr (exponent == 0.25) { // This is 3x faster than std::pow x = std::sqrt(std::sqrt(x)); - } else if constexpr (std::numeric_limits::is_iec559 && - tryUseFastPow) { - x = fastPow(x, exponent); } else { x = std::pow(x, exponent); } diff --git a/Core/include/Acts/Utilities/QuickMath.hpp b/Core/include/Acts/Utilities/QuickMath.hpp deleted file mode 100644 index 2f03ebb5e50..00000000000 --- a/Core/include/Acts/Utilities/QuickMath.hpp +++ /dev/null @@ -1,90 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 CERN for the benefit of the ACTS project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at https://mozilla.org/MPL/2.0/. - -#pragma once - -#include -#include -#include -#include - -namespace Acts { - -/// @brief Fast inverse square root function -/// Taken from https://en.wikipedia.org/wiki/Fast_inverse_square_root -/// @param x the number to calculate the inverse square root of -/// @param iterations the number of newton iterations to perform -template -constexpr T fastInverseSqrt(T x, int iterations = 1) noexcept { - static_assert(std::numeric_limits::is_iec559 && - "Only IEEE 754 is supported"); - static_assert(std::is_same_v || std::is_same_v, - "Only float and double are supported"); - using I = std::conditional_t, std::uint32_t, - std::uint64_t>; - - constexpr I magic = - std::is_same_v ? 0x5f3759df : 0x5fe6eb50c7b537a9; - - union { - T f; - I i; - } u = {x}; - - u.i = magic - (u.i >> 1); - - for (int i = 0; i < iterations; ++i) { - u.f *= 1.5f - 0.5f * x * u.f * u.f; - } - - return u.f; -} - -/// @brief Fast power function @see `std::pow` -/// Taken from -/// https://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp -/// @param a the base -/// @param b the exponent -constexpr double fastPow(double a, double b) { - // enable only on IEEE 754 - static_assert(std::numeric_limits::is_iec559); - - constexpr std::int64_t magic = 0x3FEF127F00000000; - - union { - double f; - std::int64_t i; - } u = {a}; - - u.i = static_cast(b * (u.i - magic) + magic); - - return u.f; -} - -/// @brief Fast power function more precise than @see `fastPow` -/// Taken from -/// https://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp -/// @param a the base -/// @param b the exponent -constexpr double fastPowMorePrecise(double a, double b) { - // binary exponentiation - double r = 1.0; - int exp = std::abs(static_cast(b)); - double base = b > 0 ? a : 1 / a; - while (exp != 0) { - if ((exp & 1) != 0) { - r *= base; - } - base *= base; - exp >>= 1; - } - - return r * fastPow(a, b - static_cast(b)); -} - -} // namespace Acts diff --git a/Core/src/Propagator/SympyStepper.cpp b/Core/src/Propagator/SympyStepper.cpp index 167f44dd7f5..a33e1f567bb 100644 --- a/Core/src/Propagator/SympyStepper.cpp +++ b/Core/src/Propagator/SympyStepper.cpp @@ -10,7 +10,6 @@ #include "Acts/Propagator/detail/SympyCovarianceEngine.hpp" #include "Acts/Propagator/detail/SympyJacobianEngine.hpp" -#include "Acts/Utilities/QuickMath.hpp" #include #include @@ -127,17 +126,11 @@ Result SympyStepper::stepImpl( // This is given by the order of the Runge-Kutta method constexpr double exponent = 0.25; - // Whether to use fast power function if available - constexpr bool tryUseFastPow{false}; - double x = stepTolerance / errorEstimate_; - if constexpr (exponent == 0.25 && !tryUseFastPow) { + if constexpr (exponent == 0.25) { // This is 3x faster than std::pow x = std::sqrt(std::sqrt(x)); - } else if constexpr (std::numeric_limits::is_iec559 && - tryUseFastPow) { - x = fastPow(x, exponent); } else { x = std::pow(x, exponent); } diff --git a/Tests/Benchmarks/CMakeLists.txt b/Tests/Benchmarks/CMakeLists.txt index 8906375caea..b19e67fed59 100644 --- a/Tests/Benchmarks/CMakeLists.txt +++ b/Tests/Benchmarks/CMakeLists.txt @@ -29,7 +29,6 @@ add_benchmark(SurfaceIntersection SurfaceIntersectionBenchmark.cpp) add_benchmark(RayFrustum RayFrustumBenchmark.cpp) add_benchmark(AnnulusBounds AnnulusBoundsBenchmark.cpp) add_benchmark(StraightLineStepper StraightLineStepperBenchmark.cpp) -add_benchmark(QuickMath QuickMathBenchmark.cpp) add_benchmark(SympyStepper SympyStepperBenchmark.cpp) add_benchmark(Stepper StepperBenchmark.cpp) add_benchmark(SourceLink SourceLinkBenchmark.cpp) diff --git a/Tests/Benchmarks/QuickMathBenchmark.cpp b/Tests/Benchmarks/QuickMathBenchmark.cpp deleted file mode 100644 index 0cbec176b94..00000000000 --- a/Tests/Benchmarks/QuickMathBenchmark.cpp +++ /dev/null @@ -1,92 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 CERN for the benefit of the ACTS project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at https://mozilla.org/MPL/2.0/. - -#include -#include - -#include "Acts/Tests/CommonHelpers/BenchmarkTools.hpp" -#include "Acts/Utilities/QuickMath.hpp" - -#include -#include - -namespace bdata = boost::unit_test::data; - -namespace Acts::Test { - -/// @brief Another fast power function @see `fastPow` -// Taken from -// https://martin.ankerl.com/2007/02/11/optimized-exponential-functions-for-java -/// @param a the base -/// @param b the exponent -constexpr double fastPowAnother(double a, double b) { - // enable only on IEEE 754 - static_assert(std::numeric_limits::is_iec559); - - union { - double f; - std::int64_t i; - } u = {}; - - u.i = static_cast( - 9076650 * (a - 1) / (a + 1 + 4 * std::sqrt(a)) * b + 1072632447); - u.i <<= 32; - - // result seems broken? - return u.f; -} - -// Some randomness & number crunching -const unsigned int nTests = 10; -const unsigned int nReps = 10000; - -BOOST_DATA_TEST_CASE( - benchmark_pow_25, - bdata::random( - (bdata::engine = std::mt19937(), bdata::seed = 21, - bdata::distribution = std::uniform_real_distribution(-4, 4))) ^ - bdata::xrange(nTests), - baseExp, index) { - (void)index; - - const double base = std::pow(10, baseExp); - const double exp = 0.25; - - std::cout << std::endl - << "Benchmarking base=" << base << ", exp=" << exp << "..." - << std::endl; - std::cout << "- void: " - << Acts::Test::microBenchmark([&] { return 0; }, nReps) - << std::endl; - std::cout << "- std::pow: " - << Acts::Test::microBenchmark([&] { return std::pow(base, exp); }, - nReps) - << std::endl; - std::cout << "- std::exp: " - << Acts::Test::microBenchmark( - [&] { return std::exp(std::log(base) * exp); }, nReps) - << std::endl; - std::cout << "- std::sqrt: " - << Acts::Test::microBenchmark( - [&] { return std::sqrt(std::sqrt(base)); }, nReps) - << std::endl; - std::cout << "- fastPow: " - << Acts::Test::microBenchmark([&] { return fastPow(base, exp); }, - nReps) - << std::endl; - std::cout << "- fastPowMorePrecise: " - << Acts::Test::microBenchmark( - [&] { return fastPowMorePrecise(base, exp); }, nReps) - << std::endl; - std::cout << "- fastPowAnother: " - << Acts::Test::microBenchmark( - [&] { return fastPowAnother(base, exp); }, nReps) - << std::endl; -} - -} // namespace Acts::Test diff --git a/Tests/UnitTests/Core/Utilities/CMakeLists.txt b/Tests/UnitTests/Core/Utilities/CMakeLists.txt index 68b4ba5a655..27b91bfbe29 100644 --- a/Tests/UnitTests/Core/Utilities/CMakeLists.txt +++ b/Tests/UnitTests/Core/Utilities/CMakeLists.txt @@ -54,7 +54,6 @@ target_compile_definitions( add_unittest(ParticleData ParticleDataTests.cpp) add_unittest(Zip ZipTests.cpp) add_unittest(TransformRange TransformRangeTests.cpp) -add_unittest(QuickMath QuickMathTests.cpp) add_unittest(TrackHelpers TrackHelpersTests.cpp) add_unittest(GraphViz GraphVizTests.cpp) diff --git a/Tests/UnitTests/Core/Utilities/QuickMathTests.cpp b/Tests/UnitTests/Core/Utilities/QuickMathTests.cpp deleted file mode 100644 index 6952ea78b14..00000000000 --- a/Tests/UnitTests/Core/Utilities/QuickMathTests.cpp +++ /dev/null @@ -1,64 +0,0 @@ -// This file is part of the ACTS project. -// -// Copyright (C) 2016 CERN for the benefit of the ACTS project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at https://mozilla.org/MPL/2.0/. - -#include -#include - -#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" -#include "Acts/Utilities/QuickMath.hpp" - -namespace bdata = boost::unit_test::data; - -const auto expDist = bdata::random( - (bdata::engine = std::mt19937{}, bdata::seed = 0, - bdata::distribution = std::uniform_real_distribution(-4, 4))); - -BOOST_AUTO_TEST_SUITE(Utilities) - -BOOST_DATA_TEST_CASE(fastInverseSqrt, expDist ^ bdata::xrange(100), exp, i) { - (void)i; - - const double x = std::pow(10, exp); - - const float fastFloat = Acts::fastInverseSqrt(static_cast(x)); - const double fastDouble = Acts::fastInverseSqrt(x); - - const double stdFloat = 1.0 / std::sqrt(static_cast(x)); - const double stdDouble = 1.0 / std::sqrt(x); - - CHECK_CLOSE_REL(stdFloat, fastFloat, 0.01); - CHECK_CLOSE_REL(stdDouble, fastDouble, 0.01); -} - -BOOST_DATA_TEST_CASE(fastPow, expDist ^ expDist ^ bdata::xrange(100), baseExp, - exp, i) { - (void)i; - - const double base = std::pow(10, baseExp); - - const double fast = Acts::fastPow(base, exp); - const double fastMorePrecise = Acts::fastPowMorePrecise(base, exp); - - const double std = std::pow(base, exp); - - CHECK_CLOSE_REL(fast, std, 0.15); - CHECK_CLOSE_REL(fastMorePrecise, std, 0.1); -} - -// BOOST_AUTO_TEST_CASE(fastPowChart) { -// std::cout << "a ref obs" << std::endl; -// for (double aExp = -4; aExp <= 4; aExp += 0.01) { -// double a = std::pow(10, aExp); -// double ref = std::pow(a, 0.25); -// double obs = Acts::fastPow(a, 0.25); - -// std::cout << a << " " << ref << " " << obs << std::endl; -// } -// } - -BOOST_AUTO_TEST_SUITE_END()