Skip to content

Commit

Permalink
Merge pull request #4041 from pleroy/GenerateTheTables
Browse files Browse the repository at this point in the history 
A unit test for generating accurate tables
  • Loading branch information
@pleroy
pleroy authored Jul 20, 2024
2 parents 4312fcf + 679ec9a commit fefff88
Show file tree
Hide file tree
Showing 7 changed files with 352 additions and 49 deletions.
1 change: 1 addition & 0 deletions Principia.sln
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,7 @@ EndProject
Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Solution Items", "Solution Items", "{8EF5AB3D-F63B-48B9-80A0-B789A830413C}"
ProjectSection(SolutionItems) = preProject
.editorconfig = .editorconfig
principia.reg = principia.reg
EndProjectSection
EndProject
Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "astronomy", "astronomy\astronomy.vcxproj", "{5EEA5210-0FA9-4B44-8466-C40B22D24E8E}"
Expand Down
24 changes: 23 additions & 1 deletion functions/accurate_table_generator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,11 +33,13 @@ std::vector<cpp_rational> GalExhaustiveMultisearch(

// Searches in an interval of radius |T / N| centered on |starting_argument|.
// The |polynomials| must be the degree-2 Taylor approximations of the
// |functions|. The argument and function values must be within [1/2, 1[.
// |functions| and the |remainders| must be upper bounds on the remainder of the
// Taylor series. The argument and function values must be within [1/2, 1[.
template<std::int64_t zeroes>
absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousSearch(
std::array<AccurateFunction, 2> const& functions,
std::array<AccuratePolynomial<cpp_rational, 2>, 2> const& polynomials,
std::array<AccurateFunction, 2> const& remainders,
cpp_rational const& starting_argument,
std::int64_t N,
std::int64_t T);
Expand All @@ -49,24 +51,44 @@ template<std::int64_t zeroes>
absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousFullSearch(
std::array<AccurateFunction, 2> const& functions,
std::array<AccuratePolynomial<cpp_rational, 2>, 2> const& polynomials,
std::array<AccurateFunction, 2> const& remainders,
cpp_rational const& starting_argument);

// Same as above, but performs searches in parallel using the corresponding
// |polynomials|, |remainders|, and |starting_arguments|. Returns the results
// in the same order as the parameters.
template<std::int64_t zeroes>
std::vector<absl::StatusOr<cpp_rational>>
StehléZimmermannSimultaneousMultisearch(
std::array<AccurateFunction, 2> const& functions,
std::vector<std::array<AccuratePolynomial<cpp_rational, 2>, 2>> const&
polynomials,
std::vector<std::array<AccurateFunction, 2>> const& remainders,
std::vector<cpp_rational> const& starting_arguments);

// Same as above, but instead of accumulating all the results and returning them
// in a vector, it runs |callback| each time a computation is complete. The
// |index| indicates to which parameters the result corresponds.
template<std::int64_t zeroes>
void StehléZimmermannSimultaneousStreamingMultisearch(
std::array<AccurateFunction, 2> const& functions,
std::vector<std::array<AccuratePolynomial<cpp_rational, 2>, 2>> const&
polynomials,
std::vector<std::array<AccurateFunction, 2>> const& remainders,
std::vector<cpp_rational> const& starting_arguments,
std::function<void(/*index=*/std::int64_t,
absl::StatusOr<cpp_rational>)> const& callback);

} // namespace internal

using internal::AccurateFunction;
using internal::AccuratePolynomial;
using internal::GalExhaustiveMultisearch;
using internal::GalExhaustiveSearch;
using internal::StehléZimmermannSimultaneousFullSearch;
using internal::StehléZimmermannSimultaneousMultisearch;
using internal::StehléZimmermannSimultaneousSearch;
using internal::StehléZimmermannSimultaneousStreamingMultisearch;

} // namespace _accurate_table_generator
} // namespace functions
Expand Down
170 changes: 135 additions & 35 deletions functions/accurate_table_generator_body.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,17 +8,18 @@
#include <limits>
#include <memory>
#include <thread>
#include <utility>
#include <vector>

#include "base/bits.hpp"
#include "base/for_all_of.hpp"
#include "base/status_utilities.hpp" // 🧙 For RETURN_IF_ERROR.
#include "base/tags.hpp"
#include "base/thread_pool.hpp"
#include "geometry/interval.hpp"
#include "glog/logging.h"
#include "numerics/fixed_arrays.hpp"
#include "numerics/lattices.hpp"
#include "numerics/matrix_computations.hpp"
#include "numerics/matrix_views.hpp"
#include "quantities/elementary_functions.hpp"
#include "quantities/quantities.hpp"
Expand All @@ -35,12 +36,11 @@ using namespace principia::base::_thread_pool;
using namespace principia::geometry::_interval;
using namespace principia::numerics::_fixed_arrays;
using namespace principia::numerics::_lattices;
using namespace principia::numerics::_matrix_computations;
using namespace principia::numerics::_matrix_views;
using namespace principia::quantities::_elementary_functions;
using namespace principia::quantities::_quantities;

constexpr std::int64_t ε_computation_points = 16;
constexpr std::int64_t T_max = 16;

template<int rows, int columns>
FixedMatrix<cpp_int, rows, columns> ToInt(
Expand Down Expand Up @@ -91,6 +91,50 @@ bool AllFunctionValuesHaveDesiredZeroes(
});
}

// This is essentially the same as Gal's exhaustive search, but with the
// normalization done for [SZ05].
absl::StatusOr<std::int64_t> StehléZimmermannExhaustiveSearch(
std::array<AccurateFunction, 2> const& F,
std::int64_t const M,
std::int64_t const T) {
VLOG(2) << "Exhaustive search with T = " << T;
for (std::int64_t t = 0; t <= T; ++t) {
{
bool found = true;
for (auto const& Fᵢ : F) {
auto const Fᵢ_t = Fᵢ(t);
auto const Fᵢ_t_cmod_1 = Fᵢ_t - round(Fᵢ_t);
VLOG(2) << "Fi(t) cmod 1 = " << Fᵢ_t_cmod_1;
if (M * abs(Fᵢ_t_cmod_1) >= 1) {
found = false;
break;
}
}
if (found) {
VLOG(2) << "t = " << t;
return t;
}
}
if (t > 0) {
bool found = true;
for (auto const& Fᵢ : F) {
auto const Fᵢ_minus_t = Fᵢ(-t);
auto const Fᵢ_minus_t_cmod_1 = Fᵢ_minus_t - round(Fᵢ_minus_t);
VLOG(2) << "Fi(-t) cmod 1 = " << Fᵢ_minus_t_cmod_1;
if (M * abs(Fᵢ_minus_t_cmod_1) >= 1) {
found = false;
break;
}
}
if (found) {
VLOG(2) << "t = " << -t;
return -t;
}
}
}
return absl::NotFoundError("Not enough zeroes");
}

template<std::int64_t zeroes>
cpp_rational GalExhaustiveSearch(std::vector<AccurateFunction> const& functions,
cpp_rational const& starting_argument) {
Expand Down Expand Up @@ -145,6 +189,7 @@ template<std::int64_t zeroes>
absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousSearch(
std::array<AccurateFunction, 2> const& functions,
std::array<AccuratePolynomial<cpp_rational, 2>, 2> const& polynomials,
std::array<AccurateFunction, 2> const& remainders,
cpp_rational const& starting_argument,
std::int64_t const N,
std::int64_t const T) {
Expand All @@ -162,21 +207,32 @@ absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousSearch(
return N * functions[i](starting_argument + t / N);
};
}

// If the interval is small enough, we don't use the Stehlé-Zimmermann
// algorithm. Instead, we use an exhaustive search. Note that this may yield
// a better solution, because if there is one in the interval, it is sure to
// find it, whereas Stehlé-Zimmermann may miss it.
if (T <= T_max) {
auto const status_or_t = StehléZimmermannExhaustiveSearch(F, M, T);
RETURN_IF_ERROR(status_or_t);

return starting_argument + cpp_rational(status_or_t.value(), N);
}

AccuratePolynomial<cpp_rational, 1> const shift_and_rescale(
{starting_argument, cpp_rational(1, N)});
for (std::int64_t i = 0; i < polynomials.size(); ++i) {
P[i] = N * Compose(polynomials[i], shift_and_rescale);
}

// Step 2: compute ε. We don't care too much about its accuracy because in
// general the largest error is on the boundary of the domain, and anyway ε
// has virtually no incidence on the value of Mʹ.
auto const T_increment = cpp_rational(T, ε_computation_points);
// 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;
for (std::int64_t i = 0; i < 2; ++i) {
for (cpp_rational t = -T; t <= T; t += T_increment) {
ε = std::max(ε, abs(F[i](t) - static_cast<cpp_bin_float_50>((*P[i])(t))));
}
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)));
}
VLOG(2) << "ε = " << ε;

Expand Down Expand Up @@ -254,7 +310,6 @@ absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousSearch(
};

static constexpr std::int64_t dimension = 3;
FixedMatrix<cpp_rational, 5, dimension> w;
for (std::int64_t i = 0; i < dimension; ++i) {
auto const& vᵢ = *v[i];
VLOG(2) << "v[" << i << "] = " << vᵢ;
Expand Down Expand Up @@ -316,6 +371,7 @@ template<std::int64_t zeroes>
absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousFullSearch(
std::array<AccurateFunction, 2> const& functions,
std::array<AccuratePolynomial<cpp_rational, 2>, 2> const& polynomials,
std::array<AccurateFunction, 2> const& remainders,
cpp_rational const& starting_argument) {
constexpr std::int64_t M = 1LL << zeroes;
constexpr std::int64_t N = 1LL << std::numeric_limits<double>::digits;
Expand All @@ -334,35 +390,41 @@ absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousFullSearch(
auto const argument_scale = exp2(-argument_exponent);
cpp_rational const scaled_argument = starting_argument * argument_scale;

std::array<double, 2> function_scales;
std::array<AccurateFunction, 2> scaled_functions;
std::array<AccurateFunction, 2> scaled_remainders;
for (std::int64_t i = 0; i < scaled_functions.size(); ++i) {
std::int64_t function_exponent;
auto const function_mantissa =
frexp(functions[i](starting_argument), &function_exponent);
CHECK_NE(0, function_mantissa);
auto const function_scale = exp2(-function_exponent);
scaled_functions[i] = [argument_scale, function_scale, &functions, i](
cpp_rational const& argument) {
return functions[i](argument / argument_scale) * function_scale;
function_scales[i] = exp2(-function_exponent);
scaled_functions[i] = [argument_scale,
function_scale = function_scales[i],
&functions,
i](cpp_rational const& argument) {
return function_scale * functions[i](argument / argument_scale);
};
scaled_remainders[i] = [argument_scale,
function_scale = function_scales[i],
&remainders,
i](cpp_rational const& argument) {
return function_scale * remainders[i](argument / argument_scale);
};
}

auto build_scaled_polynomial =
[argument_scale, &starting_argument](
double const function_scale,
AccuratePolynomial<cpp_rational, 2> const& polynomial) {
std::int64_t polynomial_exponent;
auto const polynomial_mantissa =
frexp(static_cast<cpp_bin_float_50>(polynomial(starting_argument)),
&polynomial_exponent);
CHECK_NE(0, polynomial_mantissa);
return exp2(-polynomial_exponent) *
return function_scale *
Compose(polynomial,
AccuratePolynomial<cpp_rational, 1>(
{0, 1 / argument_scale}));
};
std::array<AccuratePolynomial<cpp_rational, 2>, 2> const scaled_polynomials{
build_scaled_polynomial(polynomials[0]),
build_scaled_polynomial(polynomials[1])};
build_scaled_polynomial(function_scales[0], polynomials[0]),
build_scaled_polynomial(function_scales[1], polynomials[1])};

// We construct intervals above and below |scaled_argument| and search for
// solutions on each side alternatively.
Expand All @@ -380,6 +442,7 @@ absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousFullSearch(
auto const status_or_solution =
StehléZimmermannSimultaneousSearch<zeroes>(scaled_functions,
scaled_polynomials,
scaled_remainders,
high_interval.midpoint(),
N,
T);
Expand Down Expand Up @@ -416,6 +479,7 @@ absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousFullSearch(
auto const status_or_solution =
StehléZimmermannSimultaneousSearch<zeroes>(scaled_functions,
scaled_polynomials,
scaled_remainders,
low_interval.midpoint(),
N,
T);
Expand Down Expand Up @@ -462,24 +526,60 @@ StehléZimmermannSimultaneousMultisearch(
std::array<AccurateFunction, 2> const& functions,
std::vector<std::array<AccuratePolynomial<cpp_rational, 2>, 2>> const&
polynomials,
std::vector<std::array<AccurateFunction, 2>> const& remainders,
std::vector<cpp_rational> const& starting_arguments) {
ThreadPool<absl::StatusOr<cpp_rational>> search_pool(
std::thread::hardware_concurrency());
std::vector<absl::StatusOr<cpp_rational>> result;
result.resize(starting_arguments.size());
StehléZimmermannSimultaneousStreamingMultisearch<zeroes>(
functions,
polynomials,
remainders,
starting_arguments,
[&result](std::int64_t const index,
absl::StatusOr<cpp_rational> status_or_final_argument) {
result[index] = std::move(status_or_final_argument);
});
return result;
}

std::vector<std::future<absl::StatusOr<cpp_rational>>> futures;
template<std::int64_t zeroes>
void StehléZimmermannSimultaneousStreamingMultisearch(
std::array<AccurateFunction, 2> const& functions,
std::vector<std::array<AccuratePolynomial<cpp_rational, 2>, 2>> const&
polynomials,
std::vector<std::array<AccurateFunction, 2>> const& remainders,
std::vector<cpp_rational> const& starting_arguments,
std::function<void(/*index=*/std::int64_t,
absl::StatusOr<cpp_rational>)> const& callback) {
ThreadPool<void> search_pool(std::thread::hardware_concurrency());

std::vector<std::future<void>> futures;
for (std::int64_t i = 0; i < starting_arguments.size(); ++i) {
futures.push_back(
search_pool.Add([i, &functions, &polynomials, &starting_arguments]() {
return StehléZimmermannSimultaneousFullSearch<zeroes>(
functions, polynomials[i], starting_arguments[i]);
}));
futures.push_back(search_pool.Add([i,
&callback,
&functions,
&polynomials,
&remainders,
&starting_arguments]() {
auto const& starting_argument = starting_arguments[i];
LOG(INFO) << "Starting search around " << starting_argument;
auto status_or_final_argument =
StehléZimmermannSimultaneousFullSearch<zeroes>(
functions, polynomials[i], remainders[i], starting_argument);
if (status_or_final_argument.ok()) {
LOG(INFO) << "Finished search around " << starting_argument
<< ", found " << status_or_final_argument.value();
} else {
LOG(WARNING) << "Search around " << starting_argument << " failed with"
<< status_or_final_argument.status();
}
callback(i, std::move(status_or_final_argument));
}));
}

std::vector<absl::StatusOr<cpp_rational>> results;
for (auto& future : futures) {
results.push_back(future.get());
future.wait();
}
return results;
}

} // namespace internal
Expand Down
Loading

0 comments on commit fefff88

Please sign in to comment.