Prerequisites for the implementation of the Stehlé-Zimmermann algorithm

boost_multiprecision.props

@@ -6,6 +6,7 @@
   <ItemDefinitionGroup>
     <ClCompile>
       <AdditionalIncludeDirectories>$(SolutionDir)..\Boost\multiprecision\include;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <PreprocessorDefinitions>BOOST_MP_STANDALONE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
     </ClCompile>
   </ItemDefinitionGroup>
   <ItemGroup />

functions/functions.vcxproj

@@ -4,11 +4,6 @@
     <ProjectGuid>{7cca653c-2e8f-4ffd-9e9f-bee590f3efab}</ProjectGuid>
     <RootNamespace>functions</RootNamespace>
   </PropertyGroup>
-  <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
-    <ClCompile>
-      <PreprocessorDefinitions>BOOST_MP_STANDALONE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
-    </ClCompile>
-  </ItemDefinitionGroup>
   <ItemGroup>
     <ClInclude Include="accurate_table_generator.hpp" />
     <ClInclude Include="accurate_table_generator_body.hpp" />

numerics/polynomial_in_monomial_basis.hpp

@@ -151,6 +151,7 @@ class PolynomialInMonomialBasis : public Polynomial<Value_, Argument_> {
   constexpr int degree() const override;
   bool is_zero() const override;
 
+  Coefficients const& coefficients() const;
   Argument const& origin() const;
 
   // Returns a copy of this polynomial adjusted to the given origin.
@@ -355,10 +356,10 @@ operator-(Value const& left,
 // Application monoid.
 
 template<typename LValue, typename RValue,
-         typename Argument, int ldegree_, int rdegree_>
-constexpr PolynomialInMonomialBasis<LValue, Argument, ldegree_ * rdegree_>
+         typename RArgument, int ldegree_, int rdegree_>
+constexpr PolynomialInMonomialBasis<LValue, RArgument, ldegree_ * rdegree_>
 Compose(PolynomialInMonomialBasis<LValue, RValue, ldegree_> const& left,
-        PolynomialInMonomialBasis<RValue, Argument, rdegree_> const& right);
+        PolynomialInMonomialBasis<RValue, RArgument, rdegree_> const& right);
 
 // Returns a scalar polynomial obtained by pointwise inner product of two
 // vector-valued polynomials.

numerics/polynomial_in_monomial_basis_body.hpp

@@ -12,6 +12,7 @@
 #include "absl/strings/str_cat.h"
 #include "absl/strings/str_join.h"
 #include "base/not_constructible.hpp"
+#include "boost/multiprecision/number.hpp"
 #include "geometry/cartesian_product.hpp"
 #include "geometry/serialization.hpp"
 #include "numerics/combinatorics.hpp"
@@ -24,6 +25,7 @@ namespace numerics {
 namespace _polynomial_in_monomial_basis {
 namespace internal {
 
+using namespace boost::multiprecision;
 using namespace principia::base::_not_constructible;
 using namespace principia::geometry::_cartesian_product;
 using namespace principia::geometry::_serialization;
@@ -444,6 +446,13 @@ is_zero() const {
   return coefficients_ == Coefficients{};
 }
 
+template<typename Value_, typename Argument_, int degree_>
+typename PolynomialInMonomialBasis<Value_, Argument_, degree_>::
+Coefficients const&
+PolynomialInMonomialBasis<Value_, Argument_, degree_>::coefficients() const {
+  return coefficients_;
+}
+
 template<typename Value_, typename Argument_, int degree_>
 Argument_ const& PolynomialInMonomialBasis<Value_, Argument_, degree_>::
 origin() const {
@@ -508,17 +517,19 @@ PolynomialInMonomialBasis<Value_, Argument_, degree_>::WithEvaluator() && {
 template<typename Value_, typename Argument_, int degree_>
 void PolynomialInMonomialBasis<Value_, Argument_, degree_>::
     WriteToMessage(not_null<serialization::Polynomial*> message) const {
-  message->set_degree(degree_);
-  auto* const extension =
-      message->MutableExtension(
-          serialization::PolynomialInMonomialBasis::extension);
-  TupleSerializer<Coefficients, 0>::WriteToMessage(coefficients_, extension);
-  DoubleOrQuantityOrPointOrMultivectorSerializer<
-      Argument,
-      serialization::PolynomialInMonomialBasis>::WriteToMessage(origin_,
-                                                                extension);
-  Evaluator<Value, Difference<Argument>, degree_>::WriteToMessage(
-      extension->mutable_evaluator(), evaluator_);
+  // No serialization for Boost types.
+  if constexpr (!is_number<Value>::value) {
+    message->set_degree(degree_);
+    auto* const extension = message->MutableExtension(
+        serialization::PolynomialInMonomialBasis::extension);
+    TupleSerializer<Coefficients, 0>::WriteToMessage(coefficients_, extension);
+    DoubleOrQuantityOrPointOrMultivectorSerializer<
+        Argument,
+        serialization::PolynomialInMonomialBasis>::WriteToMessage(origin_,
+                                                                  extension);
+    Evaluator<Value, Difference<Argument>, degree_>::WriteToMessage(
+        extension->mutable_evaluator(), evaluator_);
+  }
 }
 
 template<typename Value_, typename Argument_, int degree_>
@@ -750,19 +761,21 @@ operator-(Value const& left,
 }
 
 template<typename LValue, typename RValue,
-         typename Argument, int ldegree_, int rdegree_>
-constexpr PolynomialInMonomialBasis<LValue, Argument, ldegree_ * rdegree_>
+         typename RArgument, int ldegree_, int rdegree_>
+constexpr PolynomialInMonomialBasis<LValue, RArgument, ldegree_ * rdegree_>
 Compose(PolynomialInMonomialBasis<LValue, RValue, ldegree_> const& left,
-        PolynomialInMonomialBasis<RValue, Argument, rdegree_> const& right) {
+        PolynomialInMonomialBasis<RValue, RArgument, rdegree_> const& right) {
+  using LArgument = RValue;
   using LCoefficients =
-      typename PolynomialInMonomialBasis<LValue, RValue, ldegree_>::
+      typename PolynomialInMonomialBasis<LValue, LArgument, ldegree_>::
           Coefficients;
   using RCoefficients =
-      typename PolynomialInMonomialBasis<RValue, Argument, rdegree_>::
+      typename PolynomialInMonomialBasis<RValue, RArgument, rdegree_>::
           Coefficients;
-  return PolynomialInMonomialBasis<LValue, Argument, ldegree_ * rdegree_>(
+  auto const left_at_origin = left.AtOrigin(LArgument{});
+  return PolynomialInMonomialBasis<LValue, RArgument, ldegree_ * rdegree_>(
       TupleComposition<LCoefficients, RCoefficients>::Compose(
-          left.coefficients_, right.coefficients_),
+          left_at_origin.coefficients_, right.coefficients_),
       right.origin_);
 }
 

numerics/polynomial_in_monomial_basis_test.cpp

@@ -2,6 +2,8 @@
 
 #include <tuple>
 
+#include "boost/multiprecision/cpp_bin_float.hpp"
+#include "boost/multiprecision/cpp_int.hpp"
 #include "geometry/frame.hpp"
 #include "geometry/grassmann.hpp"
 #include "geometry/instant.hpp"
@@ -28,6 +30,7 @@ namespace principia {
 namespace numerics {
 
 using ::testing::Eq;
+using namespace boost::multiprecision;
 using namespace principia::geometry::_frame;
 using namespace principia::geometry::_grassmann;
 using namespace principia::geometry::_instant;
@@ -277,18 +280,14 @@ TEST_F(PolynomialInMonomialBasisTest, Affine) {
   }
 }
 
-// Compose contains a fold expression which fails to compile in Clang because of
-// https://bugs.llvm.org/show_bug.cgi?id=30590.  That bug will be fixed post-
-// 11.0.0.  Since we don't use Compose as of this writing, and working around
-// the bug would be hard, we ifdef out the test.
-#if PRINCIPIA_COMPILER_MSVC
 TEST_F(PolynomialInMonomialBasisTest, Monoid) {
-  using P0 = PolynomialInMonomialBasis<Current, Temperature, 0>;
+  using P1A = PolynomialInMonomialBasis<Current, Temperature, 1>;
   using P2A = PolynomialInMonomialBasis<Temperature, Instant, 2>;
   using P2V = PolynomialInMonomialBasis<Temperature, Time, 2>;
   using P3 = PolynomialInMonomialBasis<Current, Temperature, 3>;
   Instant const t0;
-  P0 const p0(std::tuple{9 * Ampere});
+  P1A const p1a({2 * Ampere,
+                 -4 * Ampere / Kelvin}, 3 * Kelvin);
   P2A const p2a({1 * Kelvin,
                  3 * Kelvin / Second,
                  -8 * Kelvin / Second / Second}, t0);
@@ -332,11 +331,10 @@ TEST_F(PolynomialInMonomialBasisTest, Monoid) {
     EXPECT_THAT(actual_v, AlmostEquals(-46396 * Ampere, 0));
   }
   {
-    auto const actual = Compose(p0, p2a)(t0);
-    EXPECT_THAT(actual, AlmostEquals(9 * Ampere, 0));
+    auto const actual = Compose(p1a, p2a)(t0 + 1 * Second);
+    EXPECT_THAT(actual, AlmostEquals(30 * Ampere, 0));
   }
 }
-#endif
 
 TEST_F(PolynomialInMonomialBasisTest, PointwiseInnerProduct) {
   P2V::Coefficients const coefficients({
@@ -450,6 +448,21 @@ TEST_F(PolynomialInMonomialBasisTest, EvaluateLinear) {
   EXPECT_THAT(estrin_light.EvaluateDerivative(1729 * Second), Eq(SpeedOfLight));
 }
 
+// Check that polynomials may be used with Boost multiprecision types.
+TEST_F(PolynomialInMonomialBasisTest, Boost) {
+  using P2i = PolynomialInMonomialBasis<cpp_int, cpp_int, 2>;
+  P2i const p2i({1, 3, -8});
+  EXPECT_EQ(p2i(3), -62);
+
+  using P2r = PolynomialInMonomialBasis<cpp_rational, cpp_rational, 2>;
+  P2r const p2r({1, 3, -8});
+  EXPECT_EQ(p2r(cpp_rational(1, 3)), cpp_rational(10, 9));
+
+  using P2f = PolynomialInMonomialBasis<cpp_bin_float_50, cpp_bin_float_50, 2>;
+  P2f const p2f({1, 3, -8});
+  EXPECT_EQ(p2f(cpp_bin_float_50("1.25")), cpp_bin_float_50("-7.75"));
+}
+
 // Check that polynomials may be serialized.
 TEST_F(PolynomialInMonomialBasisTest, Serialization) {
   {

quantities/concepts.hpp

@@ -15,19 +15,21 @@ using namespace base::_traits;
 using namespace quantities::_quantities;
 
 // TODO(egg): additive_group should subsume affine, but we use it there.
+// We use |convertible_to| here because we want this concept to work with
+// Boost multiprecision types which heavily use implicit conversions.
 template<typename G>
 concept additive_group = requires(G x, G y, int n) {
   G{};
-  { +x } -> std::same_as<G>;
-  { -x } -> std::same_as<G>;
-  { x + y } -> std::same_as<G>;
-  { x - y } -> std::same_as<G>;
-  { x += y } -> std::same_as<G&>;
-  { x -= y } -> std::same_as<G&>;
+  { +x } -> std::convertible_to<G>;
+  { -x } -> std::convertible_to<G>;
+  { x + y } -> std::convertible_to<G>;
+  { x - y } -> std::convertible_to<G>;
+  { x += y } -> std::convertible_to<G&>;
+  { x -= y } -> std::convertible_to<G&>;
   // An abelian group is a ℤ-module; we require the corresponding operations.
-  { n * x } -> std::same_as<G>;
-  { x * n } -> std::same_as<G>;
-  { x *= n } -> std::same_as<G&>;
+  { n * x } -> std::convertible_to<G>;
+  { x * n } -> std::convertible_to<G>;
+  { x *= n } -> std::convertible_to<G&>;
 };
 
 // A set acted upon simply transitively by an additive group.

quantities/elementary_functions.hpp

@@ -1,5 +1,6 @@
 #pragma once
 
+#include "boost/multiprecision/number.hpp"
 #include "quantities/named_quantities.hpp"
 #include "quantities/quantities.hpp"
 
@@ -8,6 +9,7 @@ namespace quantities {
 namespace _elementary_functions {
 namespace internal {
 
+using namespace boost::multiprecision;
 using namespace principia::quantities::_named_quantities;
 using namespace principia::quantities::_quantities;
 
@@ -84,6 +86,11 @@ Angle ArcTanh(double x);
 // |previous_angle|.
 Angle UnwindFrom(Angle const& previous_angle, Angle const& α);
 
+// Only dimensionless quantities can be rounded.
+template<typename Q>
+  requires is_number<Q>::value || std::floating_point<Q>
+Q Round(Q const& x);
+
 }  // namespace internal
 
 using internal::Abs;
@@ -104,6 +111,7 @@ using internal::Mod;
 using internal::NextDown;
 using internal::NextUp;
 using internal::Pow;
+using internal::Round;
 using internal::Sin;
 using internal::Sinh;
 using internal::Sqrt;

quantities/elementary_functions_body.hpp

@@ -5,30 +5,45 @@
 #include <pmmintrin.h>
 
 #include <cmath>
-#include <type_traits>
 
+#include "boost/multiprecision/cpp_int.hpp"
 #include "numerics/cbrt.hpp"
 #include "numerics/fma.hpp"
 #include "numerics/next.hpp"
+#include "quantities/concepts.hpp"
 #include "quantities/si.hpp"
 
 namespace principia {
 namespace quantities {
 namespace _elementary_functions {
 namespace internal {
 
+using namespace boost::multiprecision;
 using namespace principia::numerics::_cbrt;
 using namespace principia::numerics::_fma;
 using namespace principia::numerics::_next;
+using namespace principia::quantities::_concepts;
 using namespace principia::quantities::_si;
 
 template<typename Q1, typename Q2>
 Product<Q1, Q2> FusedMultiplyAdd(Q1 const& x,
                                  Q2 const& y,
                                  Product<Q1, Q2> const& z) {
-  return si::Unit<Product<Q1, Q2>> *
-         numerics::_fma::FusedMultiplyAdd(
-             x / si::Unit<Q1>, y / si::Unit<Q2>, z / si::Unit<Product<Q1, Q2>>);
+  if constexpr (is_number<Q1>::value || is_number<Q2>::value) {
+    if constexpr ((number_category<Q1>::value == number_kind_integer ||
+                   number_category<Q1>::value == number_kind_rational) &&
+                  (number_category<Q2>::value == number_kind_integer ||
+                   number_category<Q2>::value == number_kind_rational)) {
+      return x * y + z;
+    } else {
+      return fma(x, y, z);
+    }
+  } else {
+    return si::Unit<Product<Q1, Q2>> *
+            numerics::_fma::FusedMultiplyAdd(x / si::Unit<Q1>,
+                                             y / si::Unit<Q2>,
+                                             z / si::Unit<Product<Q1, Q2>>);
+  }
 }
 
 template<typename Q1, typename Q2>
@@ -59,8 +74,13 @@ Product<Q1, Q2> FusedNegatedMultiplySubtract(Q1 const& x,
 }
 
 template<typename Q>
-FORCE_INLINE(inline) Q Abs(Q const& quantity) {
-  return si::Unit<Q> * std::abs(quantity / si::Unit<Q>);
+FORCE_INLINE(inline)
+Q Abs(Q const& x) {
+  if constexpr (is_number<Q>::value) {
+    return abs(x);
+  } else {
+    return si::Unit<Q> * std::abs(x / si::Unit<Q>);
+  }
 }
 
 template<typename Q>
@@ -143,7 +163,16 @@ inline constexpr double Pow<3>(double x) {
 
 template<int exponent, typename Q>
 constexpr Exponentiation<Q, exponent> Pow(Q const& x) {
-  return si::Unit<Exponentiation<Q, exponent>> * Pow<exponent>(x / si::Unit<Q>);
+  if constexpr (number_category<Q>::value == number_kind_rational) {
+    // It seems that Boost does not define |pow| for |cpp_rational|.
+    return cpp_rational(pow(numerator(x), exponent),
+                        pow(denominator(x), exponent));
+  } else if constexpr (is_number<Q>::value) {
+    return pow(x, exponent);
+  } else {
+    return si::Unit<Exponentiation<Q, exponent>> *
+           Pow<exponent>(x / si::Unit<Q>);
+  }
 }
 
 inline double Sin(Angle const& α) {
@@ -175,6 +204,17 @@ Angle ArcTan(Quantity<D> const& y, Quantity<D> const& x) {
   return ArcTan(y / si::Unit<Quantity<D>>, x / si::Unit<Quantity<D>>);
 }
 
+template<typename Q>
+  requires is_number<Q>::value || std::floating_point<Q>
+Q Round(Q const& x) {
+  if constexpr (is_number<Q>::value) {
+    // TODO(phl): This is clunky.  Use |divide_qr| or something.
+    return static_cast<Q>(round(static_cast<cpp_bin_float_50>(x)));
+  } else {
+    return std::round(x);
+  }
+}
+
 inline double Sinh(Angle const& α) {
   return std::sinh((α / Radian));
 }