Merge pull request #4109 from pleroy/SinCos

A simplistic implementation of Sin and Cos over [0, π/4]

functions/functions.vcxproj

@@ -17,6 +17,7 @@
     <ClCompile Include="cos.cc" />
     <ClCompile Include="multiprecision.cpp" />
     <ClCompile Include="sin.cc" />
+    <ClCompile Include="sin_cos_test.cpp" />
     <ClCompile Include="std_accuracy_test.cpp" />
   </ItemGroup>
   <Import Project="$(SolutionDir)principia.props" />

functions/functions.vcxproj.filters

@@ -36,6 +36,9 @@
     <ClCompile Include="sin.cc">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="sin_cos_test.cpp">
+      <Filter>Test Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="multiprecision.hpp">

functions/sin_cos_test.cpp

@@ -0,0 +1,85 @@
+#include "numerics/sin_cos.hpp"
+
+#include <algorithm>
+#include <limits>
+#include <random>
+
+#include "boost/multiprecision/cpp_int.hpp"
+#include "functions/multiprecision.hpp"
+#include "glog/logging.h"
+#include "gtest/gtest.h"
+#include "numerics/next.hpp"
+#include "quantities/numbers.hpp"
+
+// This test lives in `functions` to avoid pulling `boost` into `numerics`.
+namespace principia {
+namespace numerics {
+namespace _sin_cos {
+
+using namespace boost::multiprecision;
+using namespace principia::numerics::_next;
+
+class SinCosTest : public ::testing::Test {};
+
+TEST_F(SinCosTest, Random) {
+  std::mt19937_64 random(42);
+  // TODO(phl): Negative angles.
+  std::uniform_real_distribution<> uniformly_at(0, π / 4);
+
+  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;
+
+#if _DEBUG
+  static constexpr std::int64_t iterations = 100;
+#else
+  static constexpr std::int64_t iterations = 400'000;
+#endif
+
+  for (std::int64_t i = 0; i < iterations; ++i) {
+    double const principia_argument = uniformly_at(random);
+    auto const boost_argument = cpp_rational(principia_argument);
+    {
+      auto const boost_sin =
+          functions::_multiprecision::Sin(boost_argument);
+      double const principia_sin = Sin(principia_argument);
+      auto const sin_error =
+          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;
+      }
+    }
+    {
+      auto const boost_cos =
+          functions::_multiprecision::Cos(boost_argument);
+      double const principia_cos = Cos(principia_argument);
+      auto const cos_error =
+          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;
+      }
+    }
+  }
+
+  // This implementation is not quite correctly rounded, but not far from it.
+  EXPECT_LE(max_sin_ulps_error, 0.500002);
+  EXPECT_LE(max_cos_ulps_error, 0.500002);
+
+  LOG(ERROR) << "Sin error: " << max_sin_ulps_error << std::setprecision(25)
+             << " ulps for argument: " << worst_sin_argument
+             << " value: " << Sin(worst_sin_argument);
+  LOG(ERROR) << "Cos error: " << max_cos_ulps_error << std::setprecision(25)
+             << " ulps for argument: " << worst_cos_argument
+             << " value: " << Cos(worst_cos_argument);
+}
+
+}  // namespace _sin_cos
+}  // namespace numerics
+}  // namespace principia

mathematica/accurate_tables.wl

@@ -31,7 +31,7 @@
 (*group=5,*)
 (*i=IntegerPart[N[x,100]],*)
 (*f=FractionalPart[N[x,100]]},*)
-(*  StringJoin[If[signed,If[x<0,"-","+"],""],"0x",IntegerString[i,16],".",StringRiffle[StringPartition[IntegerString[Round[f*16^(group*groups)],16,group*groups],UpTo[groups]],"'"]]]*)
+(*  StringJoin[If[signed,If[x<0,"-","+"],""],"0x",IntegerString[i,16],".",StringRiffle[StringPartition[IntegerString[Round[f*16^(group*groups)],16,group*groups],UpTo[groups]],"'"],"p0"]]*)
 
 
 (* ::Text:: *)
@@ -43,6 +43,7 @@
 (*Module[{indices=Flatten[Position[ListQ/@Table[sincos`accurateTables[i],{i,0,500}],True]]-1,min,max,width},min=Min[indices];max=Max[indices];width=Ceiling[Log10[max]];*)
 (*"#pragma once*)
 (**)
+(*#include <array>*)
 (*#include <limits>*)
 (**)
 (*namespace principia {*)
@@ -56,8 +57,8 @@
 (*  double cos_x;*)
 (*};*)
 (**)
-(*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"]<>"}};*)
+(*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"]<>"}}};*)
 (**)
 (*}  // namespace internal*)
 (**)