Rationalize the handling of Boost multiprecision types

geometry/interval_body.hpp

@@ -26,12 +26,7 @@ bool Interval<T>::empty() const {
 
 template<typename T>
 T Interval<T>::midpoint() const {
-  if constexpr (is_number<T>::value) {
-    DCHECK_GE(max, min);
-    return min + measure() / 2;
-  } else {
-    return max >= min ? min + measure() / 2 : min + NaN<Difference<T>>;
-  }
+  return max >= min ? min + measure() / 2 : min + NaN<Difference<T>>;
 }
 
 template<typename T>

numerics/elementary_functions.hpp

@@ -1,6 +1,6 @@
 #pragma once
 
-#include "boost/multiprecision/number.hpp"
+#include "quantities/dimensions.hpp"
 #include "quantities/named_quantities.hpp"
 #include "quantities/quantities.hpp"
 
@@ -9,7 +9,7 @@ namespace numerics {
 namespace _elementary_functions {
 namespace internal {
 
-using namespace boost::multiprecision;
+using namespace principia::quantities::_dimensions;
 using namespace principia::quantities::_named_quantities;
 using namespace principia::quantities::_quantities;
 
@@ -81,9 +81,7 @@ 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>
+template<dimensionless Q>
 Q Round(Q const& x);
 
 }  // namespace internal

numerics/elementary_functions_body.hpp

@@ -2,82 +2,124 @@
 
 #include "numerics/elementary_functions.hpp"
 
-#include <pmmintrin.h>
-
 #include <cmath>
 
 #include "boost/multiprecision/cpp_bin_float.hpp"
-#include "boost/multiprecision/cpp_int.hpp"
 #include "numerics/cbrt.hpp"
 #include "numerics/fma.hpp"
 #include "numerics/next.hpp"
+#include "quantities/cantor.hpp"
+#include "quantities/concepts.hpp"
 #include "quantities/si.hpp"
 
 namespace principia {
 namespace numerics {
 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::_cantor;
+using namespace principia::quantities::_concepts;
 using namespace principia::quantities::_si;
 
+// For these types there is no concern about performance or accuracy, but
+// `FusedMultiplyAdd` and friends need to exist.
+namespace noncritical {
+
+template<boost_cpp_bin_float Q1, boost_cpp_bin_float Q2>
+Product<Q1, Q2> FusedMultiplyAdd(Q1 const& x,
+                                 Q2 const& y,
+                                 Product<Q1, Q2> const& z) {
+  return fma(x, y, z);
+}
+
+template<countable Q1, countable Q2>
+Product<Q1, Q2> FusedMultiplyAdd(Q1 const& x,
+                                 Q2 const& y,
+                                 Product<Q1, Q2> const& z) {
+  return x * y + z;
+}
+
+template<boost_cpp_number Q>
+Q Abs(Q const& x) {
+  return abs(x);
+}
+
+template<boost_cpp_number Q>
+Q Round(Q const& x) {
+  // TODO(phl): This is clunky.  Use `divide_qr` or something.
+  return static_cast<Q>(round(static_cast<cpp_bin_float_50>(x)));
+}
+
+}  // namespace noncritical
+
 template<typename Q1, typename Q2>
 Product<Q1, Q2> FusedMultiplyAdd(Q1 const& x,
                                  Q2 const& y,
                                  Product<Q1, Q2> const& z) {
-  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 {
+  if constexpr (convertible_to_quantity<Q1> && convertible_to_quantity<Q2>) {
     return si::Unit<Product<Q1, Q2>> *
-            numerics::_fma::FusedMultiplyAdd(x / si::Unit<Q1>,
-                                             y / si::Unit<Q2>,
-                                             z / si::Unit<Product<Q1, Q2>>);
+           numerics::_fma::FusedMultiplyAdd(x / si::Unit<Q1>,
+                                            y / si::Unit<Q2>,
+                                            z / si::Unit<Product<Q1, Q2>>);
+  } else {
+    return noncritical::FusedMultiplyAdd(x, y, z);
   }
 }
 
 template<typename Q1, typename Q2>
 Product<Q1, Q2> FusedMultiplySubtract(Q1 const& x,
                                       Q2 const& y,
                                       Product<Q1, Q2> const& z) {
-  return si::Unit<Product<Q1, Q2>> *
-         numerics::_fma::FusedMultiplySubtract(
-             x / si::Unit<Q1>, y / si::Unit<Q2>, z / si::Unit<Product<Q1, Q2>>);
+  if constexpr (convertible_to_quantity<Q1> && convertible_to_quantity<Q2>) {
+    return si::Unit<Product<Q1, Q2>> *
+           numerics::_fma::FusedMultiplySubtract(x / si::Unit<Q1>,
+                                                 y / si::Unit<Q2>,
+                                                 z / si::Unit<Product<Q1, Q2>>);
+  } else {
+    return noncritical::FusedMultiplyAdd(x, y, -z);
+  }
 }
 
 template<typename Q1, typename Q2>
 Product<Q1, Q2> FusedNegatedMultiplyAdd(Q1 const& x,
                                         Q2 const& y,
                                         Product<Q1, Q2> const& z) {
-  return si::Unit<Product<Q1, Q2>> *
-         numerics::_fma::FusedNegatedMultiplyAdd(
-             x / si::Unit<Q1>, y / si::Unit<Q2>, z / si::Unit<Product<Q1, Q2>>);
+  if constexpr (convertible_to_quantity<Q1> && convertible_to_quantity<Q2>) {
+    return si::Unit<Product<Q1, Q2>> * numerics::_fma::FusedNegatedMultiplyAdd(
+                                           x / si::Unit<Q1>,
+                                           y / si::Unit<Q2>,
+                                           z / si::Unit<Product<Q1, Q2>>);
+  } else {
+    return noncritical::FusedMultiplyAdd(-x, y, z);
+  }
 }
 
 template<typename Q1, typename Q2>
 Product<Q1, Q2> FusedNegatedMultiplySubtract(Q1 const& x,
                                              Q2 const& y,
                                              Product<Q1, Q2> const& z) {
-  return si::Unit<Product<Q1, Q2>> *
-         numerics::_fma::FusedNegatedMultiplySubtract(
-             x / si::Unit<Q1>, y / si::Unit<Q2>, z / si::Unit<Product<Q1, Q2>>);
+  if constexpr (convertible_to_quantity<Q1> && convertible_to_quantity<Q2>) {
+    return si::Unit<Product<Q1, Q2>> *
+           numerics::_fma::FusedNegatedMultiplySubtract(
+               x / si::Unit<Q1>,
+               y / si::Unit<Q2>,
+               z / si::Unit<Product<Q1, Q2>>);
+  } else {
+    return noncritical::FusedMultiplyAdd(-x, y, -z);
+  }
 }
 
 template<typename Q>
 FORCE_INLINE(inline)
 Q Abs(Q const& x) {
-  if constexpr (is_number<Q>::value) {
-    return abs(x);
-  } else {
+  if constexpr (convertible_to_quantity<Q>) {
     return si::Unit<Q> * std::abs(x / si::Unit<Q>);
+  } else {
+    return noncritical::Abs(x);
   }
 }
 
@@ -146,14 +188,14 @@ 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>
+template<dimensionless 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 {
+  if constexpr (std::floating_point<Q>) {
     return std::round(x);
+  } else if constexpr (std::integral<Q>) {
+    return x;
+  } else {
+    return noncritical::Round(x);
   }
 }
 

numerics/elementary_functions_test.cpp

@@ -1,6 +1,7 @@
 #include "numerics/elementary_functions.hpp"
 
 #include "base/cpuid.hpp"
+#include "boost/multiprecision/cpp_int.hpp"
 #include "gtest/gtest.h"
 #include "numerics/fma.hpp"
 #include "quantities/astronomy.hpp"
@@ -17,6 +18,7 @@ namespace quantities {
 
 using ::testing::Eq;
 using ::testing::Lt;
+using namespace boost::multiprecision;
 using namespace principia::base::_cpuid;
 using namespace principia::numerics::_elementary_functions;
 using namespace principia::numerics::_fma;
@@ -32,6 +34,15 @@ using namespace principia::testing_utilities::_vanishes_before;
 class ElementaryFunctionsTest : public testing::Test {};
 
 TEST_F(ElementaryFunctionsTest, FMA) {
+  EXPECT_EQ(cpp_int(11), FusedMultiplyAdd(cpp_int(2), cpp_int(3), cpp_int(5)));
+  EXPECT_EQ(cpp_rational(11, 2),
+            FusedMultiplyAdd(
+                cpp_rational(2, 1), cpp_rational(3, 2), cpp_rational(5, 2)));
+  EXPECT_EQ(cpp_bin_float_50("11.0"),
+            FusedMultiplyAdd(cpp_bin_float_50("2.0"),
+                             cpp_bin_float_50("3.0"),
+                             cpp_bin_float_50("5.0")));
+
   if (!CanEmitFMAInstructions || !CPUIDFeatureFlag::FMA.IsSet()) {
     GTEST_SKIP() << "Cannot test FMA on a machine without FMA";
   }

numerics/polynomial_in_monomial_basis_body.hpp

@@ -12,26 +12,26 @@
 #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"
 #include "numerics/elementary_functions.hpp"
 #include "numerics/quadrature.hpp"
+#include "quantities/cantor.hpp"
 #include "quantities/quantities.hpp"
 
 namespace principia {
 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;
 using namespace principia::numerics::_combinatorics;
 using namespace principia::numerics::_elementary_functions;
 using namespace principia::numerics::_quadrature;
+using namespace principia::quantities::_cantor;
 using namespace principia::quantities::_quantities;
 
 // A helper for changing the origin of a monomial (x - x₁)ⁿ.  It computes the
@@ -518,7 +518,7 @@ template<typename Value_, typename Argument_, int degree_>
 void PolynomialInMonomialBasis<Value_, Argument_, degree_>::
     WriteToMessage(not_null<serialization::Polynomial*> message) const {
   // No serialization for Boost types.
-  if constexpr (!is_number<Value>::value) {
+  if constexpr (!boost_cpp_number<Value>) {
     message->set_degree(degree_);
     auto* const extension = message->MutableExtension(
         serialization::PolynomialInMonomialBasis::extension);

quantities/cantor.hpp

@@ -0,0 +1,59 @@
+#pragma once
+
+#include <concepts>
+
+#include "base/traits.hpp"
+#include "boost/multiprecision/number.hpp"
+
+namespace principia {
+namespace quantities {
+namespace _cantor {
+namespace internal {
+
+using namespace boost::multiprecision;
+using namespace principia::base::_traits;
+
+// The `boost_cpp_` concepts should be used sparingly, and only in places where
+// the Boost multiprecision API differs from ours or from the standard C++ API.
+// At any rate, the multiprecision traits should never be used directly.
+
+template<typename T>
+concept boost_cpp_number =
+    (is_number<T>::value || is_number_expression<T>::value) &&
+    number_category<T>::value != number_kind_unknown;
+
+template<typename T>
+concept boost_cpp_bin_float =
+    boost_cpp_number<T> &&
+    number_category<T>::value == number_kind_floating_point;
+
+template<typename T>
+concept boost_cpp_int =
+    boost_cpp_number<T> && number_category<T>::value == number_kind_integer;
+
+template<typename T>
+concept boost_cpp_rational =
+    boost_cpp_number<T> && number_category<T>::value == number_kind_rational;
+
+template<typename T>
+concept discrete = std::integral<T> || boost_cpp_int<T>;
+
+template<typename T>
+concept countable = discrete<T> || boost_cpp_rational<T>;
+
+template<typename T>
+concept continuum = std::floating_point<T> || boost_cpp_bin_float<T>;
+
+}  // namespace internal
+
+using internal::boost_cpp_bin_float;
+using internal::boost_cpp_int;
+using internal::boost_cpp_number;
+using internal::boost_cpp_rational;
+using internal::continuum;
+using internal::countable;
+using internal::discrete;
+
+}  // namespace _cantor
+}  // namespace quantities
+}  // namespace principia

quantities/cantor_test.cpp

@@ -0,0 +1,47 @@
+#include "quantities/cantor.hpp"
+
+#include "boost/multiprecision/cpp_bin_float.hpp"
+#include "boost/multiprecision/cpp_int.hpp"
+#include "gtest/gtest.h"
+
+namespace principia {
+namespace quantities {
+
+using namespace boost::multiprecision;
+using namespace principia::quantities::_cantor;
+
+TEST(Cantor, CppBinFloat) {
+  static_assert(boost_cpp_bin_float<cpp_bin_float_50>);
+  static_assert(boost_cpp_bin_float<cpp_bin_float_100>);
+  static_assert(!boost_cpp_bin_float<cpp_int>);
+  static_assert(!boost_cpp_bin_float<cpp_rational>);
+  static_assert(!boost_cpp_bin_float<double>);
+  static_assert(!boost_cpp_bin_float<int>);
+}
+
+TEST(Cantor, Discrete) {
+  static_assert(discrete<int>);
+  static_assert(discrete<cpp_int>);
+  static_assert(!discrete<cpp_rational>);
+  static_assert(!discrete<double>);
+  static_assert(!discrete<cpp_bin_float_50>);
+}
+
+TEST(Cantor, Countable) {
+  static_assert(countable<int>);
+  static_assert(countable<cpp_int>);
+  static_assert(countable<cpp_rational>);
+  static_assert(!countable<double>);
+  static_assert(!countable<cpp_bin_float_50>);
+}
+
+TEST(Cantor, Continuum) {
+  static_assert(continuum<double>);
+  static_assert(continuum<cpp_bin_float_50>);
+  static_assert(!continuum<int>);
+  static_assert(!continuum<cpp_int>);
+  static_assert(!continuum<cpp_rational>);
+}
+
+}  // namespace quantities
+}  // namespace principia