Stehlé-Zimmermann multisearch

Principia.sln

@@ -682,8 +682,8 @@ Global
 		{7CCA653C-2E8F-4FFD-9E9F-BEE590F3EFAB}.Release KSP 1.7.3|x86.Build.0 = Release|Win32
 		{7CCA653C-2E8F-4FFD-9E9F-BEE590F3EFAB}.Release_LLVM|Any CPU.ActiveCfg = Release|x64
 		{7CCA653C-2E8F-4FFD-9E9F-BEE590F3EFAB}.Release_LLVM|Any CPU.Build.0 = Release|x64
-		{7CCA653C-2E8F-4FFD-9E9F-BEE590F3EFAB}.Release_LLVM|x64.ActiveCfg = Release|x64
-		{7CCA653C-2E8F-4FFD-9E9F-BEE590F3EFAB}.Release_LLVM|x64.Build.0 = Release|x64
+		{7CCA653C-2E8F-4FFD-9E9F-BEE590F3EFAB}.Release_LLVM|x64.ActiveCfg = Release_LLVM|x64
+		{7CCA653C-2E8F-4FFD-9E9F-BEE590F3EFAB}.Release_LLVM|x64.Build.0 = Release_LLVM|x64
 		{7CCA653C-2E8F-4FFD-9E9F-BEE590F3EFAB}.Release_LLVM|x86.ActiveCfg = Release|Win32
 		{7CCA653C-2E8F-4FFD-9E9F-BEE590F3EFAB}.Release_LLVM|x86.Build.0 = Release|Win32
 		{7CCA653C-2E8F-4FFD-9E9F-BEE590F3EFAB}.Release|Any CPU.ActiveCfg = Release|x64

functions/accurate_table_generator.hpp

@@ -31,32 +31,41 @@ std::vector<cpp_rational> GalExhaustiveMultisearch(
     std::vector<AccurateFunction> const& functions,
     std::vector<cpp_rational> const& starting_arguments);
 
-// Searches in an interval of radius |T / N| centered on |near_argument|.  The
-// |polynomials| must be the degree-2 Taylor approximations of the |functions|.
-// The argument and function values must be within [1/2, 1[.
+// 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[.
 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,
-    cpp_rational const& near_argument,
+    cpp_rational const& starting_argument,
     std::int64_t N,
     std::int64_t T);
 
-// Performs a search around |near_argument| to find a solution, automatically
-// adjusting the interval over which the search happens.  The argument and
-// function values must be nonzero.
+// Performs a search around |starting_argument| to find a solution,
+// automatically adjusting the interval over which the search happens.  The
+// argument and function values must be nonzero.
 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,
-    cpp_rational const& near_argument);
+    cpp_rational const& starting_argument);
+
+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<cpp_rational> const& starting_arguments);
 
 }  // namespace internal
 
 using internal::AccuratePolynomial;
 using internal::GalExhaustiveMultisearch;
 using internal::GalExhaustiveSearch;
 using internal::StehléZimmermannSimultaneousFullSearch;
+using internal::StehléZimmermannSimultaneousMultisearch;
 using internal::StehléZimmermannSimultaneousSearch;
 
 }  // namespace _accurate_table_generator

functions/accurate_table_generator_body.hpp

@@ -138,7 +138,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,
-    cpp_rational const& near_argument,
+    cpp_rational const& starting_argument,
     std::int64_t const N,
     std::int64_t const T) {
   // This implementation follows [SZ05], section 3.1.
@@ -150,13 +150,13 @@ absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousSearch(
   std::array<AccurateFunction, 2> F;
   std::array<std::optional<AccuratePolynomial<cpp_rational, 2>>, 2> P;
   for (std::int64_t i = 0; i < functions.size(); ++i) {
-    F[i] = [&functions, i, N, &near_argument](cpp_rational const& t) {
+    F[i] = [&functions, i, N, &starting_argument](cpp_rational const& t) {
       // Here |t| <= |T|.
-      return N * functions[i](near_argument + t / N);
+      return N * functions[i](starting_argument + t / N);
     };
   }
   AccuratePolynomial<cpp_rational, 1> const shift_and_rescale(
-      {near_argument, cpp_rational(1, N)});
+      {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);
   }
@@ -303,14 +303,14 @@ absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousSearch(
     }
   }
 
-  return near_argument + t₀ / N;
+  return starting_argument + t₀ / N;
 }
 
 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,
-    cpp_rational const& near_argument) {
+    cpp_rational const& starting_argument) {
   constexpr std::int64_t M = 1LL << zeroes;
   constexpr std::int64_t N = 1LL << std::numeric_limits<double>::digits;
 
@@ -322,17 +322,17 @@ absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousFullSearch(
   // Scale the argument, functions, and polynomials to lie within [1/2, 1[.
   // TODO(phl): Handle an argument that is exactly a power of 2.
   std::int64_t argument_exponent;
-  auto const argument_mantissa =
-      frexp(static_cast<cpp_bin_float_50>(near_argument), &argument_exponent);
+  auto const argument_mantissa = frexp(
+      static_cast<cpp_bin_float_50>(starting_argument), &argument_exponent);
   CHECK_NE(0, argument_mantissa);
   auto const argument_scale = exp2(-argument_exponent);
-  cpp_rational const scaled_argument = near_argument * argument_scale;
+  cpp_rational const scaled_argument = starting_argument * argument_scale;
 
   std::array<AccurateFunction, 2> scaled_functions;
   for (std::int64_t i = 0; i < scaled_functions.size(); ++i) {
     std::int64_t function_exponent;
     auto const function_mantissa =
-        frexp(functions[i](near_argument), &function_exponent);
+        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](
@@ -342,11 +342,11 @@ absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousFullSearch(
   }
 
   auto build_scaled_polynomial =
-      [argument_scale,
-       &near_argument](AccuratePolynomial<cpp_rational, 2> const& polynomial) {
+      [argument_scale, &starting_argument](
+          AccuratePolynomial<cpp_rational, 2> const& polynomial) {
         std::int64_t polynomial_exponent;
         auto const polynomial_mantissa =
-            frexp(static_cast<cpp_bin_float_50>(polynomial(near_argument)),
+            frexp(static_cast<cpp_bin_float_50>(polynomial(starting_argument)),
                   &polynomial_exponent);
         CHECK_NE(0, polynomial_mantissa);
         return exp2(-polynomial_exponent) *
@@ -430,6 +430,32 @@ absl::StatusOr<cpp_rational> StehléZimmermannSimultaneousFullSearch(
   }
 }
 
+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<cpp_rational> const& starting_arguments) {
+  ThreadPool<absl::StatusOr<cpp_rational>> search_pool(
+      std::thread::hardware_concurrency());
+
+  std::vector<std::future<absl::StatusOr<cpp_rational>>> 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]);
+        }));
+  }
+
+  std::vector<absl::StatusOr<cpp_rational>> results;
+  for (auto& future : futures) {
+    results.push_back(future.get());
+  }
+  return results;
+}
+
 }  // namespace internal
 }  // namespace _accurate_table_generator
 }  // namespace functions

functions/accurate_table_generator_test.cpp

@@ -20,6 +20,7 @@ namespace principia {
 namespace functions {
 namespace _accurate_table_generator {
 
+using ::testing::AnyOf;
 using ::testing::Eq;
 using ::testing::Lt;
 using ::testing::SizeIs;
@@ -78,7 +79,7 @@ TEST_F(AccurateTableGeneratorTest, GalSinCos5) {
   }
 }
 
-TEST_F(AccurateTableGeneratorTest, GalSinCos5Multisearch) {
+TEST_F(AccurateTableGeneratorTest, GalMultisearchSinCos5) {
   static constexpr std::int64_t index_begin = 17;
   static constexpr std::int64_t index_end = 100;
   std::vector<cpp_rational> starting_arguments;
@@ -243,6 +244,59 @@ TEST_F(AccurateTableGeneratorTest, StehléZimmermannFullSinCos15WithScaling) {
   }
 }
 
+TEST_F(AccurateTableGeneratorTest, StehléZimmermannMultisearchSinCos15) {
+  FLAGS_v = 0;
+  google::LogToStderr();
+  static constexpr std::int64_t index_begin = 17;
+  static constexpr std::int64_t index_end = 100;
+  std::vector<cpp_rational> starting_arguments;
+  std::vector<std::array<AccuratePolynomial<cpp_rational, 2>, 2>> polynomials;
+  for (std::int64_t i = index_begin; i < index_end; ++i) {
+    auto const x₀ = i / 128.0;
+    AccuratePolynomial<cpp_rational, 2> const sin_taylor2(
+        {cpp_rational(Sin(x₀)),
+         cpp_rational(Cos(x₀)),
+         -cpp_rational(Sin(x₀)) / 2},
+        x₀);
+    AccuratePolynomial<cpp_rational, 2> const cos_taylor2(
+        {cpp_rational(Cos(x₀)),
+         -cpp_rational(Sin(x₀)),
+         -cpp_rational(Cos(x₀) / 2)},
+        x₀);
+    starting_arguments.push_back(x₀);
+    polynomials.push_back({sin_taylor2, cos_taylor2});
+  }
+  auto const xs = StehléZimmermannSimultaneousMultisearch<15>(
+      {Sin, Cos}, polynomials, starting_arguments);
+  EXPECT_THAT(xs, SizeIs(index_end - index_begin));
+  for (std::int64_t i = 0; i < xs.size(); ++i) {
+    CHECK_OK(xs[i].status());
+    auto const& x = *xs[i];
+    EXPECT_THAT(static_cast<double>(x),
+                RelativeErrorFrom((i + index_begin) / 128.0, Lt(1.3e-7)));
+    {
+      std::string const mathematica = ToMathematica(Sin(x),
+                                                    /*express_in=*/std::nullopt,
+                                                    /*base=*/2);
+      std::string_view mantissa = mathematica;
+      CHECK(absl::ConsumePrefix(&mantissa, "Times[2^^"));
+      EXPECT_THAT(mantissa.substr(53, 15),
+                  AnyOf(Eq("00000""00000""00000"),
+                        Eq("11111""11111""11111")));
+    }
+    {
+      std::string const mathematica = ToMathematica(Cos(x),
+                                                    /*express_in=*/std::nullopt,
+                                                    /*base=*/2);
+      std::string_view mantissa = mathematica;
+      CHECK(absl::ConsumePrefix(&mantissa, "Times[2^^"));
+      EXPECT_THAT(mantissa.substr(53, 15),
+                  AnyOf(Eq("00000""00000""00000"),
+                        Eq("11111""11111""11111")));
+    }
+  }
+}
+
 #endif
 
 }  // namespace _accurate_table_generator