Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Logs #376

Merged
merged 7 commits into from
Jan 15, 2024
Merged

Logs #376

Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions include/boost/decimal/cmath.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@
#include <boost/decimal/detail/cmath/log.hpp>
#include <boost/decimal/detail/cmath/log1p.hpp>
#include <boost/decimal/detail/cmath/log10.hpp>
#include <boost/decimal/detail/cmath/log2.hpp>
#include <boost/decimal/detail/cmath/modf.hpp>
#include <boost/decimal/detail/cmath/nearbyint.hpp>
#include <boost/decimal/detail/cmath/next.hpp>
Expand Down
97 changes: 58 additions & 39 deletions include/boost/decimal/detail/cmath/log.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,11 +12,67 @@
#include <boost/decimal/fwd.hpp>
#include <boost/decimal/detail/type_traits.hpp>
#include <boost/decimal/detail/concepts.hpp>
#include <boost/decimal/detail/emulated128.hpp>
#include <boost/decimal/numbers.hpp>

namespace boost {
namespace decimal {

namespace detail {

#if (defined(__clang__) && (__clang__ < 6))
# pragma clang diagnostic push
# pragma clang diagnostic ignored "-Wmissing-braces"
#endif

template <typename T>
static constexpr std::array<T, 11> small_coefficient_table {
// (1,) 12, 80, 448, 2304, 11264, 53248, 245760, 1114112, 4980736, 22020096, 96468992, ...
// See also Sloane's A058962 at: https://oeis.org/A058962

// See also
// Series[Log[(1 + (z/2))/(1 - (z/2))], {z, 0, 23}]
// Or at Wolfram Alpha: https://www.wolframalpha.com/input?i=Series%5BLog%5B%281+%2B+%28z%2F2%29%29%2F%281+-+%28z%2F2%29%29%5D%2C+%7Bz%2C+0%2C+21%7D%5D

// TODO(ckormanyos)
// Consider also Pade expansion of Log[1 + z/2].
// PadeApproximant[Log[1 + z/2], {z, 6, 6}]
// FullSimplify[%].

T { UINT64_C(833333333333333333), -18 - 1 }, // * z^3
T { UINT64_C(125000000000000000), -18 - 1 }, // * z^5
T { UINT64_C(223214285714285714), -18 - 2 }, // * z^7
T { UINT64_C(434027777777777778), -18 - 3 }, // * z^9
T { UINT64_C(887784090909090909), -18 - 4 }, // * z^11
T { UINT64_C(187800480769230769), -18 - 4 }, // * z^13
T { UINT64_C(406901041666666667), -18 - 5 }, // * z^15
T { UINT64_C(897575827205882353), -18 - 6 }, // * z^17
T { UINT64_C(200773540296052632), -18 - 6 }, // * z^19
T { UINT64_C(454130626860119048), -18 - 7 }, // * z^21
T { UINT64_C(103660251783288043), -18 - 7 } // * z^23
};

template <typename T>
static constexpr std::array<T, 11> large_coefficient_table {
T{uint128{UINT64_C(451750905202293), UINT64_C(9484758194528277842)}, -35}, // z^3
T{uint128{UINT64_C(67762635780344), UINT64_C(500376525493764096)}, -35}, // z^5
T{uint128{UINT64_C(121004706750614), UINT64_C(6164027816584450626)}, -36}, // z^7
T{uint128{UINT64_C(235286929792861), UINT64_C(3787056721709964394)}, -37}, // z^9
T{uint128{UINT64_C(481268720030852), UINT64_C(8584740752302634068)}, -38}, // z^11
T{uint128{UINT64_C(101806844621911), UINT64_C(1816002851448634124)}, -38}, // z^13
T{uint128{UINT64_C(220581496680807), UINT64_C(7009130190423632548)}, -40}, // z^15
T{uint128{UINT64_C(486576830913545), UINT64_C(12748560115094843630)}, -41}, // z^17
T{uint128{UINT64_C(108839554283293), UINT64_C(2123490654414259032)}, -41}, // z^19
T{uint128{UINT64_C(246184706116972), UINT64_C(9634423737622849438)}, -42}, // z^21
T{uint128{UINT64_C(56194335091917), UINT64_C(11823550152479764302)}, -42} // z^23
};

#if (defined(__clang__) && (__clang__ < 6))
# pragma clang diagnostic pop
#endif

} //namespace detail

template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
constexpr auto log(T x) noexcept -> std::enable_if_t<detail::is_decimal_floating_point_v<T>, T> // NOLINT(misc-no-recursion)
{
Expand Down Expand Up @@ -76,49 +132,12 @@ constexpr auto log(T x) noexcept -> std::enable_if_t<detail::is_decimal_floating
--exp2val;
}

using coefficient_array_type = std::array<T, static_cast<std::size_t>(UINT8_C(11))>;

#if (defined(__clang__) && (__clang__ < 6))
# pragma clang diagnostic push
# pragma clang diagnostic ignored "-Wmissing-braces"
#endif

constexpr auto coefficient_table =
coefficient_array_type
{
// (1,) 12, 80, 448, 2304, 11264, 53248, 245760, 1114112, 4980736, 22020096, 96468992, ...
// See also Sloane's A058962 at: https://oeis.org/A058962

// See also
// Series[Log[(1 + (z/2))/(1 - (z/2))], {z, 0, 23}]
// Or at Wolfram Alpha: https://www.wolframalpha.com/input?i=Series%5BLog%5B%281+%2B+%28z%2F2%29%29%2F%281+-+%28z%2F2%29%29%5D%2C+%7Bz%2C+0%2C+21%7D%5D

// TODO(ckormanyos)
// Consider also Pade expansion of Log[1 + z/2].
// PadeApproximant[Log[1 + z/2], {z, 6, 6}]
// FullSimplify[%].

T { UINT64_C(833333333333333333), -18 - 1 }, // * z^3
T { UINT64_C(125000000000000000), -18 - 1 }, // * z^5
T { UINT64_C(223214285714285714), -18 - 2 }, // * z^7
T { UINT64_C(434027777777777778), -18 - 3 }, // * z^9
T { UINT64_C(887784090909090909), -18 - 4 }, // * z^11
T { UINT64_C(187800480769230769), -18 - 4 }, // * z^13
T { UINT64_C(406901041666666667), -18 - 5 }, // * z^15
T { UINT64_C(897575827205882353), -18 - 6 }, // * z^17
T { UINT64_C(200773540296052632), -18 - 6 }, // * z^19
T { UINT64_C(454130626860119048), -18 - 7 }, // * z^21
T { UINT64_C(103660251783288043), -18 - 7 } // * z^23
};

#if (defined(__clang__) && (__clang__ < 6))
# pragma clang diagnostic pop
#endif

const auto s = (g - one) / (g + one);
const auto z = s + s;
const auto zsq = z * z;

auto& coefficient_table = std::is_same<T, decimal128>::value ? detail::large_coefficient_table<T> : detail::small_coefficient_table<T>;

auto rit = coefficient_table.crbegin() + static_cast<std::size_t>((sizeof(T) == 4U) ? 5U : 0U);

result = *rit;
Expand Down
30 changes: 30 additions & 0 deletions include/boost/decimal/detail/cmath/log2.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
// Copyright 2024 Matt Borland
// Copyright 2024 Christopher Kormanyos
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt

#ifndef BOOST_DECIMAL_DETAIL_CMATH_LOG2_HPP
#define BOOST_DECIMAL_DETAIL_CMATH_LOG2_HPP

#include <cmath>
#include <type_traits>

#include <boost/decimal/fwd.hpp> // NOLINT(llvm-include-order)
#include <boost/decimal/detail/type_traits.hpp>
#include <boost/decimal/detail/concepts.hpp>
#include <boost/decimal/numbers.hpp>
#include <boost/decimal/detail/cmath/log.hpp>

namespace boost {
namespace decimal {

template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
constexpr auto log2(T x) noexcept -> std::enable_if_t<detail::is_decimal_floating_point_v<T>, T>
{
return log(x) / numbers::ln2_v<T>;
}

} //namespace decimal
} //namespace boost

#endif //BOOST_DECIMAL_DETAIL_CMATH_LOG2_HPP
83 changes: 80 additions & 3 deletions test/test_cmath.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,10 @@
#include <cmath>


#if !defined(BOOST_DECIMAL_REDUCE_TEST_DEPTH)
static constexpr auto N = static_cast<std::size_t>(512U); // Number of trials
#if !defined(BOOST_DECIMAL_REDUCE_TEST_DEPTH) && !defined(_WIN32)
static constexpr auto N = static_cast<std::size_t>(128U); // Number of trials
#else
static constexpr auto N = static_cast<std::size_t>(512U >> 4U); // Number of trials
static constexpr auto N = static_cast<std::size_t>(128U >> 4U); // Number of trials
#endif

static std::mt19937_64 rng(42);
Expand Down Expand Up @@ -1556,6 +1556,72 @@ void test_nan()
}
#endif

template <typename Dec>
void test_log2()
{
std::uniform_real_distribution<float> dist(-0.5F, 0.5F);

constexpr auto max_iter {std::is_same<Dec, decimal128>::value ? N / 4 : N};
for (std::size_t n {}; n < max_iter; ++n)
{
const auto val1 {dist(rng)};
Dec d1 {val1};

auto ret_val {std::log2(val1)};
auto ret_dec {static_cast<float>(log2(d1))};

if (!std::isfinite(ret_val) && !std::isfinite(ret_dec))
{
continue;
}

const auto distance {std::fabs(boost::math::float_distance(ret_val, ret_dec))};
if (!BOOST_TEST(distance < 100))
{
// LCOV_EXCL_START
std::cerr << "Val 1: " << val1
<< "\nDec 1: " << d1
<< "\nRet val: " << ret_val
<< "\nRet dec: " << ret_dec
<< "\nDist: " << distance << std::endl;
// LCOV_EXCL_STOP
}
}
}

template <typename Dec>
void test_log10()
{
std::uniform_real_distribution<float> dist(-0.5F, 0.5F);

constexpr auto max_iter {std::is_same<Dec, decimal128>::value ? N / 4 : N};
for (std::size_t n {}; n < max_iter; ++n)
{
const auto val1 {dist(rng)};
Dec d1 {val1};

auto ret_val {std::log10(val1)};
auto ret_dec {static_cast<float>(log10(d1))};

if (!std::isfinite(ret_val) && !std::isfinite(ret_dec))
{
continue;
}

const auto distance {std::fabs(boost::math::float_distance(ret_val, ret_dec))};
if (!BOOST_TEST(distance < 100))
{
// LCOV_EXCL_START
std::cerr << "Val 1: " << val1
<< "\nDec 1: " << d1
<< "\nRet val: " << ret_val
<< "\nRet dec: " << ret_dec
<< "\nDist: " << distance << std::endl;
// LCOV_EXCL_STOP
}
}
}

int main()
{
test_fmax<decimal32>();
Expand Down Expand Up @@ -1686,5 +1752,16 @@ int main()
test_nan<decimal128>();
#endif

test_log2<decimal32>();
test_log2<decimal64>();

test_log10<decimal32>();
test_log10<decimal64>();

#if !defined(BOOST_DECIMAL_REDUCE_TEST_DEPTH)
test_log2<decimal128>();
test_log10<decimal128>();
#endif

return boost::report_errors();
}