add BitSet::fromMask constructor, and have Printable concept also req…

…uire countDigits be implemented (special casing `double`, where we have a print method that does not use it)

include/Containers/BitSets.hpp

@@ -102,6 +102,7 @@ template <Collection T = math::Vector<uint64_t, 1>> struct BitSet {
     return unsigned(((N + usize - 1) >> ushift));
   }
   constexpr explicit BitSet(ptrdiff_t N) : data{numElementsNeeded(N), 0} {}
+  static constexpr auto fromMask(U u) -> BitSet { return BitSet{T{u}}; }
   constexpr void resizeData(ptrdiff_t N) {
     if constexpr (CanResize<T>) data.resize(N);
     else invariant(N <= std::ssize(data));

include/Math/Array.hpp

@@ -68,9 +68,50 @@ template <class T, Dimension S,
           class A = alloc::Mallocator<utils::compressed_t<T>>>
 struct ManagedArray;
 
+template <std::integral T> static constexpr auto maxPow10() -> size_t {
+  if constexpr (sizeof(T) == 1) return 3;
+  else if constexpr (sizeof(T) == 2) return 5;
+  else if constexpr (sizeof(T) == 4) return 10;
+  else if constexpr (std::signed_integral<T>) return 19;
+  else return 20;
+}
+
+template <std::unsigned_integral T> constexpr auto countDigits(T x) {
+  std::array<T, maxPow10<T>() + 1> powers;
+  powers[0] = 0;
+  powers[1] = 10;
+  for (ptrdiff_t i = 2; i < std::ssize(powers); i++)
+    powers[i] = powers[i - 1] * 10;
+  std::array<T, sizeof(T) * 8 + 1> bits;
+  if constexpr (sizeof(T) == 8) {
+    bits = {1,  1,  1,  1,  2,  2,  2,  3,  3,  3,  4,  4,  4,  4,  5,  5,  5,
+            6,  6,  6,  7,  7,  7,  7,  8,  8,  8,  9,  9,  9,  10, 10, 10, 10,
+            11, 11, 11, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16,
+            16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 19, 20};
+  } else if constexpr (sizeof(T) == 4) {
+    bits = {1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4,  5,  5, 5,
+            6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10};
+  } else if constexpr (sizeof(T) == 2) {
+    bits = {1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5};
+  } else if constexpr (sizeof(T) == 1) {
+    bits = {1, 1, 1, 1, 2, 2, 2, 3, 3};
+  }
+  T digits = bits[8 * sizeof(T) - std::countl_zero(x)];
+  return std::make_signed_t<T>(digits - (x < powers[digits - 1]));
+}
+template <std::signed_integral T> constexpr auto countDigits(T x) -> T {
+  using U = std::make_unsigned_t<T>;
+  if (x == std::numeric_limits<T>::min()) return T(sizeof(T) == 8 ? 20 : 11);
+  return countDigits<U>(U(std::abs(x))) + T{x < 0};
+}
+constexpr auto countDigits(Rational x) -> ptrdiff_t {
+  ptrdiff_t num = countDigits(x.numerator);
+  return (x.denominator == 1) ? num : num + countDigits(x.denominator) + 2;
+}
 template <typename T>
-concept Printable = requires(std::ostream &os, T x) {
+concept Printable = std::same_as<T, double> || requires(std::ostream &os, T x) {
   { os << x } -> std::convertible_to<std::ostream &>;
+  { countDigits(x) };
 };
 static_assert(Printable<int64_t>);
 void print_obj(std::ostream &os, Printable auto x) { os << x; };
@@ -1816,46 +1857,6 @@ inline auto operator<<(std::ostream &os, const T &A) -> std::ostream & {
   else B << A.t();
   return printVector(os, B);
 }
-template <std::integral T> static constexpr auto maxPow10() -> size_t {
-  if constexpr (sizeof(T) == 1) return 3;
-  else if constexpr (sizeof(T) == 2) return 5;
-  else if constexpr (sizeof(T) == 4) return 10;
-  else if constexpr (std::signed_integral<T>) return 19;
-  else return 20;
-}
-
-template <std::unsigned_integral T> constexpr auto countDigits(T x) {
-  std::array<T, maxPow10<T>() + 1> powers;
-  powers[0] = 0;
-  powers[1] = 10;
-  for (ptrdiff_t i = 2; i < std::ssize(powers); i++)
-    powers[i] = powers[i - 1] * 10;
-  std::array<T, sizeof(T) * 8 + 1> bits;
-  if constexpr (sizeof(T) == 8) {
-    bits = {1,  1,  1,  1,  2,  2,  2,  3,  3,  3,  4,  4,  4,  4,  5,  5,  5,
-            6,  6,  6,  7,  7,  7,  7,  8,  8,  8,  9,  9,  9,  10, 10, 10, 10,
-            11, 11, 11, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16,
-            16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 19, 20};
-  } else if constexpr (sizeof(T) == 4) {
-    bits = {1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4,  5,  5, 5,
-            6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10};
-  } else if constexpr (sizeof(T) == 2) {
-    bits = {1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5};
-  } else if constexpr (sizeof(T) == 1) {
-    bits = {1, 1, 1, 1, 2, 2, 2, 3, 3};
-  }
-  T digits = bits[8 * sizeof(T) - std::countl_zero(x)];
-  return std::make_signed_t<T>(digits - (x < powers[digits - 1]));
-}
-template <std::signed_integral T> constexpr auto countDigits(T x) {
-  using U = std::make_unsigned_t<T>;
-  if (x == std::numeric_limits<T>::min()) return T(sizeof(T) == 8 ? 20 : 11);
-  return countDigits<U>(U(std::abs(x))) + T{x < 0};
-}
-constexpr auto countDigits(Rational x) -> ptrdiff_t {
-  ptrdiff_t num = countDigits(x.numerator);
-  return (x.denominator == 1) ? num : num + countDigits(x.denominator) + 2;
-}
 /// \brief Returns the maximum number of digits per column of a matrix.
 constexpr auto getMaxDigits(PtrMatrix<Rational> A) -> Vector<ptrdiff_t> {
   ptrdiff_t M = ptrdiff_t(A.numRow());

include/Math/ExpDual.hpp

@@ -5,8 +5,15 @@ namespace poly::math {
 
 template <int l = 8> constexpr auto smax(auto x, auto y, auto z) {
   double m = std::max(std::max(value(x), value(y)), value(z));
-  constexpr double f = l, i = 1 / f;
-  return m + log(exp(f * (x - m)) + exp(f * (y - m)) + exp(f * (z - m))) * i;
+  static constexpr double f = l, i = 1 / f;
+  return m + i * log(exp(f * (x - m)) + exp(f * (y - m)) + exp(f * (z - m)));
+}
+template <int l = 8> constexpr auto smax(auto w, auto x, auto y, auto z) {
+  double m =
+    std::max(std::max(value(w), value(y)), std::max(value(x), value(z)));
+  static constexpr double f = l, i = 1 / f;
+  return m + i * log(exp(f * (w - m)) + exp(f * (x - m)) + exp(f * (y - m)) +
+                     exp(f * (z - m)));
 }
 
 } // namespace poly::math

test/bitset_test.cpp

@@ -84,8 +84,9 @@ TEST(FixedSizeBitSetTest, BasicAssertions) {
 }
 // NOLINTNEXTLINE(modernize-use-trailing-return-type)
 TEST(FixedSizeSmallBitSetTest, BasicAssertions) {
-  static_assert(sizeof(BitSet<std::array<uint16_t, 1>>) == 2);
-  BitSet<std::array<uint16_t, 1>> bs;
+  using SB = BitSet<std::array<uint16_t, 1>>;
+  static_assert(sizeof(SB) == 2);
+  SB bs;
   bs[4] = true;
   bs[10] = true;
   bs[7] = true;
@@ -98,4 +99,5 @@ TEST(FixedSizeSmallBitSetTest, BasicAssertions) {
   EXPECT_EQ(sv[1], 5);
   EXPECT_EQ(sv[2], 7);
   EXPECT_EQ(sv[3], 10);
+  EXPECT_EQ(SB::fromMask(1200).data[0], 1200);
 }