From 78e5ea385e884d1f0b17705ce4f5a5f37b45a525 Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 20 Jan 2024 12:25:21 +0100 Subject: [PATCH 1/3] Explicit control of max_degree. --- numerics/approximation.hpp | 5 +++- numerics/approximation_body.hpp | 44 ++++++++++++++++++++++----------- numerics/approximation_test.cpp | 16 ++++++------ physics/apsides_body.hpp | 11 +++++---- 4 files changed, 47 insertions(+), 29 deletions(-) diff --git a/numerics/approximation.hpp b/numerics/approximation.hpp index 9b8762405c..d54e4fd8ff 100644 --- a/numerics/approximation.hpp +++ b/numerics/approximation.hpp @@ -16,7 +16,10 @@ using namespace principia::quantities::_named_quantities; template using Value = std::invoke_result_t; -template +// Returns a Чебышёв polynomial approximant of f over +// [lower_bound, upper_bound]. Stops if the absolute error is estimated to be +// below |max_error| or if |max_degree| has been reached. +template ЧебышёвSeries, Argument> ЧебышёвPolynomialInterpolant( Function const& f, Argument const& lower_bound, diff --git a/numerics/approximation_body.hpp b/numerics/approximation_body.hpp index bac7a2babd..2bd85bb362 100644 --- a/numerics/approximation_body.hpp +++ b/numerics/approximation_body.hpp @@ -5,6 +5,7 @@ #include #include +#include "base/tags.hpp" #include "geometry/barycentre_calculator.hpp" #include "numerics/fixed_arrays.hpp" #include "quantities/elementary_functions.hpp" @@ -15,14 +16,30 @@ namespace numerics { namespace _approximation { namespace internal { +using namespace principia::base::_tags; using namespace principia::geometry::_barycentre_calculator; using namespace principia::numerics::_fixed_arrays; using namespace principia::quantities::_elementary_functions; using namespace principia::quantities::_si; -constexpr std::int64_t max_чебышёв_degree = 500; +// Compute the interpolation matrix and cache it in a static variable. +template +FixedMatrix const& ЧебышёвInterpolationMatrix() { + static FixedMatrix const ℐ = []() { + FixedMatrix ℐ(uninitialized); + for (std::int64_t j = 0; j <= N; ++j) { + double const pⱼ = j == 0 || j == N ? 2 : 1; + for (std::int64_t k = 0; k <= N; ++k) { + double const pₖ = k == 0 || k == N ? 2 : 1; + ℐ(j, k) = 2 * Cos(π * j * k * Radian / N) / (pⱼ * pₖ * N); + } + } + return ℐ; + }(); + return ℐ; +} -template +template ЧебышёвSeries, Argument> ЧебышёвPolynomialInterpolantImplementation( Function const& f, @@ -39,7 +56,7 @@ template return 0.5 * (b - a) * Cos(π * k * Radian / N) + midpoint; }; - FixedVector, N + 1> fₖ; + FixedVector, N + 1> fₖ(uninitialized); // Reuse the previous evaluations of |f|. for (std::int64_t k = 0; k <= N / 2; ++k) { @@ -51,16 +68,9 @@ template fₖ[k] = f(чебышёв_lobato_point(k)); } - FixedMatrix ℐⱼₖ; - for (std::int64_t j = 0; j <= N; ++j) { - double const pⱼ = j == 0 || j == N ? 2 : 1; - for (std::int64_t k = 0; k <= N; ++k) { - double const pₖ = k == 0 || k == N ? 2 : 1; - ℐⱼₖ(j, k) = 2 * Cos(π * j * k * Radian / N) / (pⱼ * pₖ * N); - } - } - // Compute the coefficients of the Чебышёв polynomial. + FixedMatrix const& ℐⱼₖ = + ЧебышёвInterpolationMatrix(); auto const aⱼ = ℐⱼₖ * fₖ; // Compute an upper bound for the error, based on the previous and new @@ -73,9 +83,12 @@ template error_estimate += Abs(aⱼ[j]); } - if constexpr (2 * N <= max_чебышёв_degree) { + if constexpr (N <= max_degree) { if (error_estimate > max_error) { - return ЧебышёвPolynomialInterpolantImplementation<2 * N>( + // Note that this recursive call overflows the stack when + // max_degree >= 256. We could allocate on the heap, but then we don't + // care about very high degree polynomials. + return ЧебышёвPolynomialInterpolantImplementation<2 * N, max_degree>( f, a, b, max_error, fₖ, aⱼ); } } @@ -93,7 +106,7 @@ template coefficients, a, b); } -template +template ЧебышёвSeries, Argument> ЧебышёвPolynomialInterpolant( Function const& f, @@ -108,6 +121,7 @@ template FixedVector, 2> const aⱼ( {0.5 * (f_b + f_a), 0.5 * (f_b - f_a)}); return ЧебышёвPolynomialInterpolantImplementation( f, a, b, max_error, fₖ, aⱼ); diff --git a/numerics/approximation_test.cpp b/numerics/approximation_test.cpp index d3e82a96b5..6298512152 100644 --- a/numerics/approximation_test.cpp +++ b/numerics/approximation_test.cpp @@ -22,10 +22,10 @@ using ::testing::Lt; TEST(ApproximationTest, SinInverse) { auto const f = [](double const x) { return Sin(1 * Radian / x); }; auto const approximation = - ЧебышёвPolynomialInterpolant(f, - /*lower_bound=*/0.1, - /*upper_bound=*/10, - /*max_error=*/1e-6); + ЧебышёвPolynomialInterpolant<128>(f, + /*lower_bound=*/0.1, + /*upper_bound=*/10.0, + /*max_error=*/1e-6); EXPECT_EQ(128, approximation.degree()); for (double x = 0.1; x < 10.0; x += 0.01) { EXPECT_THAT(approximation.Evaluate(x), @@ -36,10 +36,10 @@ TEST(ApproximationTest, SinInverse) { TEST(ApproximationTest, Exp) { auto const f = [](double const x) { return std::exp(x); }; auto const approximation = - ЧебышёвPolynomialInterpolant(f, - /*lower_bound=*/0.01, - /*upper_bound=*/3, - /*max_error=*/1e-6); + ЧебышёвPolynomialInterpolant<128>(f, + /*lower_bound=*/0.01, + /*upper_bound=*/3.0, + /*max_error=*/1e-6); EXPECT_EQ(16, approximation.degree()); for (double x = 0.01; x < 3; x += 0.01) { EXPECT_THAT(approximation.Evaluate(x), diff --git a/physics/apsides_body.hpp b/physics/apsides_body.hpp index 61a448faaa..c56620113c 100644 --- a/physics/apsides_body.hpp +++ b/physics/apsides_body.hpp @@ -322,11 +322,12 @@ ComputeFirstCollision( }; // Interpolate the height above the terrain using a Чебышёв polynomial. - auto const чебышёв_interpolant = ЧебышёвPolynomialInterpolant( - height_above_terrain_at_time, - interval.min, - interval.max, - reference_body.mean_radius() * max_error_relative_to_radius); + auto const чебышёв_interpolant = + ЧебышёвPolynomialInterpolant( + height_above_terrain_at_time, + interval.min, + interval.max, + reference_body.mean_radius() * max_error_relative_to_radius); auto const& real_roots = чебышёв_interpolant.RealRoots( max_error_relative_to_radius); if (real_roots.empty()) { From 15a75d52488c08632906d29af5df4080e0d0609d Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 20 Jan 2024 12:45:27 +0100 Subject: [PATCH 2/3] Return the error estimate. --- numerics/approximation.hpp | 6 ++++-- numerics/approximation_body.hpp | 22 ++++++++++++++-------- numerics/approximation_test.cpp | 15 +++++++++++++-- 3 files changed, 31 insertions(+), 12 deletions(-) diff --git a/numerics/approximation.hpp b/numerics/approximation.hpp index d54e4fd8ff..9643eeecf8 100644 --- a/numerics/approximation.hpp +++ b/numerics/approximation.hpp @@ -18,13 +18,15 @@ using Value = std::invoke_result_t; // Returns a Чебышёв polynomial approximant of f over // [lower_bound, upper_bound]. Stops if the absolute error is estimated to be -// below |max_error| or if |max_degree| has been reached. +// below |max_error| or if |max_degree| has been reached. If |error_estimate| +// is nonnull, it receives the estimate of the L∞ error. template ЧебышёвSeries, Argument> ЧебышёвPolynomialInterpolant( Function const& f, Argument const& lower_bound, Argument const& upper_bound, - Difference> const& max_error); + Difference> const& max_error, + Difference>* error_estimate = nullptr); } // namespace internal diff --git a/numerics/approximation_body.hpp b/numerics/approximation_body.hpp index 2bd85bb362..a56393f5c5 100644 --- a/numerics/approximation_body.hpp +++ b/numerics/approximation_body.hpp @@ -47,7 +47,8 @@ template Argument const& b, Difference> const& max_error, FixedVector, N / 2 + 1> const& previous_fₖ, - FixedVector, N / 2 + 1> const& previous_aⱼ) { + FixedVector, N / 2 + 1> const& previous_aⱼ, + Difference>* const error_estimate) { // This implementation follows [Boy13], section 4 and appendix A. auto const midpoint = Barycentre(std::pair{a, b}, std::pair{0.5, 0.5}); @@ -75,21 +76,21 @@ template // Compute an upper bound for the error, based on the previous and new // polynomials. - Difference> error_estimate{}; + Difference> current_error_estimate{}; for (std::int64_t j = 0; j <= N / 2; ++j) { - error_estimate += Abs(previous_aⱼ[j] - aⱼ[j]); + current_error_estimate += Abs(previous_aⱼ[j] - aⱼ[j]); } for (std::int64_t j = N / 2 + 1; j <= N; ++j) { - error_estimate += Abs(aⱼ[j]); + current_error_estimate += Abs(aⱼ[j]); } if constexpr (N <= max_degree) { - if (error_estimate > max_error) { + if (current_error_estimate > max_error) { // Note that this recursive call overflows the stack when // max_degree >= 256. We could allocate on the heap, but then we don't // care about very high degree polynomials. return ЧебышёвPolynomialInterpolantImplementation<2 * N, max_degree>( - f, a, b, max_error, fₖ, aⱼ); + f, a, b, max_error, fₖ, aⱼ, error_estimate); } } @@ -99,6 +100,9 @@ template // impose an unnecessary cost on the client (e.g., more costly evaluation). If // a client wants a more precise approximation, they just need to give a // smaller |max_error|. + if (error_estimate != nullptr) { + *error_estimate = current_error_estimate; + } std::vector> coefficients; std::copy(previous_aⱼ.begin(), previous_aⱼ.end(), std::back_inserter(coefficients)); @@ -112,7 +116,8 @@ template Function const& f, Argument const& lower_bound, Argument const& upper_bound, - Difference> const& max_error) { + Difference> const& max_error, + Difference>* const error_estimate) { auto const& a = lower_bound; auto const& b = upper_bound; auto const f_a = f(a); @@ -120,11 +125,12 @@ template FixedVector, 2> const fₖ({f_b, f_a}); FixedVector, 2> const aⱼ( {0.5 * (f_b + f_a), 0.5 * (f_b - f_a)}); + return ЧебышёвPolynomialInterpolantImplementation( - f, a, b, max_error, fₖ, aⱼ); + f, a, b, max_error, fₖ, aⱼ, error_estimate); } } // namespace internal diff --git a/numerics/approximation_test.cpp b/numerics/approximation_test.cpp index 6298512152..6b7ae18ee3 100644 --- a/numerics/approximation_test.cpp +++ b/numerics/approximation_test.cpp @@ -5,7 +5,9 @@ #include "quantities/elementary_functions.hpp" #include "quantities/named_quantities.hpp" #include "quantities/si.hpp" +#include "testing_utilities/is_near.hpp" #include "testing_utilities/numerics_matchers.hpp" +#include "testing_utilities/approximate_quantity.hpp" namespace principia { namespace numerics { @@ -14,6 +16,8 @@ namespace _approximation { using namespace principia::quantities::_elementary_functions; using namespace principia::quantities::_named_quantities; using namespace principia::quantities::_si; +using namespace principia::testing_utilities::_approximate_quantity; +using namespace principia::testing_utilities::_is_near; using namespace principia::testing_utilities::_numerics_matchers; using ::testing::AllOf; using ::testing::Ge; @@ -21,12 +25,16 @@ using ::testing::Lt; TEST(ApproximationTest, SinInverse) { auto const f = [](double const x) { return Sin(1 * Radian / x); }; + double error_estimate; auto const approximation = ЧебышёвPolynomialInterpolant<128>(f, /*lower_bound=*/0.1, /*upper_bound=*/10.0, - /*max_error=*/1e-6); + /*max_error=*/1e-6, + &error_estimate); EXPECT_EQ(128, approximation.degree()); + // Didn't reach the desired accuracy. + EXPECT_THAT(error_estimate, IsNear(4.9e-6_(1))); for (double x = 0.1; x < 10.0; x += 0.01) { EXPECT_THAT(approximation.Evaluate(x), AbsoluteErrorFrom(f(x), AllOf(Lt(3.0e-6), Ge(0)))); @@ -35,12 +43,15 @@ TEST(ApproximationTest, SinInverse) { TEST(ApproximationTest, Exp) { auto const f = [](double const x) { return std::exp(x); }; + double error_estimate; auto const approximation = ЧебышёвPolynomialInterpolant<128>(f, /*lower_bound=*/0.01, /*upper_bound=*/3.0, - /*max_error=*/1e-6); + /*max_error=*/1e-6, + &error_estimate); EXPECT_EQ(16, approximation.degree()); + EXPECT_THAT(error_estimate, IsNear(4.7e-14_(1))); for (double x = 0.01; x < 3; x += 0.01) { EXPECT_THAT(approximation.Evaluate(x), AbsoluteErrorFrom(f(x), AllOf(Lt(7.2e-15), Ge(0)))); From 8c9ea2e0232c9d75e9e18489bb4e4851f4b987d7 Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 20 Jan 2024 12:49:14 +0100 Subject: [PATCH 3/3] IWYUsing. --- numerics/approximation_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/numerics/approximation_test.cpp b/numerics/approximation_test.cpp index 6b7ae18ee3..b87b1639c2 100644 --- a/numerics/approximation_test.cpp +++ b/numerics/approximation_test.cpp @@ -5,9 +5,9 @@ #include "quantities/elementary_functions.hpp" #include "quantities/named_quantities.hpp" #include "quantities/si.hpp" +#include "testing_utilities/approximate_quantity.hpp" #include "testing_utilities/is_near.hpp" #include "testing_utilities/numerics_matchers.hpp" -#include "testing_utilities/approximate_quantity.hpp" namespace principia { namespace numerics {