From 7796ba7e799aa1d2943e378177140fd057366db0 Mon Sep 17 00:00:00 2001
From: jatin <jatinchowdhury18@gmail.com>
Date: Wed, 22 Nov 2023 18:50:15 -0800
Subject: [PATCH] First pass at log approximations

---
 include/math_approx/math_approx.hpp    |  1 +
 include/math_approx/src/log_approx.hpp | 88 ++++++++++++++++++++++++++
 include/math_approx/src/pow_approx.hpp |  2 -
 tools/plotter/plotter.cpp              | 11 ++--
 4 files changed, 95 insertions(+), 7 deletions(-)
 create mode 100644 include/math_approx/src/log_approx.hpp

diff --git a/include/math_approx/math_approx.hpp b/include/math_approx/math_approx.hpp
index 3652539..9a3f06f 100644
--- a/include/math_approx/math_approx.hpp
+++ b/include/math_approx/math_approx.hpp
@@ -10,3 +10,4 @@ namespace math_approx
 #include "src/sigmoid_approx.hpp"
 #include "src/trig_approx.hpp"
 #include "src/pow_approx.hpp"
+#include "src/log_approx.hpp"
diff --git a/include/math_approx/src/log_approx.hpp b/include/math_approx/src/log_approx.hpp
new file mode 100644
index 0000000..4861a12
--- /dev/null
+++ b/include/math_approx/src/log_approx.hpp
@@ -0,0 +1,88 @@
+#pragma once
+
+#include "basic_math.hpp"
+#include "pow_approx.hpp"
+
+namespace math_approx
+{
+namespace log_detail
+{
+    /** approximation for log2(x), optimized on the range [1, 2] */
+    template <typename T, int order>
+    constexpr T log2_approx (T x)
+    {
+        static_assert (order >= 3 && order <= 6);
+        using S = scalar_of_t<T>;
+
+        const auto x_sq = x * x;
+        if constexpr (order == 3)
+        {
+            const auto x_2_3 = (S) -1.05974531422 + (S) 0.159220010975 * x;
+            const auto x_0_1 = (S) -2.16417056258 + (S) 3.06469586582 * x;
+            return x_0_1 + x_2_3 * x_sq;
+        }
+        else if constexpr (order == 4)
+        {
+            const auto x_3_4 = (S) 0.649709537672 + (S) -0.0821303550902 * x;
+            const auto x_1_2 = (S) 4.08637809379 + (S) -2.13412984371 * x;
+            const auto x_1_2_3_4 = x_1_2 + x_3_4 * x_sq;
+            return (S) -2.51982743265 + x_1_2_3_4 * x;
+        }
+        else if constexpr (order == 5)
+        {
+            const auto x_4_5 = (S) -0.419319345483 + (S) 0.0451488402558 * x;
+            const auto x_2_3 = (S) -3.56885211615 + (S) 1.64139451414 * x;
+            const auto x_0_1 = (S) -2.80534277658 + (S) 5.10697088382 * x;
+            const auto x_2_3_4_5 = x_2_3 + x_4_5 * x_sq;
+            return x_0_1 + x_2_3_4_5 * x_sq;
+        }
+        else if constexpr (order == 6)
+        {
+            const auto x_5_6 = (S) 0.276834061071 + (S) -0.0258400886535 * x;
+            const auto x_3_4 = (S) 3.30388341157 + (S) -1.27446900713 * x;
+            const auto x_1_2 = (S) 6.12708086513 + (S) -5.36371998242 * x;
+            const auto x_3_4_5_6 = x_3_4 + x_5_6 * x_sq;
+            const auto x_1_2_3_4_5_6 = x_1_2 + x_3_4_5_6 * x_sq;
+            return (S) -3.04376925958 + x_1_2_3_4_5_6 * x;
+        }
+        else
+        {
+            return {};
+        }
+    }
+}
+
+/** approximation for log(Base, x) (32-bit) */
+template <typename Base, int order>
+float log (float x)
+{
+    const auto vi = reinterpret_cast<int32_t&> (x); // NOSONAR
+    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); // NOSONAR
+
+    static constexpr auto log2_base_r = 1.0f / Base::log2_base;
+    return log2_base_r * ((float) e + log_detail::log2_approx<float, order> (vf));
+}
+
+/** approximation for log(x) (64-bit) */
+template <typename Base, int order>
+double log (double x)
+{
+    const auto vi = reinterpret_cast<int64_t&> (x); // NOSONAR
+    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); // NOSONAR
+
+    static constexpr auto log2_base_r = 1.0 / Base::log2_base;
+    return log2_base_r * ((double) e + log_detail::log2_approx<double, order> (vf));
+}
+
+template <int order, typename T>
+T log (T x)
+{
+    return log<pow_detail::BaseE<scalar_of_t<T>>, order> (x);
+}
+}
diff --git a/include/math_approx/src/pow_approx.hpp b/include/math_approx/src/pow_approx.hpp
index 718049e..dde3418 100644
--- a/include/math_approx/src/pow_approx.hpp
+++ b/include/math_approx/src/pow_approx.hpp
@@ -2,8 +2,6 @@
 
 #include "basic_math.hpp"
 
-// Everything in this file is still WIP!
-
 namespace math_approx
 {
 namespace pow_detail
diff --git a/tools/plotter/plotter.cpp b/tools/plotter/plotter.cpp
index e4ca51f..879533b 100644
--- a/tools/plotter/plotter.cpp
+++ b/tools/plotter/plotter.cpp
@@ -59,15 +59,16 @@ void plot_function (std::span<const float> all_floats,
 int main()
 {
     plt::figure();
-    const auto range = std::make_pair (-10.0f, 10.0f);
+    const auto range = std::make_pair (0.01f, 10.0f);
     static constexpr auto tol = 1.0e-2f;
 
     const auto all_floats = test_helpers::all_32_bit_floats (range.first, range.second, tol);
-    const auto y_exact = test_helpers::compute_all (all_floats, FLOAT_FUNC(std::exp2));
+    const auto y_exact = test_helpers::compute_all (all_floats, FLOAT_FUNC(std::log));
 
-    plot_ulp_error (all_floats, y_exact, FLOAT_FUNC(math_approx::exp2<4>), "Exp2-4");
-    plot_ulp_error (all_floats, y_exact, FLOAT_FUNC(math_approx::exp2<5>), "Exp2-5");
-    plot_ulp_error (all_floats, y_exact, FLOAT_FUNC(math_approx::exp2<6>), "Exp2-6");
+    // plot_ulp_error (all_floats, y_exact, FLOAT_FUNC(math_approx::exp2<4>), "Exp2-4");
+    // plot_ulp_error (all_floats, y_exact, FLOAT_FUNC(math_approx::exp2<5>), "Exp2-5");
+    // plot_ulp_error (all_floats, y_exact, FLOAT_FUNC(math_approx::exp2<6>), "Exp2-6");
+    plot_error (all_floats, y_exact, FLOAT_FUNC (math_approx::log<6>), "Log-6");
 
     plt::legend ({ { "loc", "upper right" } });
     plt::xlim (range.first, range.second);