-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
5653c28
commit 449a1c5
Showing
8 changed files
with
286 additions
and
89 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,79 @@ | ||
| #pragma once | ||
|
|
||
| #include "log_approx.hpp" | ||
|
|
||
| namespace math_approx | ||
| { | ||
| struct AsinhLog2Provider | ||
| { | ||
| /** approximation for log2(x), optimized on the range [1, 2], to be used within an asinh(x) computation */ | ||
| template <typename T, int order, bool /*C1_continuous*/> | ||
| static constexpr T log2_approx (T x) | ||
| { | ||
| static_assert (order >= 3 && order <= 5); | ||
| using S = scalar_of_t<T>; | ||
|
|
||
| const auto x_sq = x * x; | ||
| if constexpr (order == 3) | ||
| { | ||
| const auto x_2_3 = (S) -1.21535595794871 + (S) 0.194363894384581 * x; | ||
| const auto x_0_1 = (S) -2.26452854958994 + (S) 3.28552061315407 * x; | ||
| return x_0_1 + x_2_3 * x_sq; | ||
| } | ||
| else if constexpr (order == 4) | ||
| { | ||
| const auto x_3_4 = (S) 0.770443387059628 + (S) -0.102652345633016 * x; | ||
| const auto x_1_2 = (S) 4.33013912645867 + (S) -2.39448588379361 * x; | ||
| const auto x_1_2_3_4 = x_1_2 + x_3_4 * x_sq; | ||
| return (S) -2.60344428409168 + x_1_2_3_4 * x; | ||
| } | ||
| else if constexpr (order == 5) | ||
| { | ||
| const auto x_4_5 = (S) -0.511946284688366 + (S) 0.0578217518982235 * x; | ||
| const auto x_2_3 = (S) -3.94632584968643 + (S) 1.90796087279737 * x; | ||
| const auto x_0_1 = (S) -2.87748189127908 + (S) 5.36997140095829 * 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 | ||
| { | ||
| return {}; | ||
| } | ||
| } | ||
| }; | ||
|
|
||
| /** | ||
| * Approximation of asinh(x) in the full range, using identity | ||
| * asinh(x) = log(x + sqrt(x^2 + 1)). | ||
| * | ||
| * Orders 6 and 7 use an additional Newton-Raphson iteration, | ||
| * but for most cases the accuracy improvement is not worth | ||
| * the additional cost (when compared to the performance and | ||
| * accuracy achieved by the STL implementation). | ||
| */ | ||
| template <int order, typename T> | ||
| T asinh (T x) | ||
| { | ||
| using S = scalar_of_t<T>; | ||
| using std::abs; | ||
| using std::sqrt; | ||
| #if defined(XSIMD_HPP) | ||
| using xsimd::abs; | ||
| using xsimd::sqrt; | ||
| #endif | ||
|
|
||
| const auto sign = select (x > (S) 0, (T) (S) 1, select (x < (S) 0, (T) (S) -1, (T) (S) 0)); | ||
| x = abs (x); | ||
|
|
||
| const auto log_arg = x + sqrt (x * x + (S) 1); | ||
| auto y = log<pow_detail::BaseE<scalar_of_t<T>>, std::min (order, 5), false, AsinhLog2Provider> (log_arg); | ||
|
|
||
| if constexpr (order > 5) | ||
| { | ||
| const auto exp_y = math_approx::exp<order - 1> (y); | ||
| y -= (exp_y - log_arg) / exp_y; | ||
| } | ||
|
|
||
| return sign * y; | ||
| } | ||
| } // namespace math_approx |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,54 @@ | ||
| #include "test_helpers.hpp" | ||
| #include "catch2/catch_template_test_macros.hpp" | ||
|
|
||
| #include <catch2/catch_test_macros.hpp> | ||
| #include <iostream> | ||
|
|
||
| #include <math_approx/math_approx.hpp> | ||
|
|
||
| template <typename T = float> | ||
| void test_approx (const auto& all_floats, const auto& y_exact, auto&& f_approx, float err_bound) | ||
| { | ||
| const auto y_approx = test_helpers::compute_all<T> (all_floats, f_approx); | ||
| const auto error = test_helpers::compute_error<T> (y_exact, y_approx); | ||
| const auto max_error = test_helpers::abs_max<T> (error); | ||
|
|
||
| std::cout << max_error << std::endl; | ||
| REQUIRE (std::abs (max_error) < err_bound); | ||
| } | ||
|
|
||
| TEMPLATE_TEST_CASE ("Asinh Approx Test", "", float, double) | ||
| { | ||
| #if ! defined(WIN32) | ||
| const auto all_floats = test_helpers::all_32_bit_floats<TestType> (-10.0f, 10.0f, 1.0e-2f); | ||
| #else | ||
| const auto all_floats = test_helpers::all_32_bit_floats (-10.0f, 10.0f, 1.0e-1f); | ||
| #endif | ||
| const auto y_exact = test_helpers::compute_all<TestType> (all_floats, [] (auto x) | ||
| { return std::asinh (x); }); | ||
|
|
||
| SECTION ("6th-Order") | ||
| { | ||
| test_approx<TestType> (all_floats, y_exact, [] (auto x) | ||
| { return math_approx::asinh<6> (x); }, | ||
| 5.0e-7f); | ||
| } | ||
| SECTION ("5th-Order") | ||
| { | ||
| test_approx<TestType> (all_floats, y_exact, [] (auto x) | ||
| { return math_approx::asinh<5> (x); }, | ||
| 6.0e-5f); | ||
| } | ||
| SECTION ("4th-Order") | ||
| { | ||
| test_approx<TestType> (all_floats, y_exact, [] (auto x) | ||
| { return math_approx::asinh<4> (x); }, | ||
| 3.5e-4f); | ||
| } | ||
| SECTION ("3th-Order") | ||
| { | ||
| test_approx<TestType> (all_floats, y_exact, [] (auto x) | ||
| { return math_approx::asinh<3> (x); }, | ||
| 2.5e-3f); | ||
| } | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.