diff options
author | Daniel Baumann <daniel.baumann@progress-linux.org> | 2024-04-27 18:24:20 +0000 |
---|---|---|
committer | Daniel Baumann <daniel.baumann@progress-linux.org> | 2024-04-27 18:24:20 +0000 |
commit | 483eb2f56657e8e7f419ab1a4fab8dce9ade8609 (patch) | |
tree | e5d88d25d870d5dedacb6bbdbe2a966086a0a5cf /src/boost/libs/math/test/test_autodiff_2.cpp | |
parent | Initial commit. (diff) | |
download | ceph-upstream.tar.xz ceph-upstream.zip |
Adding upstream version 14.2.21.upstream/14.2.21upstream
Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>
Diffstat (limited to 'src/boost/libs/math/test/test_autodiff_2.cpp')
-rw-r--r-- | src/boost/libs/math/test/test_autodiff_2.cpp | 512 |
1 files changed, 512 insertions, 0 deletions
diff --git a/src/boost/libs/math/test/test_autodiff_2.cpp b/src/boost/libs/math/test/test_autodiff_2.cpp new file mode 100644 index 00000000..378a9bb1 --- /dev/null +++ b/src/boost/libs/math/test/test_autodiff_2.cpp @@ -0,0 +1,512 @@ +// Copyright Matthew Pulver 2018 - 2019. +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// https://www.boost.org/LICENSE_1_0.txt) + +#include "test_autodiff.hpp" + +BOOST_AUTO_TEST_SUITE(test_autodiff_2) + +BOOST_AUTO_TEST_CASE_TEMPLATE(one_over_one_plus_x_squared, T, all_float_types) { + constexpr std::size_t m = 4; + const T cx(1); + auto f = make_fvar<T, m>(cx); + // f = 1 / ((f *= f) += 1); + f *= f; + f += T(1); + f = f.inverse(); + BOOST_CHECK_EQUAL(f.derivative(0u), 0.5); + BOOST_CHECK_EQUAL(f.derivative(1u), -0.5); + BOOST_CHECK_EQUAL(f.derivative(2u), 0.5); + BOOST_CHECK_EQUAL(f.derivative(3u), 0); + BOOST_CHECK_EQUAL(f.derivative(4u), -3); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(exp_test, T, all_float_types) { + using std::exp; + constexpr std::size_t m = 4; + const T cx = 2.0; + const auto x = make_fvar<T, m>(cx); + auto y = exp(x); + for (auto i : boost::irange(m + 1)) { + // std::cout.precision(100); + // std::cout << "y.derivative("<<i<<") = " << y.derivative(i) << ", + // std::exp(cx) = " << std::exp(cx) << std::endl; + BOOST_CHECK_CLOSE_FRACTION(y.derivative(i), exp(cx), + std::numeric_limits<T>::epsilon()); + } +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(pow, T, bin_float_types) { + const T eps = 201 * std::numeric_limits<T>::epsilon(); // percent + using std::log; + using std::pow; + constexpr std::size_t m = 5; + constexpr std::size_t n = 4; + const T cx = 2.0; + const T cy = 3.0; + const auto x = make_fvar<T, m>(cx); + const auto y = make_fvar<T, m, n>(cy); + auto z0 = pow(x, cy); + BOOST_CHECK_EQUAL(z0.derivative(0u), pow(cx, cy)); + BOOST_CHECK_EQUAL(z0.derivative(1u), cy * pow(cx, cy - 1)); + BOOST_CHECK_EQUAL(z0.derivative(2u), cy * (cy - 1) * pow(cx, cy - 2)); + BOOST_CHECK_EQUAL(z0.derivative(3u), + cy * (cy - 1) * (cy - 2) * pow(cx, cy - 3)); + BOOST_CHECK_EQUAL(z0.derivative(4u), 0u); + BOOST_CHECK_EQUAL(z0.derivative(5u), 0u); + auto z1 = pow(cx, y); + BOOST_CHECK_CLOSE(z1.derivative(0u, 0u), pow(cx, cy), eps); + for (auto j : boost::irange(std::size_t(1), n + 1)) { + BOOST_CHECK_CLOSE(z1.derivative(0u, j), pow(log(cx), j) * pow(cx, cy), eps); + } + + for (auto i : boost::irange(std::size_t(1), m + 1)) { + for (auto j : boost::irange(n + 1)) { + BOOST_CHECK_EQUAL(z1.derivative(i, j), 0); + } + } + + const auto z2 = pow(x, y); + + for (auto j : boost::irange(n + 1)) { + BOOST_CHECK_CLOSE(z2.derivative(0u, j), pow(cx, cy) * pow(log(cx), j), eps); + } + for (auto j : boost::irange(n + 1)) { + BOOST_CHECK_CLOSE(z2.derivative(1u, j), + pow(cx, cy - 1) * pow(log(cx), static_cast<int>(j) - 1) * + (cy * log(cx) + j), + eps); + } + BOOST_CHECK_CLOSE(z2.derivative(2u, 0u), pow(cx, cy - 2) * cy * (cy - 1), + eps); + BOOST_CHECK_CLOSE(z2.derivative(2u, 1u), + pow(cx, cy - 2) * (cy * (cy - 1) * log(cx) + 2 * cy - 1), + eps); + for (auto j : boost::irange(std::size_t(2u), n + 1)) { + BOOST_CHECK_CLOSE(z2.derivative(2u, j), + pow(cx, cy - 2) * pow(log(cx), j - 2) * + (j * (2 * cy - 1) * log(cx) + (j - 1) * j + + (cy - 1) * cy * pow(log(cx), 2)), + eps); + } + BOOST_CHECK_CLOSE(z2.derivative(2u, 4u), + pow(cx, cy - 2) * pow(log(cx), 2) * + (4 * (2 * cy - 1) * log(cx) + (4 - 1) * 4 + + (cy - 1) * cy * pow(log(cx), 2)), + eps); +} + +// TODO Tests around x=0 or y=0: pow(x,y) +BOOST_AUTO_TEST_CASE_TEMPLATE(pow2, T, bin_float_types) { + const T eps = 4000 * std::numeric_limits<T>::epsilon(); // percent + using std::pow; + constexpr std::size_t m = 5; + constexpr std::size_t n = 5; + const T cx = 2; + const T cy = 5 / 2.0; + const auto x = make_fvar<T, m>(cx); + const auto y = make_fvar<T, 0, n>(cy); + const auto z = pow(x, y); + using namespace boost::math::constants; + // Mathematica: Export["pow.csv", Flatten@Table[ Simplify@D[x^y,{x,i},{y,j}] + // /. {x->2, y->5/2}, + // { i, 0, 5 }, { j, 0, 5 } ] ] + // sed -rf pow.sed < pow.csv + // with pow.sed script: + // s/Log\[2\]\^([0-9]+)/pow(ln_two<T>(),\1)/g + // s/Log\[2\]/ln_two<T>()/g + // s/Sqrt\[2\]/root_two<T>()/g + // s/[0-9]\/[0-9]+/\0.0/g + // s/^"/static_cast<T>(/ + // s/"$/),/ + const T mathematica[]{ + static_cast<T>(4 * root_two<T>()), + static_cast<T>(4 * root_two<T>() * ln_two<T>()), + static_cast<T>(4 * root_two<T>() * pow(ln_two<T>(), 2)), + static_cast<T>(4 * root_two<T>() * pow(ln_two<T>(), 3)), + static_cast<T>(4 * root_two<T>() * pow(ln_two<T>(), 4)), + static_cast<T>(4 * root_two<T>() * pow(ln_two<T>(), 5)), + static_cast<T>(5 * root_two<T>()), + static_cast<T>(2 * root_two<T>() * (1 + (5 * ln_two<T>()) / 2)), + static_cast<T>(2 * root_two<T>() * ln_two<T>() * + (2 + (5 * ln_two<T>()) / 2)), + static_cast<T>(2 * root_two<T>() * pow(ln_two<T>(), 2) * + (3 + (5 * ln_two<T>()) / 2)), + static_cast<T>(2 * root_two<T>() * pow(ln_two<T>(), 3) * + (4 + (5 * ln_two<T>()) / 2)), + static_cast<T>(2 * root_two<T>() * pow(ln_two<T>(), 4) * + (5 + (5 * ln_two<T>()) / 2)), + static_cast<T>(15 / (2 * root_two<T>())), + static_cast<T>(root_two<T>() * (4 + (15 * ln_two<T>()) / 4)), + static_cast<T>(root_two<T>() * + (2 + 8 * ln_two<T>() + (15 * pow(ln_two<T>(), 2)) / 4)), + static_cast<T>(root_two<T>() * ln_two<T>() * + (6 + 12 * ln_two<T>() + (15 * pow(ln_two<T>(), 2)) / 4)), + static_cast<T>(root_two<T>() * pow(ln_two<T>(), 2) * + (12 + 16 * ln_two<T>() + (15 * pow(ln_two<T>(), 2)) / 4)), + static_cast<T>(root_two<T>() * pow(ln_two<T>(), 3) * + (20 + 20 * ln_two<T>() + (15 * pow(ln_two<T>(), 2)) / 4)), + static_cast<T>(15 / (8 * root_two<T>())), + static_cast<T>((23 / 4.0 + (15 * ln_two<T>()) / 8) / root_two<T>()), + static_cast<T>( + (9 + (23 * ln_two<T>()) / 2 + (15 * pow(ln_two<T>(), 2)) / 8) / + root_two<T>()), + static_cast<T>((6 + 27 * ln_two<T>() + (69 * pow(ln_two<T>(), 2)) / 4 + + (15 * pow(ln_two<T>(), 3)) / 8) / + root_two<T>()), + static_cast<T>( + (ln_two<T>() * (24 + 54 * ln_two<T>() + 23 * pow(ln_two<T>(), 2) + + (15 * pow(ln_two<T>(), 3)) / 8)) / + root_two<T>()), + static_cast<T>((pow(ln_two<T>(), 2) * + (60 + 90 * ln_two<T>() + (115 * pow(ln_two<T>(), 2)) / 4 + + (15 * pow(ln_two<T>(), 3)) / 8)) / + root_two<T>()), + static_cast<T>(-15 / (32 * root_two<T>())), + static_cast<T>((-1 - (15 * ln_two<T>()) / 16) / (2 * root_two<T>())), + static_cast<T>((7 - 2 * ln_two<T>() - (15 * pow(ln_two<T>(), 2)) / 16) / + (2 * root_two<T>())), + static_cast<T>((24 + 21 * ln_two<T>() - 3 * pow(ln_two<T>(), 2) - + (15 * pow(ln_two<T>(), 3)) / 16) / + (2 * root_two<T>())), + static_cast<T>((24 + 96 * ln_two<T>() + 42 * pow(ln_two<T>(), 2) - + 4 * pow(ln_two<T>(), 3) - + (15 * pow(ln_two<T>(), 4)) / 16) / + (2 * root_two<T>())), + static_cast<T>( + (ln_two<T>() * + (120 + 240 * ln_two<T>() + 70 * pow(ln_two<T>(), 2) - + 5 * pow(ln_two<T>(), 3) - (15 * pow(ln_two<T>(), 4)) / 16)) / + (2 * root_two<T>())), + static_cast<T>(45 / (128 * root_two<T>())), + static_cast<T>((9 / 16.0 + (45 * ln_two<T>()) / 32) / + (4 * root_two<T>())), + static_cast<T>((-25 / 2.0 + (9 * ln_two<T>()) / 8 + + (45 * pow(ln_two<T>(), 2)) / 32) / + (4 * root_two<T>())), + static_cast<T>((-15 - (75 * ln_two<T>()) / 2 + + (27 * pow(ln_two<T>(), 2)) / 16 + + (45 * pow(ln_two<T>(), 3)) / 32) / + (4 * root_two<T>())), + static_cast<T>((60 - 60 * ln_two<T>() - 75 * pow(ln_two<T>(), 2) + + (9 * pow(ln_two<T>(), 3)) / 4 + + (45 * pow(ln_two<T>(), 4)) / 32) / + (4 * root_two<T>())), + static_cast<T>((120 + 300 * ln_two<T>() - 150 * pow(ln_two<T>(), 2) - + 125 * pow(ln_two<T>(), 3) + + (45 * pow(ln_two<T>(), 4)) / 16 + + (45 * pow(ln_two<T>(), 5)) / 32) / + (4 * root_two<T>()))}; + std::size_t k = 0; + for (auto i : boost::irange(m + 1)) { + for (auto j : boost::irange(n + 1)) { + BOOST_CHECK_CLOSE(z.derivative(i, j), mathematica[k++], eps); + } + } +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(sqrt_test, T, all_float_types) { + using std::pow; + using std::sqrt; + constexpr std::size_t m = 5; + const T cx = 4.0; + auto x = make_fvar<T, m>(cx); + auto y = sqrt(x); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(0u), sqrt(cx), + std::numeric_limits<T>::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(1u), 0.5 * pow(cx, -0.5), + std::numeric_limits<T>::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(2u), -0.5 * 0.5 * pow(cx, -1.5), + std::numeric_limits<T>::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(3u), 0.5 * 0.5 * 1.5 * pow(cx, -2.5), + std::numeric_limits<T>::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(4u), + -0.5 * 0.5 * 1.5 * 2.5 * pow(cx, -3.5), + std::numeric_limits<T>::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(5u), + 0.5 * 0.5 * 1.5 * 2.5 * 3.5 * pow(cx, -4.5), + std::numeric_limits<T>::epsilon()); + x = make_fvar<T, m>(0); + y = sqrt(x); + // std::cout << "sqrt(0) = " << y << std::endl; // (0,inf,-inf,inf,-inf,inf) + BOOST_CHECK_EQUAL(y.derivative(0u), 0); + for (auto i : boost::irange(std::size_t(1), m + 1)) { + BOOST_CHECK_EQUAL(y.derivative(i), (i % 2 == 1 ? 1 : -1) * + std::numeric_limits<T>::infinity()); + } +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(log_test, T, all_float_types) { + using std::log; + using std::pow; + constexpr std::size_t m = 5; + const T cx = 2.0; + auto x = make_fvar<T, m>(cx); + auto y = log(x); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(0u), log(cx), + std::numeric_limits<T>::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(1u), 1 / cx, + std::numeric_limits<T>::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(2u), -1 / pow(cx, 2), + std::numeric_limits<T>::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(3u), 2 / pow(cx, 3), + std::numeric_limits<T>::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(4u), -6 / pow(cx, 4), + std::numeric_limits<T>::epsilon()); + BOOST_CHECK_CLOSE_FRACTION(y.derivative(5u), 24 / pow(cx, 5), + std::numeric_limits<T>::epsilon()); + x = make_fvar<T, m>(0); + y = log(x); + // std::cout << "log(0) = " << y << std::endl; // log(0) = + // depth(1)(-inf,inf,-inf,inf,-inf,inf) + for (auto i : boost::irange(m + 1)) { + BOOST_CHECK_EQUAL(y.derivative(i), (i % 2 == 1 ? 1 : -1) * + std::numeric_limits<T>::infinity()); + } +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(ylogx, T, all_float_types) { + using std::log; + using std::pow; + const T eps = 100 * std::numeric_limits<T>::epsilon(); // percent + constexpr std::size_t m = 5; + constexpr std::size_t n = 4; + const T cx = 2.0; + const T cy = 3.0; + const auto x = make_fvar<T, m>(cx); + const auto y = make_fvar<T, m, n>(cy); + auto z = y * log(x); + BOOST_CHECK_EQUAL(z.derivative(0u, 0u), cy * log(cx)); + BOOST_CHECK_EQUAL(z.derivative(0u, 1u), log(cx)); + BOOST_CHECK_EQUAL(z.derivative(0u, 2u), 0); + BOOST_CHECK_EQUAL(z.derivative(0u, 3u), 0); + BOOST_CHECK_EQUAL(z.derivative(0u, 4u), 0); + for (auto i : boost::irange(1u, static_cast<unsigned>(m + 1))) { + BOOST_CHECK_CLOSE(z.derivative(i, 0u), + pow(-1, i - 1) * boost::math::factorial<T>(i - 1) * cy / + pow(cx, i), + eps); + BOOST_CHECK_CLOSE( + z.derivative(i, 1u), + pow(-1, i - 1) * boost::math::factorial<T>(i - 1) / pow(cx, i), eps); + for (auto j : boost::irange(std::size_t(2), n + 1)) { + BOOST_CHECK_EQUAL(z.derivative(i, j), 0u); + } + } + auto z1 = exp(z); + // RHS is confirmed by + // https://www.wolframalpha.com/input/?i=D%5Bx%5Ey,%7Bx,2%7D,%7By,4%7D%5D+%2F.+%7Bx-%3E2.0,+y-%3E3.0%7D + BOOST_CHECK_CLOSE(z1.derivative(2u, 4u), + pow(cx, cy - 2) * pow(log(cx), 2) * + (4 * (2 * cy - 1) * log(cx) + (4 - 1) * 4 + + (cy - 1) * cy * pow(log(cx), 2)), + eps); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(frexp_test, T, all_float_types) { + using std::exp2; + using std::frexp; + constexpr std::size_t m = 3; + const T cx = 3.5; + const auto x = make_fvar<T, m>(cx); + int exp, testexp; + auto y = frexp(x, &exp); + BOOST_CHECK_EQUAL(y.derivative(0u), frexp(cx, &testexp)); + BOOST_CHECK_EQUAL(exp, testexp); + BOOST_CHECK_EQUAL(y.derivative(1u), exp2(-exp)); + BOOST_CHECK_EQUAL(y.derivative(2u), 0); + BOOST_CHECK_EQUAL(y.derivative(3u), 0); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(ldexp_test, T, all_float_types) { + BOOST_MATH_STD_USING + using boost::multiprecision::ldexp; + constexpr auto m = 3u; + const T cx = 3.5; + const auto x = make_fvar<T, m>(cx); + constexpr auto exponent = 3; + auto y = ldexp(x, exponent); + BOOST_CHECK_EQUAL(y.derivative(0u), ldexp(cx, exponent)); + BOOST_CHECK_EQUAL(y.derivative(1u), exp2(exponent)); + BOOST_CHECK_EQUAL(y.derivative(2u), 0); + BOOST_CHECK_EQUAL(y.derivative(3u), 0); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(cos_and_sin, T, bin_float_types) { + using std::cos; + using std::sin; + const T eps = 200 * std::numeric_limits<T>::epsilon(); // percent + constexpr std::size_t m = 5; + const T cx = boost::math::constants::third_pi<T>(); + const auto x = make_fvar<T, m>(cx); + auto cos5 = cos(x); + BOOST_CHECK_CLOSE(cos5.derivative(0u), cos(cx), eps); + BOOST_CHECK_CLOSE(cos5.derivative(1u), -sin(cx), eps); + BOOST_CHECK_CLOSE(cos5.derivative(2u), -cos(cx), eps); + BOOST_CHECK_CLOSE(cos5.derivative(3u), sin(cx), eps); + BOOST_CHECK_CLOSE(cos5.derivative(4u), cos(cx), eps); + BOOST_CHECK_CLOSE(cos5.derivative(5u), -sin(cx), eps); + auto sin5 = sin(x); + BOOST_CHECK_CLOSE(sin5.derivative(0u), sin(cx), eps); + BOOST_CHECK_CLOSE(sin5.derivative(1u), cos(cx), eps); + BOOST_CHECK_CLOSE(sin5.derivative(2u), -sin(cx), eps); + BOOST_CHECK_CLOSE(sin5.derivative(3u), -cos(cx), eps); + BOOST_CHECK_CLOSE(sin5.derivative(4u), sin(cx), eps); + BOOST_CHECK_CLOSE(sin5.derivative(5u), cos(cx), eps); + // Test Order = 0 for codecov + auto cos0 = cos(make_fvar<T, 0>(cx)); + BOOST_CHECK_CLOSE(cos0.derivative(0u), cos(cx), eps); + auto sin0 = sin(make_fvar<T, 0>(cx)); + BOOST_CHECK_CLOSE(sin0.derivative(0u), sin(cx), eps); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(acos_test, T, bin_float_types) { + const T eps = 300 * std::numeric_limits<T>::epsilon(); // percent + using std::acos; + using std::pow; + using std::sqrt; + constexpr std::size_t m = 5; + const T cx = 0.5; + auto x = make_fvar<T, m>(cx); + auto y = acos(x); + BOOST_CHECK_CLOSE(y.derivative(0u), acos(cx), eps); + BOOST_CHECK_CLOSE(y.derivative(1u), -1 / sqrt(1 - cx * cx), eps); + BOOST_CHECK_CLOSE(y.derivative(2u), -cx / pow(1 - cx * cx, 1.5), eps); + BOOST_CHECK_CLOSE(y.derivative(3u), + -(2 * cx * cx + 1) / pow(1 - cx * cx, 2.5), eps); + BOOST_CHECK_CLOSE(y.derivative(4u), + -3 * cx * (2 * cx * cx + 3) / pow(1 - cx * cx, 3.5), eps); + BOOST_CHECK_CLOSE(y.derivative(5u), + -(24 * (cx * cx + 3) * cx * cx + 9) / pow(1 - cx * cx, 4.5), + eps); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(acosh_test, T, bin_float_types) { + const T eps = 300 * std::numeric_limits<T>::epsilon(); // percent + using boost::math::acosh; + constexpr std::size_t m = 5; + const T cx = 2; + auto x = make_fvar<T, m>(cx); + auto y = acosh(x); + // BOOST_CHECK_EQUAL(y.derivative(0) == acosh(cx)); // FAILS! acosh(2) is + // overloaded for integral types + BOOST_CHECK_CLOSE(y.derivative(0u), acosh(static_cast<T>(x)), eps); + BOOST_CHECK_CLOSE(y.derivative(1u), + 1 / boost::math::constants::root_three<T>(), eps); + BOOST_CHECK_CLOSE(y.derivative(2u), + -2 / (3 * boost::math::constants::root_three<T>()), eps); + BOOST_CHECK_CLOSE(y.derivative(3u), + 1 / boost::math::constants::root_three<T>(), eps); + BOOST_CHECK_CLOSE(y.derivative(4u), + -22 / (9 * boost::math::constants::root_three<T>()), eps); + BOOST_CHECK_CLOSE(y.derivative(5u), + 227 / (27 * boost::math::constants::root_three<T>()), + 2 * eps); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(asin_test, T, bin_float_types) { + const T eps = 300 * std::numeric_limits<T>::epsilon(); // percent + using std::asin; + using std::pow; + using std::sqrt; + constexpr std::size_t m = 5; + const T cx = 0.5; + auto x = make_fvar<T, m>(cx); + auto y = asin(x); + BOOST_CHECK_CLOSE(y.derivative(0u), asin(static_cast<T>(x)), eps); + BOOST_CHECK_CLOSE(y.derivative(1u), 1 / sqrt(1 - cx * cx), eps); + BOOST_CHECK_CLOSE(y.derivative(2u), cx / pow(1 - cx * cx, 1.5), eps); + BOOST_CHECK_CLOSE(y.derivative(3u), (2 * cx * cx + 1) / pow(1 - cx * cx, 2.5), + eps); + BOOST_CHECK_CLOSE(y.derivative(4u), + 3 * cx * (2 * cx * cx + 3) / pow(1 - cx * cx, 3.5), eps); + BOOST_CHECK_CLOSE(y.derivative(5u), + (24 * (cx * cx + 3) * cx * cx + 9) / pow(1 - cx * cx, 4.5), + eps); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(asin_infinity, T, all_float_types) { + const T eps = 100 * std::numeric_limits<T>::epsilon(); // percent + constexpr std::size_t m = 5; + auto x = make_fvar<T, m>(1); + auto y = asin(x); + // std::cout << "asin(1) = " << y << std::endl; // + // depth(1)(1.5707963267949,inf,inf,-nan,-nan,-nan) + BOOST_CHECK_CLOSE(y.derivative(0u), boost::math::constants::half_pi<T>(), + eps); // MacOS is not exact + BOOST_CHECK_EQUAL(y.derivative(1u), std::numeric_limits<T>::infinity()); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(asin_derivative, T, bin_float_types) { + const T eps = 300 * std::numeric_limits<T>::epsilon(); // percent + using std::pow; + using std::sqrt; + constexpr std::size_t m = 4; + const T cx(0.5); + auto x = make_fvar<T, m>(cx); + auto y = T(1) - x * x; + BOOST_CHECK_EQUAL(y.derivative(0u), 1 - cx * cx); + BOOST_CHECK_EQUAL(y.derivative(1u), -2 * cx); + BOOST_CHECK_EQUAL(y.derivative(2u), -2); + BOOST_CHECK_EQUAL(y.derivative(3u), 0); + BOOST_CHECK_EQUAL(y.derivative(4u), 0); + y = sqrt(y); + BOOST_CHECK_EQUAL(y.derivative(0u), sqrt(1 - cx * cx)); + BOOST_CHECK_CLOSE(y.derivative(1u), -cx / sqrt(1 - cx * cx), eps); + BOOST_CHECK_CLOSE(y.derivative(2u), -1 / pow(1 - cx * cx, 1.5), eps); + BOOST_CHECK_CLOSE(y.derivative(3u), -3 * cx / pow(1 - cx * cx, 2.5), eps); + BOOST_CHECK_CLOSE(y.derivative(4u), + -(12 * cx * cx + 3) / pow(1 - cx * cx, 3.5), eps); + y = y.inverse(); // asin'(x) = 1 / sqrt(1-x*x). + BOOST_CHECK_CLOSE(y.derivative(0u), 1 / sqrt(1 - cx * cx), eps); + BOOST_CHECK_CLOSE(y.derivative(1u), cx / pow(1 - cx * cx, 1.5), eps); + BOOST_CHECK_CLOSE(y.derivative(2u), (2 * cx * cx + 1) / pow(1 - cx * cx, 2.5), + eps); + BOOST_CHECK_CLOSE(y.derivative(3u), + 3 * cx * (2 * cx * cx + 3) / pow(1 - cx * cx, 3.5), eps); + BOOST_CHECK_CLOSE(y.derivative(4u), + (24 * (cx * cx + 3) * cx * cx + 9) / pow(1 - cx * cx, 4.5), + eps); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(asinh_test, T, bin_float_types) { + const T eps = 300 * std::numeric_limits<T>::epsilon(); // percent + using boost::math::asinh; + constexpr std::size_t m = 5; + const T cx = 1; + auto x = make_fvar<T, m>(cx); + auto y = asinh(x); + BOOST_CHECK_CLOSE(y.derivative(0u), asinh(static_cast<T>(x)), eps); + BOOST_CHECK_CLOSE(y.derivative(1u), 1 / boost::math::constants::root_two<T>(), + eps); + BOOST_CHECK_CLOSE(y.derivative(2u), + -1 / (2 * boost::math::constants::root_two<T>()), eps); + BOOST_CHECK_CLOSE(y.derivative(3u), + 1 / (4 * boost::math::constants::root_two<T>()), eps); + BOOST_CHECK_CLOSE(y.derivative(4u), + 3 / (8 * boost::math::constants::root_two<T>()), eps); + BOOST_CHECK_CLOSE(y.derivative(5u), + -39 / (16 * boost::math::constants::root_two<T>()), eps); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(atan2_function, T, all_float_types) { + using test_constants = test_constants_t<T>; + static constexpr auto m = test_constants::order; + + test_detail::RandomSample<T> x_sampler{-2000, 2000}; + test_detail::RandomSample<T> y_sampler{-2000, 2000}; + + for (auto i : boost::irange(test_constants::n_samples)) { + std::ignore = i; + auto x = x_sampler.next(); + auto y = y_sampler.next(); + + auto autodiff_v = atan2(make_fvar<T, m>(x), make_fvar<T, m>(y)); + auto anchor_v = atan2(x, y); + BOOST_CHECK_CLOSE(autodiff_v, anchor_v, + 5000 * test_constants::pct_epsilon()); + } +} + +BOOST_AUTO_TEST_SUITE_END() |