Merge pull request #4267 from pleroy/MoreTests

Change the Sin/Cos test to use multithreading

Author: pleroy

functions/multiprecision.cpp

@@ -7,15 +7,19 @@ namespace internal {
 
 // The designers of the Boost multiprecision library, in their confusion, have
 // decided that nobody would use angles greater than about 1 / ε, so here I am,
-// using a 1000-bit rational approximation of 2π to do a first round of
+// using a 2000-bit rational approximation of 2π to do a first round of
 // reduction before calling their library.  Sigh.
 cpp_rational const two_π(
-    "74778118356356393099550700487949553821548169522945304260092193625666650672"
-    "12764554292012625975645707431952036396539900828099263122542425779444638301"
-    "713",
-    "11901307171524915725967286244795184945313912891969281885856459971598729125"
-    "90781909108254956448232159392109804270301891472482302074906413637232243084"
-    "548");
+    "85501617327051614281806889507174341018563881715359722726460240448778253362"
+    "28033574412989010393845826881894399169689030062477829223652983133581682809"
+    "07521562450510217951931573900693175198362946680952656259522439979797047648"
+    "97394953111420464708282657298132250506331749517560155388580727329016311453"
+    "71910",
+    "13608005039951911862130979398099086168163224029757063948525789956741291980"
+    "85326865353012588752613375780472488965315733291185793228273924312743154885"
+    "30992776567452143477216307689797948122382208381130766696418513975898402329"
+    "95420184160961245582170902313515516285224226723935852853290026779141303674"
+    "08461");
 
 cpp_bin_float_50 Sin(cpp_rational const& α) {
   auto const k = static_cast<cpp_int>(α / two_π);

functions/sin_cos_test.cpp

@@ -3,7 +3,9 @@
 #include <algorithm>
 #include <limits>
 #include <random>
+#include <vector>
 
+#include "base/bundle.hpp"
 #include "boost/multiprecision/cpp_int.hpp"
 #include "functions/multiprecision.hpp"
 #include "glog/logging.h"
@@ -21,6 +23,7 @@
 namespace principia {
 namespace functions_test {
 
+using namespace principia::base::_bundle;
 using namespace boost::multiprecision;
 using namespace principia::functions::_multiprecision;
 using namespace principia::numerics::_next;
@@ -30,32 +33,34 @@ namespace sin_cos = principia::numerics::_sin_cos;
 
 class SinCosTest : public ::testing::Test {
  protected:
+  struct FunctionStatistics {
+    cpp_bin_float_50 max_ulps_error = 0;
+    double worst_argument = 0;
+    cpp_bin_float_50 boost_fn_worst_argument = 0;
+    double principia_fn_worst_argument = 0;
+    std::int64_t incorrectly_rounded = 0;
+  };
+
+  struct SinCosStatistics {
+    FunctionStatistics sin;
+    FunctionStatistics cos;
+  };
+
   static void SetUpTestCase() {
     sin_cos::StaticInitialization();
   }
 
-  double a_ = 1.0;
-
-  static void RandomArgumentTest(double const lower_bound,
-                                 double const upper_bound) {
-    std::mt19937_64 random(42);
+  template<std::int64_t iterations_quantum>
+  static SinCosStatistics RandomArgumentTest(double const lower_bound,
+                                             double const upper_bound,
+                                             std::int64_t const seed) {
+    std::mt19937_64 random(seed);
     std::uniform_real_distribution<> uniformly_at(lower_bound, upper_bound);
     std::uniform_int_distribution<> uniform_sign(0, 1);
 
-    cpp_bin_float_50 max_sin_ulps_error = 0;
-    cpp_bin_float_50 max_cos_ulps_error = 0;
-    double worst_sin_argument = 0;
-    double worst_cos_argument = 0;
-    std::int64_t incorrectly_rounded_sin = 0;
-    std::int64_t incorrectly_rounded_cos = 0;
-
-#if _DEBUG
-    static constexpr std::int64_t iterations = 100;
-#else
-    static constexpr std::int64_t iterations = 250'000;
-#endif
+    SinCosStatistics s;
 
-    for (std::int64_t i = 0; i < iterations; ++i) {
+    for (std::int64_t i = 0; i < iterations_quantum; ++i) {
       double const principia_argument =
           uniformly_at(random) * ((uniform_sign(random) << 1) - 1);
       auto const boost_argument = cpp_rational(principia_argument);
@@ -66,12 +71,14 @@ class SinCosTest : public ::testing::Test {
             abs(boost_sin - static_cast<cpp_bin_float_50>(principia_sin));
         auto const ulp = NextUp(principia_sin) - principia_sin;
         auto const sin_ulps_error = sin_error / ulp;
-        if (sin_ulps_error > max_sin_ulps_error) {
-          max_sin_ulps_error = sin_ulps_error;
-          worst_sin_argument = principia_argument;
+        if (sin_ulps_error > s.sin.max_ulps_error) {
+          s.sin.max_ulps_error = sin_ulps_error;
+          s.sin.worst_argument = principia_argument;
+          s.sin.boost_fn_worst_argument = boost_sin;
+          s.sin.principia_fn_worst_argument = principia_sin;
         }
         if (sin_ulps_error > 0.5) {
-          ++incorrectly_rounded_sin;
+          ++s.sin.incorrectly_rounded;
           LOG(ERROR) << "Sin: " << sin_ulps_error << " ulps at "
                      << std::setprecision(25) << principia_argument;
         }
@@ -83,34 +90,91 @@ class SinCosTest : public ::testing::Test {
             abs(boost_cos - static_cast<cpp_bin_float_50>(principia_cos));
         auto const ulp = NextUp(principia_cos) - principia_cos;
         auto const cos_ulps_error = cos_error / ulp;
-        if (cos_ulps_error > max_cos_ulps_error) {
-          max_cos_ulps_error = cos_ulps_error;
-          worst_cos_argument = principia_argument;
+        if (cos_ulps_error > s.cos.max_ulps_error) {
+          s.cos.max_ulps_error = cos_ulps_error;
+          s.cos.worst_argument = principia_argument;
+          s.cos.boost_fn_worst_argument = boost_cos;
+          s.cos.principia_fn_worst_argument = principia_cos;
         }
         if (cos_ulps_error > 0.5) {
-          ++incorrectly_rounded_cos;
+          ++s.cos.incorrectly_rounded;
           LOG(ERROR) << "Cos: " << cos_ulps_error << " ulps at "
                      << std::setprecision(25) << principia_argument;
         }
       }
     }
 
-    // This implementation is not quite correctly rounded, but not far from it.
-    EXPECT_LE(max_sin_ulps_error, 0.5);
-    EXPECT_LE(max_cos_ulps_error, 0.5);
-    EXPECT_EQ(incorrectly_rounded_sin, 0);
-    EXPECT_EQ(incorrectly_rounded_cos, 0);
-
-    LOG(ERROR) << "Sin error: " << max_sin_ulps_error << std::setprecision(25)
-               << " ulps for argument: " << worst_sin_argument
-               << " value: " << Sin(worst_sin_argument)
-               << "; incorrectly rounded: " << std::setprecision(3)
-               << incorrectly_rounded_sin / static_cast<double>(iterations);
-    LOG(ERROR) << "Cos error: " << max_cos_ulps_error << std::setprecision(25)
-               << " ulps for argument: " << worst_cos_argument
-               << " value: " << Cos(worst_cos_argument)
-               << "; incorrectly rounded: " << std::setprecision(3)
-               << incorrectly_rounded_cos / static_cast<double>(iterations);
+    EXPECT_LE(s.sin.max_ulps_error, 0.5);
+    EXPECT_LE(s.cos.max_ulps_error, 0.5);
+    EXPECT_EQ(s.sin.incorrectly_rounded, 0);
+    EXPECT_EQ(s.cos.incorrectly_rounded, 0);
+
+    return s;
+  }
+
+  void ParallelRandomArgumentTest(double const lower_bound,
+                                  double const upper_bound) {
+#if _DEBUG
+    static constexpr std::int64_t iterations = 30'000;
+    static constexpr std::int64_t iterations_quantum = 100;
+#else
+    static constexpr std::int64_t iterations = 10'000'000;
+    static constexpr std::int64_t iterations_quantum = 10'000;
+#endif
+    static_assert(iterations % iterations_quantum == 0);
+
+    Bundle bundle;
+
+    std::vector<SinCosStatistics> statistics(iterations / iterations_quantum);
+    for (std::int64_t i = 0; i < statistics.size(); ++i) {
+      bundle.Add([i, lower_bound, upper_bound, &statistics]() {
+        statistics[i] = RandomArgumentTest<iterations_quantum>(
+            lower_bound, upper_bound, /*seed=*/i);
+        return absl::OkStatus();
+      });
+    }
+    absl::Status const status = bundle.Join();
+
+    SinCosStatistics final_statistics;
+    for (auto const& s : statistics) {
+      if (s.sin.max_ulps_error > final_statistics.sin.max_ulps_error) {
+        final_statistics.sin.max_ulps_error = s.sin.max_ulps_error;
+        final_statistics.sin.worst_argument = s.sin.worst_argument;
+        final_statistics.sin.boost_fn_worst_argument =
+            s.sin.boost_fn_worst_argument;
+        final_statistics.sin.principia_fn_worst_argument =
+            s.sin.principia_fn_worst_argument;
+      }
+      final_statistics.sin.incorrectly_rounded += s.sin.incorrectly_rounded;
+      if (s.cos.max_ulps_error > final_statistics.cos.max_ulps_error) {
+        final_statistics.cos.max_ulps_error = s.cos.max_ulps_error;
+        final_statistics.cos.worst_argument = s.cos.worst_argument;
+        final_statistics.cos.boost_fn_worst_argument =
+            s.cos.boost_fn_worst_argument;
+        final_statistics.cos.principia_fn_worst_argument =
+            s.cos.principia_fn_worst_argument;
+      }
+      final_statistics.cos.incorrectly_rounded += s.cos.incorrectly_rounded;
+    }
+
+    auto log_statistics = [](std::string_view const fn_name,
+                             FunctionStatistics const& s) {
+      LOG(ERROR) << fn_name << " error: " << s.max_ulps_error
+                 << std::setprecision(25)
+                 << " ulps for argument: " << s.worst_argument << " ("
+                 << std::hexfloat << s.worst_argument << std::defaultfloat
+                 << ") value: " << s.principia_fn_worst_argument << " ("
+                 << std::hexfloat << s.principia_fn_worst_argument
+                 << std::defaultfloat << ") vs. expected "
+                 << s.boost_fn_worst_argument << " (" << std::hexfloat
+                 << static_cast<double>(s.boost_fn_worst_argument)
+                 << std::defaultfloat << "); incorrectly rounded probability: "
+                 << std::setprecision(3)
+                 << s.incorrectly_rounded / static_cast<double>(iterations);
+    };
+
+    log_statistics("Sin", final_statistics.sin);
+    log_statistics("Cos", final_statistics.cos);
   }
 };
 
@@ -169,19 +233,20 @@ TEST_F(SinCosTest, ReduceIndex) {
 }
 
 TEST_F(SinCosTest, RandomSmall) {
-  RandomArgumentTest(0, π / 4);
+  ParallelRandomArgumentTest(0, π / 4);
 }
 
 TEST_F(SinCosTest, RandomTwoTerms) {
-  RandomArgumentTest(π / 4, 1 << 8);
+  ParallelRandomArgumentTest(π / 4, 1 << 8);
 }
 
 TEST_F(SinCosTest, RandomThreeTerms) {
-  RandomArgumentTest(1 << 8, 1 << 18);
+  ParallelRandomArgumentTest(1 << 8, 1 << 18);
 }
 
 TEST_F(SinCosTest, RandomLarge) {
-  RandomArgumentTest(1 << 18, std::numeric_limits<double>::max() / 1.0e30);
+  ParallelRandomArgumentTest(1 << 18,
+                             std::numeric_limits<double>::max() / 1.0e30);
 }
 
 // Values for which the base algorithm gives an error of 1 ULP.

mathematica/accurate_tables.wl

@@ -63,7 +63,7 @@
 (*};*)
 (**)
 (*constexpr std::array<SinCosAccurateValues, " <> ToString[max + 1] <> "> SinCosAccurateTable{{\n" <>*)
-(*StringRiffle[Join[Table["    /*"<>StringPadLeft[ToString[i],width]<>"*/{.x =     std::numeric_limits<double>::signaling_NaN(),\n"<>StringRepeat[" ",width+9]<>".sin_x = std::numeric_limits<double>::signaling_NaN(),\n"<>StringRepeat[" ",width+9]<>".cos_x = std::numeric_limits<double>::signaling_NaN()",{i,0,min-1}],Table["    /*"<>StringPadLeft[ToString[i],width]<>"*/{.x =     " <> hexFloatLiteral[sincos`accurateTables[i][[1]],False] <> ",\n"<>StringRepeat[" ",width+9]<>".sin_x = " <> hexFloatLiteral[sincos`accurateTables[i][[2]],False] <> ",\n"<>StringRepeat[" ",width+9]<>".cos_x = "<> hexFloatLiteral[sincos`accurateTables[i][[3]],False], {i,indices}]],"},\n"]<>"}}};*)
+(*StringRiffle[Join[Table["    /*"<>StringPadLeft[ToString[i],width]<>"*/{.x =     0,\n"<>StringRepeat[" ",width+9]<>".sin_x = 0,\n"<>StringRepeat[" ",width+9]<>".cos_x = 1",{i,0,min-1}],Table["    /*"<>StringPadLeft[ToString[i],width]<>"*/{.x =     " <> hexFloatLiteral[sincos`accurateTables[i][[1]],False] <> ",\n"<>StringRepeat[" ",width+9]<>".sin_x = " <> hexFloatLiteral[sincos`accurateTables[i][[2]],False] <> ",\n"<>StringRepeat[" ",width+9]<>".cos_x = "<> hexFloatLiteral[sincos`accurateTables[i][[3]],False], {i,indices}]],"},\n"]<>"}}};*)
 (**)
 (*}  // namespace internal*)
 (**)

mathematica/mathematica.vcxproj

@@ -28,6 +28,7 @@
     <ClInclude Include="mathematica_body.hpp" />
   </ItemGroup>
   <ItemGroup>
+    <None Include="accurate_tables.wl" />
     <None Include="associated_legendre_function.wl" />
     <None Include="elliptic_integrals.wl" />
     <None Include="ieee754_floating_point.wl" />

mathematica/mathematica.vcxproj.filters

@@ -102,5 +102,8 @@
     <None Include="ieee754_floating_point_test.wl">
       <Filter>Scripts</Filter>
     </None>
+    <None Include="accurate_tables.wl">
+      <Filter>Scripts</Filter>
+    </None>
   </ItemGroup>
 </Project>

numerics/accurate_tables.mathematica.h

@@ -16,9 +16,9 @@ struct SinCosAccurateValues {
 };
 
 constexpr std::array<SinCosAccurateValues, 403> SinCosAccurateTable{{
-    /*  0*/{.x =     std::numeric_limits<double>::signaling_NaN(),
-            .sin_x = std::numeric_limits<double>::signaling_NaN(),
-            .cos_x = std::numeric_limits<double>::signaling_NaN()},
+    /*  0*/{.x =     0,
+            .sin_x = 0,
+            .cos_x = 1},
     /*  1*/{.x =     0x0.007f'ffce'7f29'afe0'0000'0000p0,
             .sin_x = 0x0.007f'ffc9'29da'9bb4'0000'db97p0,
             .cos_x = 0x0.ffff'e000'196b'1000'0148'882dp0},

numerics/double_precision_body.hpp

@@ -102,7 +102,11 @@ template<typename T, typename U>
   requires convertible_to_quantity<T> && convertible_to_quantity<U>
 struct ComponentwiseComparator<T, U> {
   static bool GreaterThanOrEqualOrZero(T const& left, U const& right) {
-    return Abs(left) >= Abs(right) || left == T{} || !IsFinite(left);
+    // In the elementary functions, we use NaN to fall back to the CORE-MATH
+    // implementation.  We don't want to die because of the weird comparisons of
+    // NaNs.
+    return Abs(left) >= Abs(right) || left == T{} ||
+           !IsFinite(left) || !IsFinite(right);
   }
 };
 

numerics/numerics.vcxproj

@@ -26,6 +26,7 @@
     </ClCompile>
   </ItemDefinitionGroup>
   <ItemGroup>
+    <ClInclude Include="accurate_tables.mathematica.h" />
     <ClInclude Include="angle_reduction.hpp" />
     <ClInclude Include="angle_reduction_body.hpp" />
     <ClInclude Include="approximation.hpp" />

numerics/numerics.vcxproj.filters

@@ -272,6 +272,12 @@
     <ClInclude Include="sin_cos.hpp">
       <Filter>Header Files</Filter>
     </ClInclude>
+    <ClInclude Include="osaca.hpp">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="accurate_tables.mathematica.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="fixed_arrays_test.cpp">

numerics/sin_cos.cpp

@@ -159,15 +159,19 @@ double FusedNegatedMultiplyAdd(double const a, double const b, double const c) {
 // negative.
 template<FMAPolicy fma_policy, double e>
 double DetectDangerousRounding(double const x, double const Δx) {
+  // We don't check that `Δx` is not NaN because that's how we trigger fallback
+  // to the slow path.
+  DCHECK(x == x);
   DoublePrecision<double> const sum = QuickTwoSum(x, Δx);
   double const& value = sum.value;
   double const& error = sum.error;
   OSACA_IF(value == FusedMultiplyAdd<fma_policy>(error, e, value)) {
     return value;
   } else {
 #if _DEBUG
-    LOG(ERROR) << std::setprecision(25) << x << " " << std::hexfloat << value
-               << " " << error << " " << e;
+    LOG_IF(ERROR, value == value && error == error)
+        << std::setprecision(25) << x << " " << std::hexfloat << value << " "
+        << error << " " << e;
 #endif
     return std::numeric_limits<double>::quiet_NaN();
   }