Merge pull request #4108 from pleroy/C++Tables

Mathematica code for generating C++ tables

functions/accurate_table_generator_test.cpp

@@ -451,7 +451,8 @@ TEST_F(AccurateTableGeneratorTest, StehléZimmermannMultisearchSinCos15) {
 }
 
 TEST_F(AccurateTableGeneratorTest, DISABLED_SECULAR_SinCos18) {
-  Logger logger(TEMP_DIR / absl::StrCat("sin_cos_18.wl"),
+  static constexpr std::int64_t bits = 18;
+  Logger logger(TEMP_DIR / absl::StrCat("sin_cos_", bits, ".wl"),
                 /*make_unique=*/false);
 
   // The radius of each interval.
@@ -518,17 +519,40 @@ TEST_F(AccurateTableGeneratorTest, DISABLED_SECULAR_SinCos18) {
     remainders.push_back({remainder_sin_taylor2, remainder_cos_taylor2});
   }
 
-  StehléZimmermannSimultaneousStreamingMultisearch<18>(
+  StehléZimmermannSimultaneousStreamingMultisearch<bits>(
       {accurate_sin, accurate_cos},
       polynomials,
       remainders,
       starting_arguments,
       [i_min, &logger](std::int64_t const index,
                        absl::StatusOr<cpp_rational> status_or_x) {
         auto const& x = status_or_x.value();
-        logger.Set(
-            absl::StrCat("accurateTables[", index + i_min, "]"),
-            std::tuple{static_cast<cpp_bin_float_50>(x), Sin(x), Cos(x)});
+        auto const sin_x = Sin(x);
+        auto const cos_x = Cos(x);
+        {
+          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, bits),
+                      AnyOf(Eq("00000""00000""00000""000"),
+                            Eq("11111""11111""11111""111")));
+        }
+        {
+          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, bits),
+                      AnyOf(Eq("00000""00000""00000""000"),
+                            Eq("11111""11111""11111""111")));
+        }
+        logger.Set(absl::StrCat("accurateTables[", index + i_min, "]"),
+                   std::tuple{static_cast<cpp_bin_float_50>(x), sin_x, cos_x});
         logger.FlushAndClear();
       });
 }

mathematica/accurate_tables.wl

@@ -0,0 +1,69 @@
+(* ::Package:: *)
+
+(* ::Title:: *)
+(*Accurate Table Printing*)
+
+
+(* ::Input:: *)
+(*SetDirectory[ParentDirectory[NotebookDirectory[]]]*)
+
+
+(* ::Input:: *)
+(*Begin["sincos`"]*)
+
+
+(* ::Input:: *)
+(*<< "functions\\sin_cos_18.wl"*)
+
+
+(* ::Input:: *)
+(*End[]*)
+
+
+(* ::Text:: *)
+(*Returns a string representation of the 100-digit approximation of the given rational, using the C++ syntax for hex literals.*)
+
+
+(* ::Input:: *)
+(*hexFloatLiteral[x_Rational,signed_]:=*)
+(* With[*)
+(*  {groups=5,*)
+(*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]],"'"]]]*)
+
+
+(* ::Text:: *)
+(*Prints the C++ header defining the accurate tables.*)
+
+
+(* ::Input:: *)
+(*Export["numerics\\accurate_tables.mathematica.h", *)
+(*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 <limits>*)
+(**)
+(*namespace principia {*)
+(*namespace numerics {*)
+(*namespace _accurate_tables {*)
+(*namespace internal {*)
+(**)
+(*struct SinCosAccurateValues {*)
+(*  double x;*)
+(*  double sin_x;*)
+(*  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"]<>"}};*)
+(**)
+(*}  // namespace internal*)
+(**)
+(*using internal::SinCosAccurateTable;*)
+(**)
+(*}  // namespace _accurate_tables*)
+(*}  // namespace numerics*)
+(*}  // namespace principia*)
+(*"], "text"]*)