From 9704d5e468834cf8ce2e3b8bdbe8b3ee3d436982 Mon Sep 17 00:00:00 2001 From: Hans Goudey Date: Tue, 9 Jan 2024 18:05:12 +0100 Subject: [PATCH] BLI: Add "numbers" math header, decouple C API Adds a header that defines the same constants as the C++ 20 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 --- .../blender/blenlib/BLI_math_angle_types.hh | 26 +++--- source/blender/blenlib/BLI_math_base.hh | 6 +- source/blender/blenlib/BLI_math_matrix.hh | 10 +-- source/blender/blenlib/BLI_math_numbers.hh | 84 +++++++++++++++++++ source/blender/blenlib/BLI_math_rotation.hh | 24 +++--- source/blender/blenlib/BLI_mesh_intersect.hh | 30 ++----- source/blender/blenlib/CMakeLists.txt | 1 + source/blender/blenlib/intern/delaunay_2d.cc | 4 +- 8 files changed, 126 insertions(+), 59 deletions(-) create mode 100644 source/blender/blenlib/BLI_math_numbers.hh diff --git a/source/blender/blenlib/BLI_math_angle_types.hh b/source/blender/blenlib/BLI_math_angle_types.hh index dbdf8c14fd4..c00a4500695 100644 --- a/source/blender/blenlib/BLI_math_angle_types.hh +++ b/source/blender/blenlib/BLI_math_angle_types.hh @@ -50,7 +50,7 @@ template struct AngleRadianBase { static AngleRadianBase from_degree(const T °rees) { - return degrees * T(M_PI / 180.0); + return degrees * T(numbers::pi / 180.0); } /** Conversions. */ @@ -64,7 +64,7 @@ template 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 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 struct AngleCartesianBase { static AngleCartesianBase from_degree(const T °rees) { - return AngleCartesianBase(degrees * T(M_PI / 180.0)); + return AngleCartesianBase(degrees * T(numbers::pi / 180.0)); } /** @@ -224,7 +224,7 @@ template 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 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 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 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 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". */ diff --git a/source/blender/blenlib/BLI_math_base.hh b/source/blender/blenlib/BLI_math_base.hh index bb7b0daf168..394985c5b19 100644 --- a/source/blender/blenlib/BLI_math_base.hh +++ b/source/blender/blenlib/BLI_math_base.hh @@ -12,7 +12,7 @@ #include #include -#include "BLI_math_base.h" +#include "BLI_math_numbers.hh" #include "BLI_utildefines.h" namespace blender::math { @@ -184,7 +184,7 @@ template inline T exp(const T &x) template 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 inline T asin(const T &a) diff --git a/source/blender/blenlib/BLI_math_matrix.hh b/source/blender/blenlib/BLI_math_matrix.hh index 36219fa6664..693047ecf62 100644 --- a/source/blender/blenlib/BLI_math_matrix.hh +++ b/source/blender/blenlib/BLI_math_matrix.hh @@ -253,7 +253,7 @@ template [[nodiscard]] MatT orthogonalize(const MatT &mat, const /** * Construct a transformation that is pivoted around the given origin point. So for instance, - * from_origin_transform(from_rotation(M_PI_2), float2(0.0f, 2.0f)) + * from_origin_transform(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 @@ -985,10 +985,10 @@ MatBase from_rotation(const QuaternionBase &rotation) { using MatT = MatBase; using DoublePrecision = typename TypeTraits::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; diff --git a/source/blender/blenlib/BLI_math_numbers.hh b/source/blender/blenlib/BLI_math_numbers.hh new file mode 100644 index 00000000000..abc7518a824 --- /dev/null +++ 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 header and will be replaced when we switch. + */ + +#include + +#include "BLI_utildefines.h" + +namespace blender::math::numbers { + +template))> +inline constexpr T e_v = 2.718281828459045235360287471352662498L; + +/* log_2 e */ +template))> +inline constexpr T log2e_v = 1.442695040888963407359924681001892137L; + +/* log_10 e */ +template))> +inline constexpr T log10e_v = 0.434294481903251827651128918916605082L; + +/* pi */ +template))> +inline constexpr T pi_v = 3.141592653589793238462643383279502884L; + +/* 1/pi */ +template))> +inline constexpr T inv_pi_v = 0.318309886183790671537767526745028724L; + +/* 1/sqrt(pi) */ +template))> +inline constexpr T inv_sqrtpi_v = 0.564189583547756286948079451560772586L; + +/* log_e 2 */ +template))> +inline constexpr T ln2_v = 0.693147180559945309417232121458176568L; + +/* log_e 10 */ +template))> +inline constexpr T ln10_v = 2.302585092994045684017991454684364208L; + +/* sqrt(2) */ +template))> +inline constexpr T sqrt2_v = 1.414213562373095048801688724209698079L; + +/* sqrt(3) */ +template))> +inline constexpr T sqrt3_v = 1.732050807568877293527446341505872367L; + +/* 1/sqrt(3) */ +template))> +inline constexpr T inv_sqrt3_v = 0.577350269189625764509148780501957456L; + +/* The Euler-Mascheroni constant */ +template))> +inline constexpr T egamma_v = 0.577215664901532860606512090082402431L; + +/* The golden ratio, (1+sqrt(5))/2 */ +template))> +inline constexpr T phi_v = 1.618033988749894848204586834365638118L; + +inline constexpr double e = e_v; +inline constexpr double log2e = log2e_v; +inline constexpr double log10e = log10e_v; +inline constexpr double pi = pi_v; +inline constexpr double inv_pi = inv_pi_v; +inline constexpr double inv_sqrtpi = inv_sqrtpi_v; +inline constexpr double ln2 = ln2_v; +inline constexpr double ln10 = ln10_v; +inline constexpr double sqrt2 = sqrt2_v; +inline constexpr double sqrt3 = sqrt3_v; +inline constexpr double inv_sqrt3 = inv_sqrt3_v; +inline constexpr double egamma = egamma_v; +inline constexpr double phi = phi_v; + +} // namespace blender::math::numbers diff --git a/source/blender/blenlib/BLI_math_rotation.hh b/source/blender/blenlib/BLI_math_rotation.hh index d83f630587b..e4e1c1a00c4 100644 --- a/source/blender/blenlib/BLI_math_rotation.hh +++ b/source/blender/blenlib/BLI_math_rotation.hh @@ -407,49 +407,49 @@ template QuaternionBase to_quaternion(const CartesianBasis &rotat case map(AxisSigned::Z_POS, AxisSigned::X_POS, AxisSigned::Y_POS): return QuaternionBase{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(M_SQRT1_2), T(0), T(0), T(-M_SQRT1_2)}; + return QuaternionBase{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(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(0), T(M_SQRT1_2), T(M_SQRT1_2), T(0)}; + return QuaternionBase{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(M_SQRT1_2), T(0), T(M_SQRT1_2), T(0)}; + return QuaternionBase{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(M_SQRT1_2), T(0), T(-M_SQRT1_2), T(0)}; + return QuaternionBase{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(0), T(0), T(1), T(0)}; case map(AxisSigned::Y_POS, AxisSigned::Z_POS, AxisSigned::X_POS): return QuaternionBase{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(0), T(0), T(M_SQRT1_2), T(M_SQRT1_2)}; + return QuaternionBase{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(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(M_SQRT1_2), T(M_SQRT1_2), T(0), T(0)}; + return QuaternionBase{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(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(M_SQRT1_2), T(0), T(0), T(M_SQRT1_2)}; + return QuaternionBase{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(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(0), T(-M_SQRT1_2), T(M_SQRT1_2), T(0)}; + return QuaternionBase{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(0), T(M_SQRT1_2), T(0), T(M_SQRT1_2)}; + return QuaternionBase{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(0), T(0), T(0), T(1)}; case map(AxisSigned::Z_NEG, AxisSigned::Y_NEG, AxisSigned::X_NEG): - return QuaternionBase{T(0), T(-M_SQRT1_2), T(0), T(M_SQRT1_2)}; + return QuaternionBase{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(0), T(1), T(0), T(0)}; case map(AxisSigned::Y_NEG, AxisSigned::Z_NEG, AxisSigned::X_POS): return QuaternionBase{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(M_SQRT1_2), T(-M_SQRT1_2), T(0), T(0)}; + return QuaternionBase{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(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(0), T(0), T(-M_SQRT1_2), T(M_SQRT1_2)}; + return QuaternionBase{T(0), T(0), T(-rcp(numbers::sqrt2)), T(rcp(numbers::sqrt2))}; } } diff --git a/source/blender/blenlib/BLI_mesh_intersect.hh b/source/blender/blenlib/BLI_mesh_intersect.hh index b331968d619..fc935bcbd21 100644 --- 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; } }; diff --git a/source/blender/blenlib/CMakeLists.txt b/source/blender/blenlib/CMakeLists.txt index 855ad5ba554..605543e7731 100644 --- 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 diff --git a/source/blender/blenlib/intern/delaunay_2d.cc b/source/blender/blenlib/intern/delaunay_2d.cc index 2e3150daebb..1fff6fd347c 100644 --- 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 *cdt_state, int maxflen = 0; for (const int f : input_faces.index_range()) { - maxflen = max_ii(maxflen, input_faces[f].size()); + maxflen = std::max(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.