From d5993fb05ea7678464d18896889b402d5e51c3ed Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 28 Sep 2024 12:41:39 +0200 Subject: [PATCH 1/4] =?UTF-8?q?Compute=20=CE=B5=20and=20M=CA=B9=20using=20?= =?UTF-8?q?only=20double.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- functions/accurate_table_generator.hpp | 10 +- functions/accurate_table_generator_body.hpp | 22 ++-- functions/accurate_table_generator_test.cpp | 122 +++++++++----------- 3 files changed, 72 insertions(+), 82 deletions(-) diff --git a/functions/accurate_table_generator.hpp b/functions/accurate_table_generator.hpp index 48d203a0e4..2f49be052f 100644 --- a/functions/accurate_table_generator.hpp +++ b/functions/accurate_table_generator.hpp @@ -19,6 +19,7 @@ using namespace principia::base::_thread_pool; using namespace principia::numerics::_polynomial_in_monomial_basis; using AccurateFunction = std::function; +using ApproximateFunction = std::function; template using AccuratePolynomial = @@ -41,7 +42,7 @@ template absl::StatusOr StehléZimmermannSimultaneousSearch( std::array const& functions, std::array, 2> const& polynomials, - std::array const& remainders, + std::array const& remainders, cpp_rational const& starting_argument, std::int64_t N, std::int64_t T); @@ -54,7 +55,7 @@ template absl::StatusOr StehléZimmermannSimultaneousFullSearch( std::array const& functions, std::array, 2> const& polynomials, - std::array const& remainders, + std::array const& remainders, cpp_rational const& starting_argument, ThreadPool* search_pool = nullptr); @@ -67,7 +68,7 @@ StehléZimmermannSimultaneousMultisearch( std::array const& functions, std::vector, 2>> const& polynomials, - std::vector> const& remainders, + std::vector> const& remainders, std::vector const& starting_arguments); // Same as above, but instead of accumulating all the results and returning them @@ -78,7 +79,7 @@ void StehléZimmermannSimultaneousStreamingMultisearch( std::array const& functions, std::vector, 2>> const& polynomials, - std::vector> const& remainders, + std::vector> const& remainders, std::vector const& starting_arguments, std::function)> const& callback); @@ -87,6 +88,7 @@ void StehléZimmermannSimultaneousStreamingMultisearch( using internal::AccurateFunction; using internal::AccuratePolynomial; +using internal::ApproximateFunction; using internal::GalExhaustiveMultisearch; using internal::GalExhaustiveSearch; using internal::StehléZimmermannSimultaneousFullSearch; diff --git a/functions/accurate_table_generator_body.hpp b/functions/accurate_table_generator_body.hpp index 09b806d965..35240aa677 100644 --- a/functions/accurate_table_generator_body.hpp +++ b/functions/accurate_table_generator_body.hpp @@ -68,7 +68,7 @@ bool AllFunctionValuesHaveDesiredZeroes( struct StehléZimmermannSpecification { std::array functions; std::array, 2> polynomials; - std::array remainders; + std::array remainders; cpp_rational argument; }; @@ -136,7 +136,7 @@ StehléZimmermannSpecification ScaleToBinade01( std::array function_scales; std::array scaled_functions; - std::array scaled_remainders; + std::array scaled_remainders; for (std::int64_t i = 0; i < scaled_functions.size(); ++i) { function_scales[i] = compute_scale(functions[i](lower_bound), functions[i](upper_bound)); @@ -439,7 +439,7 @@ template absl::StatusOr StehléZimmermannSimultaneousSearch( std::array const& functions, std::array, 2> const& polynomials, - std::array const& remainders, + std::array const& remainders, cpp_rational const& starting_argument, std::int64_t const N, std::int64_t const T) { @@ -473,17 +473,17 @@ absl::StatusOr StehléZimmermannSimultaneousSearch( // Step 2: compute ε. We use the remainders provided by the clients. Note // that we could save on the number of evaluations by providing both bounds to // a single call. - cpp_bin_float_50 ε = 0; + double ε = 0; + auto const T_over_N = cpp_rational(T, N); for (std::int64_t i = 0; i < remainders.size(); ++i) { - auto const T_over_N = cpp_rational(T, N); - ε = std::max(ε, abs(N * remainders[i](starting_argument - T_over_N))); - ε = std::max(ε, abs(N * remainders[i](starting_argument + T_over_N))); + ε = std::max(ε, std::abs(N * remainders[i](starting_argument - T_over_N))); + ε = std::max(ε, std::abs(N * remainders[i](starting_argument + T_over_N))); } VLOG(3) << "ε = " << ε; // Step 3, first part: compute Mʹ and C. Give up is C is 0, which may happen // if ε is too large. - auto const Mʹ = static_cast(floor(M / (2 + 2 * M * ε))); + auto const Mʹ = static_cast(std::floor(M / (2 + 2 * M * ε))); auto const C = 3 * Mʹ; if (C == 0) { return absl::FailedPreconditionError("Error too large"); @@ -616,7 +616,7 @@ template absl::StatusOr StehléZimmermannSimultaneousFullSearch( std::array const& functions, std::array, 2> const& polynomials, - std::array const& remainders, + std::array const& remainders, cpp_rational const& starting_argument, ThreadPool* const search_pool) { // Start by scaling the specification of the search. The rest of this @@ -772,7 +772,7 @@ StehléZimmermannSimultaneousMultisearch( std::array const& functions, std::vector, 2>> const& polynomials, - std::vector> const& remainders, + std::vector> const& remainders, std::vector const& starting_arguments) { std::vector> result; result.resize(starting_arguments.size()); @@ -793,7 +793,7 @@ void StehléZimmermannSimultaneousStreamingMultisearch( std::array const& functions, std::vector, 2>> const& polynomials, - std::vector> const& remainders, + std::vector> const& remainders, std::vector const& starting_arguments, std::function)> const& callback) { diff --git a/functions/accurate_table_generator_test.cpp b/functions/accurate_table_generator_test.cpp index 626911f3f6..b8fd8e1b2e 100644 --- a/functions/accurate_table_generator_test.cpp +++ b/functions/accurate_table_generator_test.cpp @@ -145,17 +145,15 @@ TEST_F(AccurateTableGeneratorTest, StehléZimmermannSinCos15) { // Use the Lagrange form of the remainder. u may be above or below u₀. Note // that we use the fact that the functions are monotonic. - auto const remainder_sin_taylor2 = - [u₀ = cpp_rational(u₀)](cpp_rational const& u) { - auto const Δu = u - u₀; - auto const Δu³ = static_cast(Δu * Δu * Δu); - return Δu³ * (-Cos(std::min(u, u₀) / 4) / 16) / Factorial(3); + auto const remainder_sin_taylor2 = [u₀](cpp_rational const& u) { + auto const Δu = static_cast(u) - u₀; + auto const Δu³ = Δu * Δu * Δu; + return Δu³ * (-std::cos(std::min(u₀ + Δu, u₀) / 4) / 16) / Factorial(3); }; - auto const remainder_cos_taylor2 = - [u₀ = cpp_rational(u₀)](cpp_rational const& u) { - auto const Δu = u - u₀; - auto const Δu³ = static_cast(Δu * Δu * Δu); - return Δu³ * (Sin(std::max(u, u₀) / 4) / 64) / Factorial(3); + auto const remainder_cos_taylor2 = [u₀](cpp_rational const& u) { + auto const Δu = static_cast(u) - u₀; + auto const Δu³ = Δu * Δu * Δu; + return Δu³ * (std::sin(std::max(u₀ + Δu, u₀) / 4) / 64) / Factorial(3); }; auto const u = StehléZimmermannSimultaneousSearch<15>( @@ -207,17 +205,15 @@ TEST_F(AccurateTableGeneratorTest, StehléZimmermannFullSinCos5NoScaling) { -cpp_rational(Cos(u₀ / 4) / 32)}, u₀); - auto const remainder_sin_taylor2 = - [u₀ = cpp_rational(u₀)](cpp_rational const& u) { - auto const Δu = u - u₀; - auto const Δu³ = static_cast(Δu * Δu * Δu); - return Δu³ * (-Cos(std::min(u, u₀) / 4) / 16) / Factorial(3); + auto const remainder_sin_taylor2 = [u₀](cpp_rational const& u) { + auto const Δu = static_cast(u) - u₀; + auto const Δu³ = Δu * Δu * Δu; + return Δu³ * (-std::cos(std::min(u₀ + Δu, u₀) / 4) / 16) / Factorial(3); }; - auto const remainder_cos_taylor2 = - [u₀ = cpp_rational(u₀)](cpp_rational const& u) { - auto const Δu = u - u₀; - auto const Δu³ = static_cast(Δu * Δu * Δu); - return Δu³ * (Sin(std::max(u, u₀) / 4) / 64) / Factorial(3); + auto const remainder_cos_taylor2 = [u₀](cpp_rational const& u) { + auto const Δu = static_cast(u) - u₀; + auto const Δu³ = Δu * Δu * Δu; + return Δu³ * (std::sin(std::max(u₀ + Δu, u₀) / 4) / 64) / Factorial(3); }; auto const u = StehléZimmermannSimultaneousFullSearch<5>( @@ -267,17 +263,15 @@ TEST_F(AccurateTableGeneratorTest, StehléZimmermannFullSinCos15NoScaling) { -cpp_rational(Cos(u₀ / 4) / 32)}, u₀); - auto const remainder_sin_taylor2 = - [u₀ = cpp_rational(u₀)](cpp_rational const& u) { - auto const Δu = u - u₀; - auto const Δu³ = static_cast(Δu * Δu * Δu); - return Δu³ * (-Cos(std::min(u, u₀) / 4) / 16) / Factorial(3); + auto const remainder_sin_taylor2 = [u₀](cpp_rational const& u) { + auto const Δu = static_cast(u) - u₀; + auto const Δu³ = Δu * Δu * Δu; + return Δu³ * (-std::cos(std::min(u₀ + Δu, u₀) / 4) / 16) / Factorial(3); }; - auto const remainder_cos_taylor2 = - [u₀ = cpp_rational(u₀)](cpp_rational const& u) { - auto const Δu = u - u₀; - auto const Δu³ = static_cast(Δu * Δu * Δu); - return Δu³ * (Sin(std::max(u, u₀) / 4) / 64) / Factorial(3); + auto const remainder_cos_taylor2 = [u₀](cpp_rational const& u) { + auto const Δu = static_cast(u) - u₀; + auto const Δu³ = Δu * Δu * Δu; + return Δu³ * (std::sin(std::max(u₀ + Δu, u₀) / 4) / 64) / Factorial(3); }; auto const u = StehléZimmermannSimultaneousFullSearch<15>( @@ -319,17 +313,15 @@ TEST_F(AccurateTableGeneratorTest, StehléZimmermannFullSinCos15WithScaling) { -cpp_rational(Cos(x₀) / 2)}, x₀); - auto const remainder_sin_taylor2 = - [x₀ = cpp_rational(x₀)](cpp_rational const& x) { - auto const Δx = x - x₀; - auto const Δx³ = static_cast(Δx * Δx * Δx); - return -Δx³ * -Cos(std::min(x, x₀)) / Factorial(3); + auto const remainder_sin_taylor2 = [x₀](cpp_rational const& x) { + auto const Δx = static_cast(x) - x₀; + auto const Δx³ = Δx * Δx * Δx; + return -Δx³ * -std::cos(std::min(x₀ + Δx, x₀)) / Factorial(3); }; - auto const remainder_cos_taylor2 = - [x₀ = cpp_rational(x₀)](cpp_rational const& x) { - auto const Δx = x - x₀; - auto const Δx³ = static_cast(Δx * Δx * Δx); - return Δx³ * Sin(std::max(x, x₀)) / Factorial(3); + auto const remainder_cos_taylor2 = [x₀](cpp_rational const& x) { + auto const Δx = static_cast(x) - x₀; + auto const Δx³ = Δx * Δx * Δx; + return Δx³ * std::sin(std::max(x₀ + Δx, x₀)) / Factorial(3); }; auto const x = StehléZimmermannSimultaneousFullSearch<15>( @@ -365,7 +357,7 @@ TEST_F(AccurateTableGeneratorTest, StehléZimmermannMultisearchSinCos15) { static constexpr std::int64_t index_end = 100; std::vector starting_arguments; std::vector, 2>> polynomials; - std::vector> remainders; + std::vector> remainders; for (std::int64_t i = index_begin; i < index_end; ++i) { auto const x₀ = i / 128.0; AccuratePolynomial const sin_taylor2( @@ -379,18 +371,16 @@ TEST_F(AccurateTableGeneratorTest, StehléZimmermannMultisearchSinCos15) { -cpp_rational(Cos(x₀) / 2)}, x₀); - auto const remainder_sin_taylor2 = - [x₀ = cpp_rational(x₀)](cpp_rational const& x) { - auto const Δx = x - x₀; - auto const Δx³ = static_cast(Δx * Δx * Δx); - return -Δx³ * Cos(std::min(x, x₀)) / Factorial(3); - }; - auto const remainder_cos_taylor2 = - [x₀ = cpp_rational(x₀)](cpp_rational const& x) { - auto const Δx = x - x₀; - auto const Δx³ = static_cast(Δx * Δx * Δx); - return Δx³ * Sin(std::max(x, x₀)) / Factorial(3); - }; + auto const remainder_sin_taylor2 = [x₀](cpp_rational const& x) { + auto const Δx = static_cast(x) - x₀; + auto const Δx³ = Δx * Δx * Δx; + return -Δx³ * -std::cos(std::min(x₀ + Δx, x₀)) / Factorial(3); + }; + auto const remainder_cos_taylor2 = [x₀](cpp_rational const& x) { + auto const Δx = static_cast(x) - x₀; + auto const Δx³ = Δx * Δx * Δx; + return Δx³ * std::sin(std::max(x₀ + Δx, x₀)) / Factorial(3); + }; starting_arguments.push_back(x₀); polynomials.push_back({sin_taylor2, cos_taylor2}); @@ -452,8 +442,8 @@ TEST_F(AccurateTableGeneratorTest, DISABLED_SECULAR_SinCos18) { std::vector starting_arguments; std::vector, 2>> polynomials; - std::vector> remainders; - for (std::int64_t i = i_min; i <= i_max; ++i) { + std::vector> remainders; + for (std::int64_t i = 1; i <= 1; ++i) { double const x₀ = centre(i); AccuratePolynomial const sin_taylor2( {cpp_rational(Sin(x₀)), @@ -468,18 +458,16 @@ TEST_F(AccurateTableGeneratorTest, DISABLED_SECULAR_SinCos18) { // The remainders don't need to be extremely precise, so for speed // they are computed using double. - auto const remainder_sin_taylor2 = - [x₀ = static_cast(cpp_rational(x₀))](cpp_rational const& x) { - auto const Δx = static_cast(x) - x₀; - auto const Δx³ = Δx * Δx * Δx; - return -Δx³ * -std::cos(std::min(x₀ + Δx, x₀)) / Factorial(3); - }; - auto const remainder_cos_taylor2 = - [x₀ = static_cast(cpp_rational(x₀))](cpp_rational const& x) { - auto const Δx = static_cast(x) - x₀; - auto const Δx³ = Δx * Δx * Δx; - return Δx³ * std::sin(std::max(x₀ + Δx, x₀)) / Factorial(3); - }; + auto const remainder_sin_taylor2 = [x₀](cpp_rational const& x) { + auto const Δx = static_cast(x) - x₀; + auto const Δx³ = Δx * Δx * Δx; + return -Δx³ * -std::cos(std::min(x₀ + Δx, x₀)) / Factorial(3); + }; + auto const remainder_cos_taylor2 = [x₀](cpp_rational const& x) { + auto const Δx = static_cast(x) - x₀; + auto const Δx³ = Δx * Δx * Δx; + return Δx³ * std::sin(std::max(x₀ + Δx, x₀)) / Factorial(3); + }; starting_arguments.push_back(x₀); polynomials.push_back({sin_taylor2, cos_taylor2}); From 0c89ba7eede6640ac402c018ceea49b8630688be Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 28 Sep 2024 13:02:03 +0200 Subject: [PATCH 2/4] Use the basic boost sin/cos functions. --- functions/accurate_table_generator_test.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/functions/accurate_table_generator_test.cpp b/functions/accurate_table_generator_test.cpp index b8fd8e1b2e..000f573fe3 100644 --- a/functions/accurate_table_generator_test.cpp +++ b/functions/accurate_table_generator_test.cpp @@ -440,6 +440,14 @@ TEST_F(AccurateTableGeneratorTest, DISABLED_SECULAR_SinCos18) { CHECK_LT(centre(i_max) - h, π / 4); CHECK_LT(π / 4, centre(i_max) + h); + // No need for fancy angle reduction as the angles are small. + AccurateFunction const accurate_sin = [](cpp_rational const& x) { + return sin(static_cast(x)); + }; + AccurateFunction const accurate_cos = [](cpp_rational const& x) { + return cos(static_cast(x)); + }; + std::vector starting_arguments; std::vector, 2>> polynomials; std::vector> remainders; @@ -475,7 +483,7 @@ TEST_F(AccurateTableGeneratorTest, DISABLED_SECULAR_SinCos18) { } StehléZimmermannSimultaneousStreamingMultisearch<18>( - {Sin, Cos}, + {accurate_sin, accurate_cos}, polynomials, remainders, starting_arguments, From 2352e4f531f598395fd850df49e47e6da713aac5 Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 28 Sep 2024 13:17:01 +0200 Subject: [PATCH 3/4] Tune T_max. --- functions/accurate_table_generator_body.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/functions/accurate_table_generator_body.hpp b/functions/accurate_table_generator_body.hpp index 35240aa677..9201c5cfa3 100644 --- a/functions/accurate_table_generator_body.hpp +++ b/functions/accurate_table_generator_body.hpp @@ -39,7 +39,7 @@ using namespace principia::numerics::_matrix_views; using namespace principia::quantities::_elementary_functions; using namespace principia::quantities::_quantities; -constexpr std::int64_t T_max = 16; +constexpr std::int64_t T_max = 4; static_assert(T_max >= 1); template From d5c5562a75cdc8555ab79d59b1b7b5e5bda6a826 Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 28 Sep 2024 13:22:02 +0200 Subject: [PATCH 4/4] Revert testing changes. --- functions/accurate_table_generator_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/functions/accurate_table_generator_test.cpp b/functions/accurate_table_generator_test.cpp index 000f573fe3..e1996cdf9c 100644 --- a/functions/accurate_table_generator_test.cpp +++ b/functions/accurate_table_generator_test.cpp @@ -451,7 +451,7 @@ TEST_F(AccurateTableGeneratorTest, DISABLED_SECULAR_SinCos18) { std::vector starting_arguments; std::vector, 2>> polynomials; std::vector> remainders; - for (std::int64_t i = 1; i <= 1; ++i) { + for (std::int64_t i = i_min; i <= i_max; ++i) { double const x₀ = centre(i); AccuratePolynomial const sin_taylor2( {cpp_rational(Sin(x₀)),