Using bit_cast, and making more things constexpr

.github/workflows/run_tests.yml

@@ -27,6 +27,12 @@ jobs:
 #          sudo apt-get update
 #          sudo apt install libasound2-dev libcurl4-openssl-dev libx11-dev libxinerama-dev libxext-dev libfreetype6-dev libwebkit2gtk-4.0-dev libglu1-mesa-dev libsamplerate-dev
 
+      - name: Install Xcode
+        if: runner.os == 'MacOS'
+        uses: maxim-lobanov/setup-xcode@v1
+        with:
+          xcode-version: '14.3.1'
+
       - name: Get latest CMake
         uses: lukka/get-cmake@latest
 

include/math_approx/src/asinh_approx.hpp

@@ -52,7 +52,7 @@ struct AsinhLog2Provider
  * accuracy achieved by the STL implementation).
  */
 template <int order, typename T>
-T asinh (T x)
+constexpr T asinh (T x)
 {
     using S = scalar_of_t<T>;
     using std::abs;

include/math_approx/src/basic_math.hpp

@@ -24,6 +24,8 @@ using scalar_of_t = typename scalar_of<T>::type;
 template <typename T>
 T rsqrt (T x)
 {
+    // @TODO: figure out a way that we can make this method constexpr
+
     // sqrtss followed by divss... this seems to measure a bit faster than the rsqrtss plus NR iteration below
     return (T) 1 / std::sqrt (x);
 

include/math_approx/src/log_approx.hpp

@@ -103,29 +103,29 @@ namespace log_detail
 
 /** approximation for log(Base, x) (32-bit) */
 template <typename Base, int order, bool C1_continuous, typename Log2ProviderType = log_detail::Log2Provider>
-float log (float x)
+constexpr float log (float x)
 {
-    const auto vi = reinterpret_cast<int32_t&> (x);
+    const auto vi = std::bit_cast<int32_t> (x);
     const auto ex = vi & 0x7f800000;
     const auto e = (ex >> 23) - 127;
     const auto vfi = (vi - ex) | 0x3f800000;
-    const auto vf = reinterpret_cast<const float&> (vfi);
+    const auto vf = std::bit_cast<float> (vfi);
 
-    static constexpr auto log2_base_r = 1.0f / Base::log2_base;
+    constexpr auto log2_base_r = 1.0f / Base::log2_base;
     return log2_base_r * ((float) e + Log2ProviderType::template log2_approx<float, order, C1_continuous> (vf));
 }
 
 /** approximation for log(x) (64-bit) */
 template <typename Base, int order, bool C1_continuous, typename Log2ProviderType = log_detail::Log2Provider>
-double log (double x)
+constexpr double log (double x)
 {
-    const auto vi = reinterpret_cast<int64_t&> (x);
+    const auto vi = std::bit_cast<int64_t> (x);
     const auto ex = vi & 0x7ff0000000000000;
     const auto e = (ex >> 52) - 1023;
     const auto vfi = (vi - ex) | 0x3ff0000000000000;
-    const auto vf = reinterpret_cast<const double&> (vfi);
+    const auto vf = std::bit_cast<double> (vfi);
 
-    static constexpr auto log2_base_r = 1.0 / Base::log2_base;
+    constexpr auto log2_base_r = 1.0 / Base::log2_base;
     return log2_base_r * ((double) e + Log2ProviderType::template log2_approx<double, order, C1_continuous> (vf));
 }
 
@@ -134,11 +134,11 @@ double log (double x)
 template <typename Base, int order, bool C1_continuous, typename Log2ProviderType = log_detail::Log2Provider>
 xsimd::batch<float> log (xsimd::batch<float> x)
 {
-    const auto vi = reinterpret_cast<xsimd::batch<int32_t>&> (x); // NOSONAR
+    const auto vi = xsimd::bit_cast<xsimd::batch<int32_t>> (x);
     const auto ex = vi & 0x7f800000;
     const auto e = (ex >> 23) - 127;
     const auto vfi = (vi - ex) | 0x3f800000;
-    const auto vf = reinterpret_cast<const xsimd::batch<float>&> (vfi); // NOSONAR
+    const auto vf = xsimd::bit_cast<xsimd::batch<float>> (vfi);
 
     static constexpr auto log2_base_r = 1.0f / Base::log2_base;
     return log2_base_r * (xsimd::to_float (e) + Log2ProviderType::template log2_approx<xsimd::batch<float>, order, C1_continuous> (vf));
@@ -148,11 +148,11 @@ xsimd::batch<float> log (xsimd::batch<float> x)
 template <typename Base, int order, bool C1_continuous, typename Log2ProviderType = log_detail::Log2Provider>
 xsimd::batch<double> log (xsimd::batch<double> x)
 {
-    const auto vi = reinterpret_cast<xsimd::batch<int64_t>&> (x); // NOSONAR
+    const auto vi = xsimd::bit_cast<xsimd::batch<int64_t>> (x);
     const auto ex = vi & 0x7ff0000000000000;
     const auto e = (ex >> 52) - 1023;
     const auto vfi = (vi - ex) | 0x3ff0000000000000;
-    const auto vf = reinterpret_cast<const xsimd::batch<double>&> (vfi); // NOSONAR
+    const auto vf = xsimd::bit_cast<xsimd::batch<double>> (vfi);
 
     static constexpr auto log2_base_r = 1.0 / Base::log2_base;
     return log2_base_r * (xsimd::to_float (e) + Log2ProviderType::template log2_approx<xsimd::batch<double>, order, C1_continuous> (vf));
@@ -164,19 +164,19 @@ xsimd::batch<double> log (xsimd::batch<double> x)
 #endif
 
 template <int order, bool C1_continuous = false, typename T>
-T log (T x)
+constexpr T log (T x)
 {
     return log<pow_detail::BaseE<scalar_of_t<T>>, order, C1_continuous> (x);
 }
 
 template <int order, bool C1_continuous = false, typename T>
-T log2 (T x)
+constexpr T log2 (T x)
 {
     return log<pow_detail::Base2<scalar_of_t<T>>, order, C1_continuous> (x);
 }
 
 template <int order, bool C1_continuous = false, typename T>
-T log10 (T x)
+constexpr T log10 (T x)
 {
     return log<pow_detail::Base10<scalar_of_t<T>>, order, C1_continuous> (x);
 }

include/math_approx/src/polylogarithm_approx.hpp

@@ -14,7 +14,7 @@ namespace math_approx
  * improve the accuracy very much.
  */
 template <int order, typename T>
-T li2_0_half (T x)
+constexpr T li2_0_half (T x)
 {
     static_assert (order >= 1 && order <= 6);
     using S = scalar_of_t<T>;
@@ -104,13 +104,13 @@ T li2_0_half (T x)
  * improve the accuracy very much.
  */
 template <int order, int log_order = std::min (order + 2, 6), bool log_C1 = (log_order >= 5), typename T>
-T li2 (T x)
+constexpr T li2 (T x)
 {
     const auto x_r = (T) 1 / x;
     const auto x_r1 = (T) 1 / (x - (T) 1);
 
-    static constexpr auto pisq_o_6 = (T) M_PI * (T) M_PI / (T) 6;
-    static constexpr auto pisq_o_3 = (T) M_PI * (T) M_PI / (T) 3;
+    constexpr auto pisq_o_6 = (T) M_PI * (T) M_PI / (T) 6;
+    constexpr auto pisq_o_3 = (T) M_PI * (T) M_PI / (T) 3;
 
     T y, r;
     bool sign = true;

include/math_approx/src/pow_approx.hpp

@@ -139,7 +139,7 @@ namespace pow_detail
 
 /** approximation for pow(Base, x) (32-bit) */
 template <typename Base, int order, bool C1_continuous>
-float pow (float x)
+constexpr float pow (float x)
 {
     x = std::max (-126.0f, Base::log2_base * x);
 
@@ -148,12 +148,12 @@ float pow (float x)
     const auto f = x - (float) l;
     const auto vi = (l + 127) << 23;
 
-    return reinterpret_cast<const float&> (vi) * pow_detail::pow2_approx<float, order, C1_continuous> (f);
+    return std::bit_cast<float> (vi) * pow_detail::pow2_approx<float, order, C1_continuous> (f);
 }
 
 /** approximation for pow(Base, x) (64-bit) */
 template <typename Base, int order, bool C1_continuous>
-double pow (double x)
+constexpr double pow (double x)
 {
     x = std::max (-1022.0, Base::log2_base * x);
 
@@ -162,7 +162,7 @@ double pow (double x)
     const auto d = x - (double) l;
     const auto vi = (l + 1023) << 52;
 
-    return reinterpret_cast<const double&> (vi) * pow_detail::pow2_approx<double, order, C1_continuous> (d);
+    return std::bit_cast<double> (vi) * pow_detail::pow2_approx<double, order, C1_continuous> (d);
 }
 
 #if defined(XSIMD_HPP)
@@ -177,7 +177,7 @@ xsimd::batch<float> pow (xsimd::batch<float> x)
     const auto f = x - xsimd::to_float (l);
     const auto vi = (l + 127) << 23;
 
-    return reinterpret_cast<const xsimd::batch<float>&> (vi) * pow_detail::pow2_approx<xsimd::batch<float>, order, C1_continuous> (f);
+    return xsimd::bit_cast<xsimd::batch<float>> (vi) * pow_detail::pow2_approx<xsimd::batch<float>, order, C1_continuous> (f);
 }
 
 /** approximation for pow(Base, x) (64-bit SIMD) */
@@ -191,7 +191,7 @@ xsimd::batch<double> pow (xsimd::batch<double> x)
     const auto d = x - xsimd::to_float (l);
     const auto vi = (l + 1023) << 52;
 
-    return reinterpret_cast<const xsimd::batch<double>&> (vi) * pow_detail::pow2_approx<xsimd::batch<double>, order, C1_continuous> (d);
+    return xsimd::bit_cast<xsimd::batch<double>> (vi) * pow_detail::pow2_approx<xsimd::batch<double>, order, C1_continuous> (d);
 }
 #endif
 
@@ -200,19 +200,19 @@ xsimd::batch<double> pow (xsimd::batch<double> x)
 #endif
 
 template <int order, bool C1_continuous = false, typename T>
-T exp (T x)
+constexpr T exp (T x)
 {
     return pow<pow_detail::BaseE<scalar_of_t<T>>, order, C1_continuous> (x);
 }
 
 template <int order, bool C1_continuous = false, typename T>
-T exp2 (T x)
+constexpr T exp2 (T x)
 {
     return pow<pow_detail::Base2<scalar_of_t<T>>, order, C1_continuous> (x);
 }
 
 template <int order, bool C1_continuous = false, typename T>
-T exp10 (T x)
+constexpr T exp10 (T x)
 {
     return pow<pow_detail::Base10<scalar_of_t<T>>, order, C1_continuous> (x);
 }

include/math_approx/src/sigmoid_approx.hpp

@@ -9,7 +9,7 @@ namespace sigmoid_detail
     // These polynomial fits were generated from: https://www.wolframcloud.com/obj/chowdsp/Published/sigmoid_approx.nb
 
     template <typename T>
-    T sig_poly_9 (T x)
+    constexpr T sig_poly_9 (T x)
     {
         using S = scalar_of_t<T>;
         const auto x_sq = x * x;
@@ -21,7 +21,7 @@ namespace sigmoid_detail
     }
 
     template <typename T>
-    T sig_poly_7 (T x)
+    constexpr T sig_poly_7 (T x)
     {
         using S = scalar_of_t<T>;
         const auto x_sq = x * x;
@@ -32,7 +32,7 @@ namespace sigmoid_detail
     }
 
     template <typename T>
-    T sig_poly_5 (T x)
+    constexpr T sig_poly_5 (T x)
     {
         using S = scalar_of_t<T>;
         const auto x_sq = x * x;
@@ -42,7 +42,7 @@ namespace sigmoid_detail
     }
 
     template <typename T>
-    T sig_poly_3 (T x)
+    constexpr T sig_poly_3 (T x)
     {
         using S = scalar_of_t<T>;
         const auto x_sq = x * x;

include/math_approx/src/sinh_cosh_approx.hpp

@@ -11,7 +11,7 @@ namespace math_approx
 
 /** Approximation of sinh(x), using exp(x) internally */
 template <int order, typename T>
-T sinh (T x)
+constexpr T sinh (T x)
 {
     using S = scalar_of_t<T>;
     auto B = exp<order> (x);
@@ -22,7 +22,7 @@ T sinh (T x)
 
 /** Approximation of cosh(x), using exp(x) internally */
 template <int order, typename T>
-T cosh (T x)
+constexpr T cosh (T x)
 {
     using S = scalar_of_t<T>;
     auto B = exp<order> (x);
@@ -38,7 +38,7 @@ T cosh (T x)
  * For more information see the comments above.
  */
 template <int order, typename T>
-auto sinh_cosh (T x)
+constexpr auto sinh_cosh (T x)
 {
     using S = scalar_of_t<T>;
     auto B = exp<order> (x);

include/math_approx/src/tanh_approx.hpp

@@ -9,7 +9,7 @@ namespace tanh_detail
     // These polynomial fits were generated from: https://www.wolframcloud.com/obj/chowdsp/Published/tanh_approx.nb
 
     template <typename T>
-    T tanh_poly_11 (T x)
+    constexpr T tanh_poly_11 (T x)
     {
         using S = scalar_of_t<T>;
         const auto x_sq = x * x;
@@ -22,7 +22,7 @@ namespace tanh_detail
     }
 
     template <typename T>
-    T tanh_poly_9 (T x)
+    constexpr T tanh_poly_9 (T x)
     {
         using S = scalar_of_t<T>;
         const auto x_sq = x * x;
@@ -34,7 +34,7 @@ namespace tanh_detail
     }
 
     template <typename T>
-    T tanh_poly_7 (T x)
+    constexpr T tanh_poly_7 (T x)
     {
         using S = scalar_of_t<T>;
         const auto x_sq = x * x;
@@ -45,7 +45,7 @@ namespace tanh_detail
     }
 
     template <typename T>
-    T tanh_poly_5 (T x)
+    constexpr T tanh_poly_5 (T x)
     {
         using S = scalar_of_t<T>;
         const auto x_sq = x * x;
@@ -55,7 +55,7 @@ namespace tanh_detail
     }
 
     template <typename T>
-    T tanh_poly_3 (T x)
+    constexpr T tanh_poly_3 (T x)
     {
         using S = scalar_of_t<T>;
         const auto x_sq = x * x;