diff --git a/.gitignore b/.gitignore index 98b57a2..48854f7 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ .vscode/ build*/ + +.DS_Store diff --git a/README.md b/README.md index bcd86e7..228783a 100644 --- a/README.md +++ b/README.md @@ -2,4 +2,11 @@ A C++ library for math approximations. -More info coming soon! +Currently supported: + +- sin/cos +- exp/exp2/exp10 +- log/log2/log10 +- tanh +- sigmoid +- Wright-Omega function diff --git a/include/math_approx/math_approx.hpp b/include/math_approx/math_approx.hpp index 9a3f06f..c68ad8f 100644 --- a/include/math_approx/math_approx.hpp +++ b/include/math_approx/math_approx.hpp @@ -11,3 +11,4 @@ namespace math_approx #include "src/trig_approx.hpp" #include "src/pow_approx.hpp" #include "src/log_approx.hpp" +#include "src/wright_omega_approx.hpp" diff --git a/include/math_approx/src/wright_omega_approx.hpp b/include/math_approx/src/wright_omega_approx.hpp new file mode 100644 index 0000000..287e17e --- /dev/null +++ b/include/math_approx/src/wright_omega_approx.hpp @@ -0,0 +1,84 @@ +#pragma once + +#include "basic_math.hpp" + +namespace math_approx +{ +template +T wright_omega (T x) +{ + static_assert (poly_order == 3 || poly_order == 5); + + using S = scalar_of_t; + static constexpr auto E = (S) 2.7182818284590452354; + + const auto x1 = [] (T _x) + { + const auto x_sq = _x * _x; + if constexpr (poly_order == 3) + { + const auto y_2_3 = (S) 0.0534379648805832 + (S) -0.00251076420630778 * _x; + const auto y_0_1 = (S) 0.616522951065868 + (S) 0.388418422853809 * _x; + return y_0_1 + y_2_3 * x_sq; + } + else if constexpr (poly_order == 5) + { + const auto y_4_5 = (S) -0.00156418794118294 + (S) -0.00151562297325209 * _x; + const auto y_2_3 = (S) 0.0719291313363515 + (S) 0.0216881206167543 * _x; + const auto y_0_1 = (S) 0.569291529016010 + (S) 0.290890537885083 * _x; + const auto y_2_3_4_5 = y_2_3 + y_4_5 * x_sq; + return y_0_1 + y_2_3_4_5 * x_sq; + } + else + { + return T {}; + } + }(x); + const auto x2 = x - log (x) + (S) 0.32352057096397160124 * exp ((S) -0.029614177658043381316 * x); + + auto y = select (x < (S) -3, T {}, select (x < (S) E, x1, x2)); + + const auto nr_update = [] (T _x, T _y) + { + return _y - (_y - exp (_x - _y)) / (_y + (S) 1); + }; + + for (int i = 0; i < num_nr_iters; ++i) + y = nr_update (x, y); + + return y; +} + +/** + * Wright-Omega function using Stephano D'Angelo's derivation (https://www.dafx.de/paper-archive/2019/DAFx2019_paper_5.pdf) + * With `num_nr_iters == 0`, this is the fastest implementation, but the least accurate. + * With `num_nr_iters == 1`, this is faster than the other implementation with 0 iterations, and little bit more accurate. + * For more accuracy, use the other implementation with at least 1 NR iteration. + */ +template +T wright_omega_dangelo (T x) +{ + using S = scalar_of_t; + + const auto x1 = [] (T _x) + { + const auto x_sq = _x * _x; + const auto y_2_3 = (S) 4.775931364975583e-2 + (S) -1.314293149877800e-3 * _x; + const auto y_0_1 = (S) 6.313183464296682e-1 + (S) 3.631952663804445e-1 * _x; + return y_0_1 + y_2_3 * x_sq; + }(x); + const auto x2 = x - log (x); + + auto y = select (x < (S) -3.341459552768620, T {}, select (x < (S) 8, x1, x2)); + + const auto nr_update = [] (T _x, T _y) + { + return _y - (_y - exp (_x - _y)) / (_y + (S) 1); + }; + + for (int i = 0; i < num_nr_iters; ++i) + y = nr_update (x, y); + + return y; +} +} // namespace math_approx diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 622920c..8dc0101 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -33,3 +33,4 @@ setup_catch_test(sigmoid_approx_test) setup_catch_test(trig_approx_test) setup_catch_test(pow_approx_test) setup_catch_test(log_approx_test) +setup_catch_test(wright_omega_approx_test) diff --git a/test/src/reference/toms917.hpp b/test/src/reference/toms917.hpp new file mode 100644 index 0000000..69be3e0 --- /dev/null +++ b/test/src/reference/toms917.hpp @@ -0,0 +1,383 @@ +# include +# include +# include +# include +# include +# include +# include +# include + +namespace toms917 +{ +using namespace std; +// +// DBL_EPSILON, provided by , is twice the machine epsilon for +// double precision arithmetic. +// +# define TWOITERTOL DBL_EPSILON + +//****************************************************************************80 + +inline int wrightomega_ext ( complex z, complex &w, + complex &e, complex &r, complex &cond ) + +//****************************************************************************80 +// +// Purpose: +// +// WRIGHTOMEGA_EXT computes the Wright Omega function with extra information. +// +// Discussion: +// +// WRIGHTOMEGA_EXT is the extended routine for evaluating the Wright +// Omega function with the option of extracting the last update step, +// the penultimate residual and the condition number estimate. +// +// Modified: +// +// 14 May 2016 +// +// Author: +// +// Piers Lawrence, Robert Corless, David Jeffrey +// +// Reference: +// +// Piers Lawrence, Robert Corless, David Jeffrey, +// Algorithm 917: Complex Double-Precision Evaluation of the Wright Omega +// Function, +// ACM Transactions on Mathematical Software, +// Volume 38, Number 3, Article 20, April 2012, 17 pages. +// +// Parameters: +// +// Input, complex Z, value at which to evaluate Wrightomega(). +// +// Output, complex &W, the value of Wrightomega(z). +// +// Output, complex &E, the last update step in the iterative scheme. +// +// Output, complex &R, the penultimate residual, +// r_k = z - w_k - log(w_k) +// +// Output, complex &COND, the condition number estimate. +// +// Output, int WRIGHTOMEGA_EXT, error flag; +// 0, successful computation. +// nonzero, the computation failed. +// +{ + double near; + double pi = M_PI; + complex pz; + double s = 1.0; + complex t; + complex wp1; + double x; + double y; + double ympi; + double yppi; +// +// Extract real and imaginary parts of Z. +// + x = real ( z ); + y = imag ( z ); +// +// Compute if we are near the branch cuts. +// + ympi = y - pi; + yppi = y + pi; + near = 0.01; +// +// Test for floating point exceptions: +// + +// +// NaN output for NaN input. +// + if ( isnan ( x ) || isnan ( y ) ) + { + // w = complex ( ( 0.0 / 0.0 ), ( 0.0 / 0.0 ) ); + // e = complex ( 0.0, 0.0 ); + // r = complex ( 0.0, 0.0 ); + return 0; + } +// +// Signed zeros between branches. +// + else if ( isinf ( x ) && ( x < 0.0 ) && ( - pi < y ) && ( y <= pi ) ) + { + if ( fabs ( y ) <= pi / 2.0 ) + { + w = + 0.0; + } + else + { + w = - 0.0; + } + + if ( 0.0 <= y ) + { + w = w + complex ( 0.0, 0.0 ); + } + else + { + w = w + complex ( 0.0, - 1.0 * 0.0 ); + } + + e = complex ( 0.0, 0.0 ); + r = complex ( 0.0, 0.0 ); + return 0; + } +// +// Asymptotic for large z. +// + else if ( isinf ( x ) || isinf ( y ) ) + { + w = complex ( x, y ); + e = complex ( 0.0, 0.0 ); + r = complex ( 0.0, 0.0 ); + return 0; + } +// +// Test if exactly on the singular points. +// + if ( ( x == - 1.0 ) && ( fabs ( y ) == pi ) ) + { + w = complex ( - 1.0, 0.0 ); + e = complex ( 0.0, 0.0 ); + r = complex ( 0.0, 0.0 ); + return 0; + } +// +// Choose approximation based on region. +// + +// +// Region 1: upper branch point. +// Series about z=-1+Pi*I. +// + if ( ( - 2.0 < x && x <= 1.0 && 1.0 < y && y < 2.0 * pi ) ) + { + pz = conj ( sqrt ( conj ( 2.0 * ( z + complex ( 1.0, - pi ) ) ) ) ); + + w = - 1.0 + + ( complex ( 0.0, 1.0 ) + + ( 1.0 / 3.0 + + ( - 1.0 / 36.0 * complex ( 0.0, 1.0 ) + + ( 1.0 / 270.0 + 1.0 / 4320.0 * complex ( 0.0, 1.0 ) * pz ) + * pz ) * pz ) * pz ) * pz; + } +// +// Region 2: lower branch point. +// Series about z=-1-Pi*I. +// + else if ( ( - 2.0 < x && x <= 1.0 && - 2.0 * pi < y && y <- 1.0 ) ) + { + pz = conj ( sqrt ( conj ( 2.0 * ( z + 1.0 + complex ( 0.0, pi ) ) ) ) ); + + w = - 1.0 + + ( - complex ( 0.0, 1.0 ) + ( 1.0 / 3.0 + + ( 1.0 / 36.0 * complex ( 0.0, 1.0 ) + + ( 1.0 / 270.0 - 1.0 / 4320.0 * complex ( 0.0, 1.0 ) * pz ) + * pz ) * pz ) * pz ) * pz; + } +// +// Region 3: between branch cuts. +// Series: About -infinity. +// + else if ( x <= - 2.0 && - pi < y && y <= pi ) + { + pz = exp ( z ); + w = ( 1.0 + + ( - 1.0 + + ( 3.0 / 2.0 + + ( - 8.0 / 3.0 + + 125.0 / 24.0 * pz ) * pz ) * pz ) * pz ) * pz; + } +// +// Region 4: Mushroom. +// Series about z=1. +// + else if ( ( ( - 2.0 < x ) && ( x <= 1.0 ) && ( - 1.0 <= y ) && ( y <= 1.0 ) ) + || ( ( - 2.0 < x ) && ( x - 1.0 ) * ( x - 1.0 ) + y * y <= pi * pi ) ) + { + pz = z - 1.0; + w = 1.0 / 2.0 + 1.0 / 2.0 * z + + ( 1.0 / 16.0 + + ( - 1.0 / 192.0 + + ( - 1.0 / 3072.0 + 13.0 / 61440.0 * pz ) * pz ) * pz ) * pz * pz; + } +// +// Region 5: Top wing. +// Negative log series. +// + else if ( x <= - 1.05 && pi < y && y - pi <= - 0.75 * ( x + 1.0 ) ) + { + t = z - complex ( 0.0, pi ); + pz = log ( - t ); + w = ( ( 1.0 + ( - 3.0 / 2.0 + 1.0 / 3.0 * pz ) * pz ) * pz + + ( ( -1.0 + 1.0 / 2.0 * pz ) * pz + ( pz + ( - pz + t ) * t ) * t ) * t ) + / ( t * t * t ); + } +// +// Region 6: Bottom wing. +// Negative log series. +// + else if ( x <= - 1.05 && 0.75 * ( x + 1.0 ) < y + pi && y + pi <= 0.0 ) + { + t = z + complex ( 0.0, pi ); + pz = log ( - t ); + w = ( ( 1.0 + ( - 3.0 / 2.0 + 1.0 / 3.0 * pz ) * pz ) * pz + + ( ( - 1.0 + 1.0 / 2.0 * pz ) * pz + ( pz + ( - pz + t ) * t ) * t ) * t ) + / ( t * t * t ); + } +// +// Region 7: Everywhere else. +// Series solution about infinity. +// + else + { + pz = log ( z ); + w = ( ( 1.0 + ( - 3.0 / 2.0 + 1.0 / 3.0 * pz ) * pz ) * pz + + ( ( - 1.0 + 1.0 / 2.0 * pz ) * pz + ( pz + ( - pz + z ) * z ) * z ) * z ) + / ( z * z * z ); + } +// +// Regularize if near branch cuts. +/// + if ( x <= - 1.0 + near && ( fabs ( ympi ) <= near || fabs ( yppi ) <= near ) ) + { + s = - 1.0; + if ( fabs ( ympi ) <= near ) + { +// +// Recompute ympi with directed rounding. +// + fesetround ( FE_UPWARD ); + ympi = y - pi; + + if ( ympi <= 0.0 ) + { + fesetround ( FE_DOWNWARD ); + ympi = y - pi; + } + + z = complex ( x, ympi ); +// +// Return rounding to default. +// + fesetround ( FE_TONEAREST ); + } + else + { +// +// Recompute yppi with directed rounding. +// + fesetround ( FE_UPWARD ); + yppi = y + pi; + + if ( yppi <= 0.0 ) + { + fesetround ( FE_DOWNWARD ); + yppi = y + pi; + } + + z = complex ( x, yppi ); +// +// Return rounding to default. +// + fesetround ( FE_TONEAREST ); + } + } +// +// Iteration one. +// + w = s * w; + r = z - s * w - log ( w ); + wp1 = s * w + 1.0; + e = r / wp1 * ( 2.0 * wp1 * ( wp1 + 2.0 / 3.0 * r ) - r ) + / ( 2.0 * wp1 * ( wp1 + 2.0 / 3.0 * r ) - 2.0 * r ); + w = w * ( 1.0 + e ); +// +// Iteration two. +// + if ( abs ( ( 2.0 * w * w - 8.0 * w - 1.0 ) * pow ( abs ( r ), 4.0 ) ) + >= TWOITERTOL * 72.0 * pow ( abs ( wp1 ), 6.0 ) ) + { + r = z - s * w - log ( w ); + wp1 = s * w + 1.0; + e = r / wp1 * ( 2.0 * wp1 * ( wp1 + 2.0 / 3.0 * r ) - r ) + / ( 2.0 * wp1 * ( wp1 + 2.0 / 3.0 * r ) - 2.0 * r ); + w = w * ( 1.0 + e ); + } +// +// Undo regularization. +// + w = s * w; +// +// Provide condition number estimate. +// + cond = z / ( 1.0 + w ); + + return 0; +} +//****************************************************************************80 + +//****************************************************************************80 + +inline complex wrightomega ( complex z ) + +//****************************************************************************80 +// +// Purpose: +// +// WRIGHTOMEGA is the simple routine for evaluating the Wright Omega function. +// +// Discussion: +// +// This function is called by: +// +// w = wrightomega ( z ) +// +// This function makes a call to the more powerful wrightomega_ext() function. +// +// Modified: +// +// 14 May 2016 +// +// Author: +// +// Piers Lawrence, Robert Corless, David Jeffrey +// +// Reference: +// +// Piers Lawrence, Robert Corless, David Jeffrey, +// Algorithm 917: Complex Double-Precision Evaluation of the Wright Omega +// Function, +// ACM Transactions on Mathematical Software, +// Volume 38, Number 3, Article 20, April 2012, 17 pages. +// +// Parameters: +// +// Input, complex Z, the argument. +// +// Output, complex WRIGHTOMEGA, the value of the Wright Omega +// function of Z. +// +{ + complex cond; + complex e; + complex r; + complex w; + + wrightomega_ext ( z, w, e, r, cond ); + + return w; +} + +inline float wrightomega ( float z ) +{ + return (float) std::real (wrightomega (std::complex { double (z), 0.0 })); +} +} diff --git a/test/src/wright_omega_approx_test.cpp b/test/src/wright_omega_approx_test.cpp new file mode 100644 index 0000000..091d92a --- /dev/null +++ b/test/src/wright_omega_approx_test.cpp @@ -0,0 +1,109 @@ +#include "test_helpers.hpp" +#include +#include + +#include +#include "reference/toms917.hpp" + +TEST_CASE ("Wright-Omega Approx Test") +{ +#if ! defined(WIN32) + const auto all_floats = test_helpers::all_32_bit_floats (-10.0f, 30.0f, 5.0e-3f); +#else + const auto all_floats = test_helpers::all_32_bit_floats (-10.0f, 30.0f, 5.0e-1f); +#endif + const auto y_exact = test_helpers::compute_all (all_floats, [] (auto x) + { return toms917::wrightomega (x); }); + + const auto test_approx = [&all_floats, &y_exact] (auto&& f_approx, float err_bound, float rel_err_bound, uint32_t ulp_bound) + { + const auto y_approx = test_helpers::compute_all (all_floats, f_approx); + + const auto error = test_helpers::compute_error (y_exact, y_approx); + const auto rel_error = test_helpers::compute_rel_error (y_exact, y_approx); + const auto ulp_error = test_helpers::compute_ulp_error (y_exact, y_approx); + + const auto max_error = test_helpers::abs_max (error); + const auto max_rel_error = test_helpers::abs_max (rel_error); + const auto max_ulp_error = *std::max_element (ulp_error.begin(), ulp_error.end()); + + std::cout << max_error << ", " << max_rel_error << ", " << max_ulp_error << std::endl; + REQUIRE (std::abs (max_error) < err_bound); + REQUIRE (std::abs (max_rel_error) < rel_err_bound); + if (ulp_bound > 0) + REQUIRE (max_ulp_error < ulp_bound); + }; + + SECTION ("Iter-3_Poly-3_LogExp-5") + { + test_approx ([] (auto x) + { return math_approx::wright_omega<3, 3, 5> (x); }, + 2.0e-6f, + 1.5e-6f, + 20); + } + SECTION ("Iter-3_Poly-3") + { + test_approx ([] (auto x) + { return math_approx::wright_omega<3, 3> (x); }, + 4.0e-6f, + 4.5e-6f, + 70); + } + SECTION ("Iter-2_Poly-5") + { + test_approx ([] (auto x) + { return math_approx::wright_omega<2, 5> (x); }, + 7.0e-6f, + 1.5e-4f, + 0); + } + SECTION ("Iter-2_Poly-3") + { + test_approx ([] (auto x) + { return math_approx::wright_omega<2, 3> (x); }, + 1.5e-5f, + 2.0e-4f, + 0); + } + SECTION ("Iter-2_Poly-3_LogExp-3") + { + test_approx ([] (auto x) + { return math_approx::wright_omega<2, 3, 3> (x); }, + 1.0e-4f, + 3.0e-4f, + 0); + } + SECTION ("Iter-1_Poly-5") + { + test_approx ([] (auto x) + { return math_approx::wright_omega<1, 5> (x); }, + 3.0e-3f, + 5.1e-2f, + 0); + } + SECTION ("Iter-1_Poly-3") + { + test_approx ([] (auto x) + { return math_approx::wright_omega<1, 3> (x); }, + 3.5e-3f, + 5.5e-2f, + 0); + } + SECTION ("Iter-0_Poly-5") + { + test_approx ([] (auto x) + { return math_approx::wright_omega<0, 5> (x); }, + 5.5e-2f, + 2.0f, + 0); + } + SECTION ("Iter-0_Poly-3") + { + test_approx ([] (auto x) + { return math_approx::wright_omega<0, 3> (x); }, + 6.0e-2f, + 2.0f, + 0); + } +} diff --git a/tools/bench/CMakeLists.txt b/tools/bench/CMakeLists.txt index 773d33c..c96dff6 100644 --- a/tools/bench/CMakeLists.txt +++ b/tools/bench/CMakeLists.txt @@ -19,3 +19,6 @@ target_link_libraries(pow_approx_bench PRIVATE benchmark::benchmark math_approx) add_executable(log_approx_bench log_bench.cpp) target_link_libraries(log_approx_bench PRIVATE benchmark::benchmark math_approx) + +add_executable(wright_omega_approx_bench wright_omega_bench.cpp) +target_link_libraries(wright_omega_approx_bench PRIVATE benchmark::benchmark math_approx) diff --git a/tools/bench/wright_omega_bench.cpp b/tools/bench/wright_omega_bench.cpp new file mode 100644 index 0000000..9279252 --- /dev/null +++ b/tools/bench/wright_omega_bench.cpp @@ -0,0 +1,67 @@ +#include +#include + +#include "../../test/src/reference/toms917.hpp" + +static constexpr size_t N = 2000; +const auto data = [] +{ + std::vector x; + x.resize (N, 0.0f); + for (size_t i = 0; i < N; ++i) + x[i] = -10.0f + 40.0f * (float) i / (float) N; + return x; +}(); + +#define WO_BENCH(name, func) \ +void name (benchmark::State& state) \ +{ \ +for (auto _ : state) \ +{ \ +for (auto& x : data) \ +{ \ +auto y = func (x); \ +benchmark::DoNotOptimize (y); \ +} \ +} \ +} \ +BENCHMARK (name); +WO_BENCH (wright_omega_toms917, toms917::wrightomega) +WO_BENCH (wright_omega_iter3_poly3_logexp5, (math_approx::wright_omega<3, 3, 5>)) +WO_BENCH (wright_omega_iter3_poly3, (math_approx::wright_omega<3, 3>)) +WO_BENCH (wright_omega_iter2_poly5, (math_approx::wright_omega<2, 5>)) +WO_BENCH (wright_omega_iter2_poly3, (math_approx::wright_omega<2, 3>)) +WO_BENCH (wright_omega_iter2_poly3_logexp3, (math_approx::wright_omega<2, 3, 3>)) +WO_BENCH (wright_omega_iter1_poly5, (math_approx::wright_omega<1, 5>)) +WO_BENCH (wright_omega_iter1_poly3, (math_approx::wright_omega<1, 3>)) +WO_BENCH (wright_omega_iter0_poly5, (math_approx::wright_omega<0, 5>)) +WO_BENCH (wright_omega_iter0_poly3, (math_approx::wright_omega<0, 3>)) +WO_BENCH (wright_omega_dangelo2, (math_approx::wright_omega_dangelo<2>)) +WO_BENCH (wright_omega_dangelo1, (math_approx::wright_omega_dangelo<1>)) +WO_BENCH (wright_omega_dangelo0, (math_approx::wright_omega_dangelo<0>)) + +#define WO_SIMD_BENCH(name, func) \ +void name (benchmark::State& state) \ +{ \ +for (auto _ : state) \ +{ \ +for (auto& x : data) \ +{ \ +auto y = func (xsimd::broadcast (x)); \ +static_assert (std::is_same_v, decltype(y)>); \ +benchmark::DoNotOptimize (y); \ +} \ +} \ +} \ +BENCHMARK (name); +WO_SIMD_BENCH (wright_omega_simd_iter3_poly3_logexp5, (math_approx::wright_omega<3, 3, 5>)) +WO_SIMD_BENCH (wright_omega_simd_iter3_poly3, (math_approx::wright_omega<3, 3>)) +WO_SIMD_BENCH (wright_omega_simd_iter2_poly5, (math_approx::wright_omega<2, 5>)) +WO_SIMD_BENCH (wright_omega_simd_iter2_poly3, (math_approx::wright_omega<2, 3>)) +WO_SIMD_BENCH (wright_omega_simd_iter2_poly3_logexp3, (math_approx::wright_omega<2, 3, 3>)) +WO_SIMD_BENCH (wright_omega_simd_iter1_poly5, (math_approx::wright_omega<1, 5>)) +WO_SIMD_BENCH (wright_omega_simd_iter1_poly3, (math_approx::wright_omega<1, 3>)) +WO_SIMD_BENCH (wright_omega_simd_iter0_poly5, (math_approx::wright_omega<0, 5>)) +WO_SIMD_BENCH (wright_omega_simd_iter0_poly3, (math_approx::wright_omega<0, 3>)) + +BENCHMARK_MAIN(); diff --git a/tools/plotter/plotter.cpp b/tools/plotter/plotter.cpp index 879533b..9d25d44 100644 --- a/tools/plotter/plotter.cpp +++ b/tools/plotter/plotter.cpp @@ -7,6 +7,8 @@ namespace plt = matplotlibcpp; #include "../../test/src/test_helpers.hpp" +#include "../../test/src/reference/toms917.hpp" +#include "../../test/src/reference/dangelo_omega.hpp" #include template @@ -15,9 +17,9 @@ void plot_error (std::span all_floats, F_Approx&& f_approx, const std::string& name) { - const auto y_approx = test_helpers::compute_all (all_floats, f_approx); - const auto error = test_helpers::compute_error (y_exact, y_approx); - std::cout << "Max Error: " << test_helpers::abs_max (error) << std::endl; + const auto y_approx = test_helpers::compute_all (all_floats, f_approx); + const auto error = test_helpers::compute_error (y_exact, y_approx); + std::cout << "Max Error: " << test_helpers::abs_max (error) << std::endl; plt::named_plot (name, all_floats, error); } @@ -27,9 +29,9 @@ void plot_rel_error (std::span all_floats, F_Approx&& f_approx, const std::string& name) { - const auto y_approx = test_helpers::compute_all (all_floats, f_approx); - const auto rel_error = test_helpers::compute_rel_error (y_exact, y_approx); - std::cout << "Max Relative Error: " << test_helpers::abs_max (rel_error) << std::endl; + const auto y_approx = test_helpers::compute_all (all_floats, f_approx); + const auto rel_error = test_helpers::compute_rel_error (y_exact, y_approx); + std::cout << "Max Relative Error: " << test_helpers::abs_max (rel_error) << std::endl; plt::named_plot (name, all_floats, rel_error); } @@ -39,7 +41,7 @@ void plot_ulp_error (std::span all_floats, F_Approx&& f_approx, const std::string& name) { - const auto y_approx = test_helpers::compute_all (all_floats, f_approx); + const auto y_approx = test_helpers::compute_all (all_floats, f_approx); const auto ulp_error = test_helpers::compute_ulp_error (y_exact, y_approx); std::cout << "Max Relative Error: " << *std::max_element (ulp_error.begin(), ulp_error.end()) << std::endl; plt::named_plot (name, all_floats, std::vector { ulp_error.begin(), ulp_error.end() }); @@ -59,16 +61,18 @@ void plot_function (std::span all_floats, int main() { plt::figure(); - const auto range = std::make_pair (0.01f, 10.0f); - static constexpr auto tol = 1.0e-2f; + const auto range = std::make_pair (-10.0f, 30.0f); + static constexpr auto tol = 1.0e-1f; 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::log)); + const auto y_exact = test_helpers::compute_all (all_floats, FLOAT_FUNC(toms917::wrightomega)); - // 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"); + // plot_error (all_floats, y_exact, FLOAT_FUNC((math_approx::wright_omega<0>)), "W-O 0-3"); + // plot_error (all_floats, y_exact, FLOAT_FUNC((math_approx::wright_omega<0, 5>)), "W-O 0-5"); + plot_error (all_floats, y_exact, FLOAT_FUNC((math_approx::wright_omega<1>)), "W-O 1-3"); + // plot_error (all_floats, y_exact, FLOAT_FUNC((math_approx::wright_omega_dangelo<0>)), "W-O D'Angelo 0"); + plot_error (all_floats, y_exact, FLOAT_FUNC((math_approx::wright_omega_dangelo<1>)), "W-O D'Angelo 1"); + plot_error (all_floats, y_exact, FLOAT_FUNC((math_approx::wright_omega_dangelo<2>)), "W-O D'Angelo 2"); plt::legend ({ { "loc", "upper right" } }); plt::xlim (range.first, range.second);