From b7f9f747eb9a523d787ada1026ba9b0126a3bd57 Mon Sep 17 00:00:00 2001 From: jatin Date: Fri, 12 Jan 2024 01:05:27 -0800 Subject: [PATCH] Add lots of comments --- include/math_approx/src/basic_math.hpp | 10 +++++++++- .../math_approx/src/hyperbolic_trig_approx.hpp | 6 +++++- .../src/inverse_hyperbolic_trig_approx.hpp | 17 +++++++++++++---- include/math_approx/src/inverse_trig_approx.hpp | 17 +++++++++++++++++ include/math_approx/src/log_approx.hpp | 15 +++++++++++++++ .../math_approx/src/polylogarithm_approx.hpp | 2 ++ include/math_approx/src/pow_approx.hpp | 6 ++++++ include/math_approx/src/sigmoid_approx.hpp | 7 ++++++- include/math_approx/src/trig_approx.hpp | 16 ++++++++++++++-- include/math_approx/src/wright_omega_approx.hpp | 10 ++++++++++ tools/plotter/plotter.cpp | 6 +++--- 11 files changed, 100 insertions(+), 12 deletions(-) diff --git a/include/math_approx/src/basic_math.hpp b/include/math_approx/src/basic_math.hpp index 8a7a501..5f2e711 100644 --- a/include/math_approx/src/basic_math.hpp +++ b/include/math_approx/src/basic_math.hpp @@ -20,9 +20,15 @@ struct scalar_of { using type = T; }; + +/** + * When T is a scalar floating-point type, scalar_of_t is T. + * When T is a SIMD floating-point type, scalar_of_t is the corresponding scalar type. + */ template using scalar_of_t = typename scalar_of::type; +/** Inverse square root */ template T rsqrt (T x) { @@ -40,6 +46,7 @@ T rsqrt (T x) // return x * r; } +/** Function interface for the ternary operator. */ template T select (bool q, T t, T f) { @@ -53,6 +60,7 @@ struct scalar_of> using type = T; }; +/** Inverse square root */ template xsimd::batch rsqrt (xsimd::batch x) { @@ -65,6 +73,7 @@ xsimd::batch rsqrt (xsimd::batch x) return x * r; } +/** Function interface for the ternary operator. */ template xsimd::batch select (xsimd::batch_bool q, xsimd::batch t, xsimd::batch f) { @@ -94,5 +103,4 @@ inline typename std::enable_if::value, To>::type b #else using std::bit_cast; #endif - } // namespace math_approx diff --git a/include/math_approx/src/hyperbolic_trig_approx.hpp b/include/math_approx/src/hyperbolic_trig_approx.hpp index 059b921..0f1960d 100644 --- a/include/math_approx/src/hyperbolic_trig_approx.hpp +++ b/include/math_approx/src/hyperbolic_trig_approx.hpp @@ -53,7 +53,7 @@ constexpr auto sinh_cosh (T x) namespace tanh_detail { - // These polynomial fits were generated from: https://www.wolframcloud.com/obj/chowdsp/Published/tanh_approx.nb + // See notebooks/tanh_approx.nb for the derivation of these polynomials template constexpr T tanh_poly_11 (T x) @@ -111,6 +111,10 @@ namespace tanh_detail } } // namespace tanh_detail +/** + * Approximation of tanh(x), using tanh(x) ≈ p(x) / (p(x)^2 + 1), + * where p(x) is an odd polynomial fit to minimize the maxinimum relative error. + */ template T tanh (T x) { diff --git a/include/math_approx/src/inverse_hyperbolic_trig_approx.hpp b/include/math_approx/src/inverse_hyperbolic_trig_approx.hpp index ea486e5..78df2ed 100644 --- a/include/math_approx/src/inverse_hyperbolic_trig_approx.hpp +++ b/include/math_approx/src/inverse_hyperbolic_trig_approx.hpp @@ -7,6 +7,8 @@ namespace math_approx { struct AsinhLog2Provider { + // for polynomial derivations, see notebooks/asinh_approx.nb + /** approximation for log2(x), optimized on the range [1, 2], to be used within an asinh(x) computation */ template static constexpr T log2_approx (T x) @@ -76,6 +78,10 @@ constexpr T asinh (T x) return sign * y; } +/** + * Approximation of acosh(x) in the full range, using identity + * acosh(x) = log(x + sqrt(x^2 - 1)). + */ template constexpr T acosh (T x) { @@ -85,15 +91,18 @@ constexpr T acosh (T x) using xsimd::sqrt; #endif - const auto z0 = x - (S) 1; - const auto z1 = z0 + sqrt (z0 + z0 + z0 * z0); - return log1p (z1); + const auto z1 = x + sqrt (x * x - (S) 1); + return log (z1); } +/** + * Approximation of atanh(x), using identity + * atanh(x) = (1/2) log((x + 1) / (x - 1)). + */ template constexpr T atanh (T x) { using S = scalar_of_t; return (S) 0.5 * log (((S) 1 + x) / ((S) 1 - x)); } -} +} // namespace math_approx diff --git a/include/math_approx/src/inverse_trig_approx.hpp b/include/math_approx/src/inverse_trig_approx.hpp index 432c386..1fe65c7 100644 --- a/include/math_approx/src/inverse_trig_approx.hpp +++ b/include/math_approx/src/inverse_trig_approx.hpp @@ -6,6 +6,8 @@ namespace math_approx { namespace inv_trig_detail { + // for polynomial derivations, see notebooks/asin_acos_approx.nb + template constexpr T asin_kernel (T x) { @@ -66,6 +68,8 @@ namespace inv_trig_detail } } + // for polynomial derivations, see notebooks/arctan_approx.nb + template constexpr T atan_kernel (T x) { @@ -100,6 +104,11 @@ namespace inv_trig_detail } } // namespace inv_trig_detail +/** + * Approximation of asin(x) using asin(x) ≈ p(x^2) * x^3 + x for x in [0, 0.5], + * and asin(x) ≈ pi/2 - p((1-x)/2) * ((1-x)/2)^3/2 + ((1-x)/2)^1/2 for x in [0.5, 1], + * where p(x) is a polynomial fit to achieve the minimum absolute error. + */ template T asin (T x) { @@ -123,6 +132,10 @@ T asin (T x) return select (x > (S) 0, res, -res); } +/** + * Approximation of acos(x) using the same approach as asin(x), + * but with a different polynomial fit. + */ template T acos (T x) { @@ -146,6 +159,10 @@ T acos (T x) return (S) M_PI_2 - select (x > (S) 0, res, -res); } +/** + * Approximation of atan(x) using a polynomial approximation of arctan(x) on [0, 1], + * and atan(x) = pi/2 - arctan(1/x) for x > 1. + */ template T atan (T x) { diff --git a/include/math_approx/src/log_approx.hpp b/include/math_approx/src/log_approx.hpp index 9ba64d6..d58ec4e 100644 --- a/include/math_approx/src/log_approx.hpp +++ b/include/math_approx/src/log_approx.hpp @@ -9,6 +9,8 @@ namespace log_detail { struct Log2Provider { + // for polynomial derivations, see notebooks/log_approx.nb + /** approximation for log2(x), optimized on the range [1, 2] */ template static constexpr T log2_approx (T x) @@ -163,24 +165,37 @@ xsimd::batch log (xsimd::batch x) #pragma GCC diagnostic pop // end ignore strict-aliasing warnings #endif +/** + * Approximation of log(x), using + * log(x) = (1 / log2(e)) * (Exponent(x) + log2(1 + Mantissa(x)) + */ template constexpr T log (T x) { return log>, order, C1_continuous> (x); } +/** + * Approximation of log2(x), using + * log2(x) = Exponent(x) + log2(1 + Mantissa(x) + */ template constexpr T log2 (T x) { return log>, order, C1_continuous> (x); } +/** + * Approximation of log10(x), using + * log10(x) = (1 / log2(10)) * (Exponent(x) + log2(1 + Mantissa(x)) + */ template constexpr T log10 (T x) { return log>, order, C1_continuous> (x); } +/** Approximation of log(1 + x), using math_approx::log(x) */ template constexpr T log1p (T x) { diff --git a/include/math_approx/src/polylogarithm_approx.hpp b/include/math_approx/src/polylogarithm_approx.hpp index 2339436..a58d796 100644 --- a/include/math_approx/src/polylogarithm_approx.hpp +++ b/include/math_approx/src/polylogarithm_approx.hpp @@ -12,6 +12,8 @@ namespace math_approx * Orders higher than 3 are generally not recommended for * single-precision floating-point types, since they don't * improve the accuracy very much. + * + * For derivations, see notebooks/li2_approx.nb */ template constexpr T li2_0_half (T x) diff --git a/include/math_approx/src/pow_approx.hpp b/include/math_approx/src/pow_approx.hpp index 68c218d..51c5a3a 100644 --- a/include/math_approx/src/pow_approx.hpp +++ b/include/math_approx/src/pow_approx.hpp @@ -6,6 +6,8 @@ namespace math_approx { namespace pow_detail { + // for polynomial derivations, see notebooks/exp_approx.nb + /** approximation for 2^x, optimized on the range [0, 1] */ template constexpr T pow2_approx (T x) @@ -199,24 +201,28 @@ xsimd::batch pow (xsimd::batch x) #pragma GCC diagnostic pop // end ignore strict-aliasing warnings #endif +/** Approximation of exp(x), using exp(x) = 2^floor(x * log2(e)) * 2^frac(x * log2(e)) */ template constexpr T exp (T x) { return pow>, order, C1_continuous> (x); } +/** Approximation of exp2(x), using exp(x) = 2^floor(x) * 2^frac(x) */ template constexpr T exp2 (T x) { return pow>, order, C1_continuous> (x); } +/** Approximation of exp(x), using exp10(x) = 2^floor(x * log2(10)) * 2^frac(x * log2(10)) */ template constexpr T exp10 (T x) { return pow>, order, C1_continuous> (x); } +/** Approximation of exp(1) - 1, using math_approx::exp(x) */ template constexpr T expm1 (T x) { diff --git a/include/math_approx/src/sigmoid_approx.hpp b/include/math_approx/src/sigmoid_approx.hpp index ac5252d..7605f57 100644 --- a/include/math_approx/src/sigmoid_approx.hpp +++ b/include/math_approx/src/sigmoid_approx.hpp @@ -6,7 +6,7 @@ namespace math_approx { namespace sigmoid_detail { - // These polynomial fits were generated from: https://www.wolframcloud.com/obj/chowdsp/Published/sigmoid_approx.nb + // for polynomial derivations, see notebooks/sigmoid_approx.nb template constexpr T sig_poly_9 (T x) @@ -51,6 +51,11 @@ namespace sigmoid_detail } } // namespace sigmoid_detail +/** + * Approximation of sigmoid(x) := 1 / (1 + e^-x), + * using sigmoid(x) ≈ (1/2) p(x) / (p(x)^2 + 1) + (1/2), + * where p(x) is an odd polynomial fit to minimize the maxinimum relative error. + */ template T sigmoid (T x) { diff --git a/include/math_approx/src/trig_approx.hpp b/include/math_approx/src/trig_approx.hpp index d74db9e..e93257b 100644 --- a/include/math_approx/src/trig_approx.hpp +++ b/include/math_approx/src/trig_approx.hpp @@ -48,6 +48,8 @@ namespace trig_detail return select (x >= (T) 0, mod, mod + pi) - half_pi; } + // for polynomial derivations, see notebooks/sin_approx.nb + template constexpr T sin_poly_9 (T x, T x_sq) { @@ -79,6 +81,7 @@ namespace trig_detail } } // namespace sin_detail +/** Polynomial approximation of sin(x) on the range [-pi, pi] */ template constexpr T sin_mpi_pi (T x) { @@ -100,12 +103,17 @@ constexpr T sin_mpi_pi (T x) return (pi_sq - x_sq) * x_poly; } +/** Full range approximation of sin(x) */ template constexpr T sin (T x) { return sin_mpi_pi (trig_detail::fast_mod_mpi_pi (x)); } +/** + * Polynomial approximation of cos(x) on the range [-pi, pi], + * using a range-shifted approximation of sin(x). + */ template constexpr T cos_mpi_pi (T x) { @@ -136,18 +144,21 @@ constexpr T cos_mpi_pi (T x) return (pi_sq - hpmx_sq) * x_poly; } +/** Full range approximation of cos(x) */ template constexpr T cos (T x) { return cos_mpi_pi (trig_detail::fast_mod_mpi_pi (x)); } -/** Approximation of tan(x) on the range [-pi/4, pi/4] */ +/** Polynomial approximation of tan(x) on the range [-pi/4, pi/4] */ template constexpr T tan_mquarterpi_quarterpi (T x) { static_assert (order % 2 == 1 && order >= 3 && order <= 15, "Order must be an odd number within [3, 15]"); + // for polynomial derivation, see notebooks/tan_approx.nb + using S = scalar_of_t; const auto x_sq = x * x; if constexpr (order == 3) @@ -217,7 +228,8 @@ constexpr T tan_mquarterpi_quarterpi (T x) } /** - * Approximation of tan(x) on the range [-pi/2, pi/2] + * Approximation of tan(x) on the range [-pi/2, pi/2], + * using the tangent half-angle formula. * * Accuracy may suffer as x approaches ±pi/2. */ diff --git a/include/math_approx/src/wright_omega_approx.hpp b/include/math_approx/src/wright_omega_approx.hpp index b14c4f0..eecb199 100644 --- a/include/math_approx/src/wright_omega_approx.hpp +++ b/include/math_approx/src/wright_omega_approx.hpp @@ -4,6 +4,16 @@ namespace math_approx { +/** + * Approximation of the Wright-Omega function, using + * w(x) ≈ 0 for x < -3 + * w(x) ≈ p(x) for -3 <= x < e + * w(x) ≈ x - log(x) + alpha * exp(-beta * x) for x >= e, + * where p(x) is a polynomial, and alpha and beta are coefficients, + * all fit to minimize the maximum absolute error. + * + * The above fit is optionally followed by some number of Newton-Raphson iterations. + */ template constexpr T wright_omega (T x) { diff --git a/tools/plotter/plotter.cpp b/tools/plotter/plotter.cpp index cca6d0c..be727cb 100644 --- a/tools/plotter/plotter.cpp +++ b/tools/plotter/plotter.cpp @@ -61,12 +61,12 @@ void plot_function (std::span all_floats, int main() { plt::figure(); - const auto range = std::make_pair (-0.99f, 0.99f); + const auto range = std::make_pair (1.0f, 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::atanh)); - plot_error (all_floats, y_exact, FLOAT_FUNC ((math_approx::atanh<5>) ), "atanh_log-5"); + const auto y_exact = test_helpers::compute_all (all_floats, FLOAT_FUNC (std::acosh)); + plot_error (all_floats, y_exact, FLOAT_FUNC ((math_approx::acosh<5>) ), "acosh-5"); plt::legend ({ { "loc", "upper right" } }); plt::xlim (range.first, range.second);