BLI: Add "numbers" math header, decouple C API

Adds a header that defines the same constants as the C++ 20 <numbers> header. Benefits: - Decouple our C++ and C math APIs - Avoid using macros everywhere, nicer syntax - Less header parsing during compilation - Can be replaced by `std::numbers` with C++ 20 Downsides: - There are fewer numbers defined in the C++ standard header - Maybe we should just wait until we can use C++ 20 Pull Request: https://projects.blender.org/blender/blender/pulls/116805
2024-01-09 18:05:12 +01:00
parent 96f93ff42e
commit 9704d5e468
8 changed files with 126 additions and 59 deletions
--- a/source/blender/blenlib/BLI_math_angle_types.hh
+++ b/source/blender/blenlib/BLI_math_angle_types.hh
@@ -50,7 +50,7 @@ template<typename T> struct AngleRadianBase {

  static AngleRadianBase from_degree(const T &degrees)
  {
-    return degrees * T(M_PI / 180.0);
+    return degrees * T(numbers::pi / 180.0);
  }

  /** Conversions. */
@@ -64,7 +64,7 @@ template<typename T> struct AngleRadianBase {
  /* Return angle value in degree. */
  T degree() const
  {
-    return value_ * T(180.0 / M_PI);
+    return value_ * T(180.0 / numbers::pi);
  }

  /* Return angle value in radian. */
@@ -80,7 +80,7 @@ template<typename T> struct AngleRadianBase {
   */
  AngleRadianBase wrapped() const
  {
-    return math::mod_periodic(value_ + T(M_PI), T(2 * M_PI)) - T(M_PI);
+    return math::mod_periodic(value_ + T(numbers::pi), T(2 * numbers::pi)) - T(numbers::pi);
  }

  /**
@@ -197,7 +197,7 @@ template<typename T> struct AngleCartesianBase {

  static AngleCartesianBase from_degree(const T &degrees)
  {
-    return AngleCartesianBase(degrees * T(M_PI / 180.0));
+    return AngleCartesianBase(degrees * T(numbers::pi / 180.0));
  }

  /**
@@ -224,7 +224,7 @@ template<typename T> struct AngleCartesianBase {
  /* Return angle value in degree. */
  T degree() const
  {
-    return T(*this) * T(180.0 / M_PI);
+    return T(*this) * T(180.0 / numbers::pi);
  }

  /* Return angle value in radian. */
@@ -420,19 +420,19 @@ template<typename T = float> struct AngleFraction {
    const bool is_negative = numerator_ < 0;
    /* TODO jump table. */
    if (abs(numerator_) == denominator_ * 2) {
-      return is_negative ? T(-M_PI * 2) : T(M_PI * 2);
+      return is_negative ? T(-numbers::pi * 2) : T(numbers::pi * 2);
    }
    if (abs(numerator_) == denominator_) {
-      return is_negative ? T(-M_PI) : T(M_PI);
+      return is_negative ? T(-numbers::pi) : T(numbers::pi);
    }
    if (numerator_ == 0) {
      return T(0);
    }
    if (abs(numerator_) * 2 == denominator_) {
-      return is_negative ? T(-M_PI_2) : T(M_PI_2);
+      return is_negative ? T(-numbers::pi * 0.5) : T(numbers::pi * 0.5);
    }
    if (abs(numerator_) * 4 == denominator_) {
-      return is_negative ? T(-M_PI_4) : T(M_PI_4);
+      return is_negative ? T(-numbers::pi * 0.25) : T(numbers::pi * 0.25);
    }
    /* TODO(fclem): No idea if this is precise or not. Just doing something for now. */
    const int64_t number_of_pi = numerator_ / denominator_;
@@ -442,14 +442,14 @@ template<typename T = float> struct AngleFraction {
    /* TODO(fclem): This is conservative. Could find a better threshold. */
    if (slice_numerator > 0xFFFFFFFF || denominator_ > 0xFFFFFFFF) {
      /* Certainly loose precision. */
-      slice_of_pi = T(M_PI) * slice_numerator / T(denominator_);
+      slice_of_pi = T(numbers::pi) * slice_numerator / T(denominator_);
    }
    else {
      /* Pi as a fraction can be expressed as 80143857 / 25510582 with 15th digit of precision. */
      slice_of_pi = T(slice_numerator * 80143857) / T(denominator_ * 25510582);
    }
    /* If angle is inside [-pi..pi] range, `number_of_pi` is 0 and has no effect on precision. */
-    return slice_of_pi + T(M_PI) * number_of_pi;
+    return slice_of_pi + T(numbers::pi) * number_of_pi;
  }

  /** Methods. */
@@ -631,7 +631,7 @@ template<typename T = float> struct AngleFraction {
          break;
        case 1:
        case 3:
-          x = y = T(M_SQRT1_2);
+          x = y = math::rcp(T(numbers::sqrt2));
          break;
        default:
          BLI_assert_unreachable();
@@ -661,7 +661,7 @@ template<typename T = float> struct AngleFraction {
      }
      /* Resulting angle should be oscillating in [0..pi/4] range. */
      BLI_assert(a.numerator_ >= 0 && a.numerator_ <= a.denominator_ / 4);
-      T angle = T(M_PI) * (T(a.numerator_) / T(a.denominator_));
+      T angle = T(numbers::pi) * (T(a.numerator_) / T(a.denominator_));
      x = math::cos(angle);
      y = math::sin(angle);
      /* Diagonal symmetry "unfolding". */
--- a/source/blender/blenlib/BLI_math_base.hh
+++ b/source/blender/blenlib/BLI_math_base.hh
@@ -12,7 +12,7 @@
 #include <cmath>
 #include <type_traits>

-#include "BLI_math_base.h"
+#include "BLI_math_numbers.hh"
 #include "BLI_utildefines.h"

 namespace blender::math {
@@ -184,7 +184,7 @@ template<typename T> inline T exp(const T &x)
 template<typename T> inline T safe_acos(const T &a)
 {
  if (UNLIKELY(a <= T(-1))) {
-    return T(M_PI);
+    return T(numbers::pi);
  }
  else if (UNLIKELY(a >= T(1))) {
    return T(0);
@@ -207,7 +207,7 @@ inline float safe_acos_approx(float x)
   */
  const float a = std::sqrt(1.0f - m) *
                  (1.5707963267f + m * (-0.213300989f + m * (0.077980478f + m * -0.02164095f)));
-  return x < 0.0f ? float(M_PI) - a : a;
+  return x < 0.0f ? float(numbers::pi) - a : a;
 }

 template<typename T> inline T asin(const T &a)
--- a/source/blender/blenlib/BLI_math_matrix.hh
+++ b/source/blender/blenlib/BLI_math_matrix.hh
@@ -253,7 +253,7 @@ template<typename MatT> [[nodiscard]] MatT orthogonalize(const MatT &mat, const

 /**
 * Construct a transformation that is pivoted around the given origin point. So for instance,
- * from_origin_transform<MatT>(from_rotation(M_PI_2), float2(0.0f, 2.0f))
+ * from_origin_transform<MatT>(from_rotation(numbers::pi * 0.5), float2(0.0f, 2.0f))
 * will construct a transformation representing a 90 degree rotation around the point (0, 2).
 */
 template<typename MatT, typename VectorT>
@@ -985,10 +985,10 @@ MatBase<T, NumCol, NumRow> from_rotation(const QuaternionBase<T> &rotation)
 {
  using MatT = MatBase<T, NumCol, NumRow>;
  using DoublePrecision = typename TypeTraits<T>::DoublePrecision;
-  const DoublePrecision q0 = M_SQRT2 * DoublePrecision(rotation.w);
-  const DoublePrecision q1 = M_SQRT2 * DoublePrecision(rotation.x);
-  const DoublePrecision q2 = M_SQRT2 * DoublePrecision(rotation.y);
-  const DoublePrecision q3 = M_SQRT2 * DoublePrecision(rotation.z);
+  const DoublePrecision q0 = numbers::sqrt2 * DoublePrecision(rotation.w);
+  const DoublePrecision q1 = numbers::sqrt2 * DoublePrecision(rotation.x);
+  const DoublePrecision q2 = numbers::sqrt2 * DoublePrecision(rotation.y);
+  const DoublePrecision q3 = numbers::sqrt2 * DoublePrecision(rotation.z);

  const DoublePrecision qda = q0 * q1;
  const DoublePrecision qdb = q0 * q2;
--- a/source/blender/blenlib/BLI_math_numbers.hh
+++ b/source/blender/blenlib/BLI_math_numbers.hh
@@ -0,0 +1,84 @@
+/* SPDX-FileCopyrightText: 2024 Blender Authors
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ *
+ * Contains the same definitions as the C++20 <numbers> header and will be replaced when we switch.
+ */
+
+#include <type_traits>
+
+#include "BLI_utildefines.h"
+
+namespace blender::math::numbers {
+
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T e_v = 2.718281828459045235360287471352662498L;
+
+/* log_2 e */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T log2e_v = 1.442695040888963407359924681001892137L;
+
+/* log_10 e */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T log10e_v = 0.434294481903251827651128918916605082L;
+
+/* pi */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T pi_v = 3.141592653589793238462643383279502884L;
+
+/* 1/pi */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T inv_pi_v = 0.318309886183790671537767526745028724L;
+
+/* 1/sqrt(pi) */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T inv_sqrtpi_v = 0.564189583547756286948079451560772586L;
+
+/* log_e 2 */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T ln2_v = 0.693147180559945309417232121458176568L;
+
+/* log_e 10 */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T ln10_v = 2.302585092994045684017991454684364208L;
+
+/* sqrt(2) */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T sqrt2_v = 1.414213562373095048801688724209698079L;
+
+/* sqrt(3) */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T sqrt3_v = 1.732050807568877293527446341505872367L;
+
+/* 1/sqrt(3) */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T inv_sqrt3_v = 0.577350269189625764509148780501957456L;
+
+/* The Euler-Mascheroni constant */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T egamma_v = 0.577215664901532860606512090082402431L;
+
+/* The golden ratio, (1+sqrt(5))/2 */
+template<typename T, BLI_ENABLE_IF((std::is_floating_point_v<T>))>
+inline constexpr T phi_v = 1.618033988749894848204586834365638118L;
+
+inline constexpr double e = e_v<double>;
+inline constexpr double log2e = log2e_v<double>;
+inline constexpr double log10e = log10e_v<double>;
+inline constexpr double pi = pi_v<double>;
+inline constexpr double inv_pi = inv_pi_v<double>;
+inline constexpr double inv_sqrtpi = inv_sqrtpi_v<double>;
+inline constexpr double ln2 = ln2_v<double>;
+inline constexpr double ln10 = ln10_v<double>;
+inline constexpr double sqrt2 = sqrt2_v<double>;
+inline constexpr double sqrt3 = sqrt3_v<double>;
+inline constexpr double inv_sqrt3 = inv_sqrt3_v<double>;
+inline constexpr double egamma = egamma_v<double>;
+inline constexpr double phi = phi_v<double>;
+
+}  // namespace blender::math::numbers
--- a/source/blender/blenlib/BLI_math_rotation.hh
+++ b/source/blender/blenlib/BLI_math_rotation.hh
@@ -407,49 +407,49 @@ template<typename T> QuaternionBase<T> to_quaternion(const CartesianBasis &rotat
    case map(AxisSigned::Z_POS, AxisSigned::X_POS, AxisSigned::Y_POS):
      return QuaternionBase<T>{T(0.5), T(-0.5), T(-0.5), T(-0.5)};
    case map(AxisSigned::Y_NEG, AxisSigned::X_POS, AxisSigned::Z_POS):
-      return QuaternionBase<T>{T(M_SQRT1_2), T(0), T(0), T(-M_SQRT1_2)};
+      return QuaternionBase<T>{T(rcp(numbers::sqrt2)), T(0), T(0), T(-rcp(numbers::sqrt2))};
    case map(AxisSigned::Z_NEG, AxisSigned::X_POS, AxisSigned::Y_NEG):
      return QuaternionBase<T>{T(0.5), T(0.5), T(0.5), T(-0.5)};
    case map(AxisSigned::Y_POS, AxisSigned::X_POS, AxisSigned::Z_NEG):
-      return QuaternionBase<T>{T(0), T(M_SQRT1_2), T(M_SQRT1_2), T(0)};
+      return QuaternionBase<T>{T(0), T(rcp(numbers::sqrt2)), T(rcp(numbers::sqrt2)), T(0)};
    case map(AxisSigned::Z_NEG, AxisSigned::Y_POS, AxisSigned::X_POS):
-      return QuaternionBase<T>{T(M_SQRT1_2), T(0), T(M_SQRT1_2), T(0)};
+      return QuaternionBase<T>{T(rcp(numbers::sqrt2)), T(0), T(rcp(numbers::sqrt2)), T(0)};
    case map(AxisSigned::Z_POS, AxisSigned::Y_POS, AxisSigned::X_NEG):
-      return QuaternionBase<T>{T(M_SQRT1_2), T(0), T(-M_SQRT1_2), T(0)};
+      return QuaternionBase<T>{T(rcp(numbers::sqrt2)), T(0), T(-rcp(numbers::sqrt2)), T(0)};
    case map(AxisSigned::X_NEG, AxisSigned::Y_POS, AxisSigned::Z_NEG):
      return QuaternionBase<T>{T(0), T(0), T(1), T(0)};
    case map(AxisSigned::Y_POS, AxisSigned::Z_POS, AxisSigned::X_POS):
      return QuaternionBase<T>{T(0.5), T(0.5), T(0.5), T(0.5)};
    case map(AxisSigned::X_NEG, AxisSigned::Z_POS, AxisSigned::Y_POS):
-      return QuaternionBase<T>{T(0), T(0), T(M_SQRT1_2), T(M_SQRT1_2)};
+      return QuaternionBase<T>{T(0), T(0), T(rcp(numbers::sqrt2)), T(rcp(numbers::sqrt2))};
    case map(AxisSigned::Y_NEG, AxisSigned::Z_POS, AxisSigned::X_NEG):
      return QuaternionBase<T>{T(0.5), T(0.5), T(-0.5), T(-0.5)};
    case map(AxisSigned::X_POS, AxisSigned::Z_POS, AxisSigned::Y_NEG):
-      return QuaternionBase<T>{T(M_SQRT1_2), T(M_SQRT1_2), T(0), T(0)};
+      return QuaternionBase<T>{T(rcp(numbers::sqrt2)), T(rcp(numbers::sqrt2)), T(0), T(0)};
    case map(AxisSigned::Z_NEG, AxisSigned::X_NEG, AxisSigned::Y_POS):
      return QuaternionBase<T>{T(0.5), T(-0.5), T(0.5), T(0.5)};
    case map(AxisSigned::Y_POS, AxisSigned::X_NEG, AxisSigned::Z_POS):
-      return QuaternionBase<T>{T(M_SQRT1_2), T(0), T(0), T(M_SQRT1_2)};
+      return QuaternionBase<T>{T(rcp(numbers::sqrt2)), T(0), T(0), T(rcp(numbers::sqrt2))};
    case map(AxisSigned::Z_POS, AxisSigned::X_NEG, AxisSigned::Y_NEG):
      return QuaternionBase<T>{T(0.5), T(0.5), T(-0.5), T(0.5)};
    case map(AxisSigned::Y_NEG, AxisSigned::X_NEG, AxisSigned::Z_NEG):
-      return QuaternionBase<T>{T(0), T(-M_SQRT1_2), T(M_SQRT1_2), T(0)};
+      return QuaternionBase<T>{T(0), T(-rcp(numbers::sqrt2)), T(rcp(numbers::sqrt2)), T(0)};
    case map(AxisSigned::Z_POS, AxisSigned::Y_NEG, AxisSigned::X_POS):
-      return QuaternionBase<T>{T(0), T(M_SQRT1_2), T(0), T(M_SQRT1_2)};
+      return QuaternionBase<T>{T(0), T(rcp(numbers::sqrt2)), T(0), T(rcp(numbers::sqrt2))};
    case map(AxisSigned::X_NEG, AxisSigned::Y_NEG, AxisSigned::Z_POS):
      return QuaternionBase<T>{T(0), T(0), T(0), T(1)};
    case map(AxisSigned::Z_NEG, AxisSigned::Y_NEG, AxisSigned::X_NEG):
-      return QuaternionBase<T>{T(0), T(-M_SQRT1_2), T(0), T(M_SQRT1_2)};
+      return QuaternionBase<T>{T(0), T(-rcp(numbers::sqrt2)), T(0), T(rcp(numbers::sqrt2))};
    case map(AxisSigned::X_POS, AxisSigned::Y_NEG, AxisSigned::Z_NEG):
      return QuaternionBase<T>{T(0), T(1), T(0), T(0)};
    case map(AxisSigned::Y_NEG, AxisSigned::Z_NEG, AxisSigned::X_POS):
      return QuaternionBase<T>{T(0.5), T(-0.5), T(0.5), T(-0.5)};
    case map(AxisSigned::X_POS, AxisSigned::Z_NEG, AxisSigned::Y_POS):
-      return QuaternionBase<T>{T(M_SQRT1_2), T(-M_SQRT1_2), T(0), T(0)};
+      return QuaternionBase<T>{T(rcp(numbers::sqrt2)), T(-rcp(numbers::sqrt2)), T(0), T(0)};
    case map(AxisSigned::Y_POS, AxisSigned::Z_NEG, AxisSigned::X_NEG):
      return QuaternionBase<T>{T(0.5), T(-0.5), T(-0.5), T(0.5)};
    case map(AxisSigned::X_NEG, AxisSigned::Z_NEG, AxisSigned::Y_NEG):
-      return QuaternionBase<T>{T(0), T(0), T(-M_SQRT1_2), T(M_SQRT1_2)};
+      return QuaternionBase<T>{T(0), T(0), T(-rcp(numbers::sqrt2)), T(rcp(numbers::sqrt2))};
  }
 }

--- a/source/blender/blenlib/BLI_mesh_intersect.hh
+++ b/source/blender/blenlib/BLI_mesh_intersect.hh
@@ -342,42 +342,24 @@ struct BoundingBox {

  void combine(const float3 &p)
  {
-    min.x = min_ff(min.x, p.x);
-    min.y = min_ff(min.y, p.y);
-    min.z = min_ff(min.z, p.z);
-    max.x = max_ff(max.x, p.x);
-    max.y = max_ff(max.y, p.y);
-    max.z = max_ff(max.z, p.z);
+    math::min_max(p, this->min, this->max);
  }

  void combine(const double3 &p)
  {
-    min.x = min_ff(min.x, float(p.x));
-    min.y = min_ff(min.y, float(p.y));
-    min.z = min_ff(min.z, float(p.z));
-    max.x = max_ff(max.x, float(p.x));
-    max.y = max_ff(max.y, float(p.y));
-    max.z = max_ff(max.z, float(p.z));
+    math::min_max(float3(p), this->min, this->max);
  }

  void combine(const BoundingBox &bb)
  {
-    min.x = min_ff(min.x, bb.min.x);
-    min.y = min_ff(min.y, bb.min.y);
-    min.z = min_ff(min.z, bb.min.z);
-    max.x = max_ff(max.x, bb.max.x);
-    max.y = max_ff(max.y, bb.max.y);
-    max.z = max_ff(max.z, bb.max.z);
+    min = math::min(this->min, bb.min);
+    max = math::max(this->max, bb.max);
  }

  void expand(float pad)
  {
-    min.x -= pad;
-    min.y -= pad;
-    min.z -= pad;
-    max.x += pad;
-    max.y += pad;
-    max.z += pad;
+    min -= pad;
+    max += pad;
  }
 };

--- a/source/blender/blenlib/CMakeLists.txt
+++ b/source/blender/blenlib/CMakeLists.txt
@@ -292,6 +292,7 @@ set(SRC
  BLI_math_matrix.hh
  BLI_math_matrix_types.hh
  BLI_math_mpq.hh
+  BLI_math_numbers.hh
  BLI_math_quaternion.hh
  BLI_math_quaternion_types.hh
  BLI_math_rotation.h
--- a/source/blender/blenlib/intern/delaunay_2d.cc
+++ b/source/blender/blenlib/intern/delaunay_2d.cc
@@ -2214,11 +2214,11 @@ int add_face_constraints(CDT_state<T> *cdt_state,

  int maxflen = 0;
  for (const int f : input_faces.index_range()) {
-    maxflen = max_ii(maxflen, input_faces[f].size());
+    maxflen = std::max<int>(maxflen, input_faces[f].size());
  }
  /* For convenience in debugging, make face_edge_offset be a power of 10. */
  cdt_state->face_edge_offset = power_of_10_greater_equal_to(
-      max_ii(maxflen, cdt_state->face_edge_offset));
+      std::max(maxflen, cdt_state->face_edge_offset));
  /* The original_edge encoding scheme doesn't work if the following is false.
   * If we really have that many faces and that large a max face length that when multiplied
   * together the are >= INT_MAX, then the Delaunay calculation will take unreasonably long anyway.