From 65f6653abb4ae8049f09d426d0bab1be32cfec99 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Fri, 24 May 2024 06:12:08 +0200 Subject: [PATCH 01/11] Make Riemann-zeta skeleton --- include/boost/decimal/cmath.hpp | 1 + .../decimal/detail/cmath/riemann_zeta.hpp | 85 +++++++++++++++++++ test/Jamfile | 1 + test/test_cbrt.cpp | 11 +-- test/test_zeta.cpp | 61 +++++++++++++ 5 files changed, 151 insertions(+), 8 deletions(-) create mode 100644 include/boost/decimal/detail/cmath/riemann_zeta.hpp create mode 100644 test/test_zeta.cpp diff --git a/include/boost/decimal/cmath.hpp b/include/boost/decimal/cmath.hpp index acc5e9ffe..3011fb495 100644 --- a/include/boost/decimal/cmath.hpp +++ b/include/boost/decimal/cmath.hpp @@ -45,6 +45,7 @@ #include #include #include +#include #include #include #include diff --git a/include/boost/decimal/detail/cmath/riemann_zeta.hpp b/include/boost/decimal/detail/cmath/riemann_zeta.hpp new file mode 100644 index 000000000..8ed02b891 --- /dev/null +++ b/include/boost/decimal/detail/cmath/riemann_zeta.hpp @@ -0,0 +1,85 @@ +// 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_RIEMANN_ZETA_HPP +#define BOOST_DECIMAL_DETAIL_CMATH_RIEMANN_ZETA_HPP + +#include +#include +#include +#include +#include + +#ifndef BOOST_DECIMAL_BUILD_MODULE +#include +#include +#endif + +namespace boost { +namespace decimal { + +namespace detail { + +template +constexpr auto riemann_zeta_impl(T x) noexcept + BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) +{ + const auto fpc = fpclassify(x); + + constexpr T one { 1, 0 }; + + T result { }; + + if (fpc == FP_ZERO) + { + result = T { 5, -1, true }; + } + else if (fpc != FP_NORMAL) + { + if (fpc == FP_INFINITE) + { + result = (signbit(x) ? -std::numeric_limits::infinity() : one); + } + else + { + result = x; + } + } + else + { + // TODO(ckormanyos) Implement the Riemann-zeta function. + result = T { 0 }; + } + + return result; +} + +} //namespace detail + +BOOST_DECIMAL_EXPORT template +constexpr auto riemann_zeta(T x) noexcept + BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) +{ + #if BOOST_DECIMAL_DEC_EVAL_METHOD == 0 + + using evaluation_type = T; + + #elif BOOST_DECIMAL_DEC_EVAL_METHOD == 1 + + using evaluation_type = detail::promote_args_t; + + #else // BOOST_DECIMAL_DEC_EVAL_METHOD == 2 + + using evaluation_type = detail::promote_args_t; + + #endif + + return static_cast(detail::riemann_zeta_impl(static_cast(x))); +} + +} //namespace decimal +} //namespace boost + +#endif // BOOST_DECIMAL_DETAIL_CMATH_RIEMANN_ZETA_HPP diff --git a/test/Jamfile b/test/Jamfile index b066569a6..106c137c3 100644 --- a/test/Jamfile +++ b/test/Jamfile @@ -126,3 +126,4 @@ run test_tgamma.cpp ; run test_to_chars.cpp ; run test_to_string.cpp ; run test_type_traits.cpp ; +run test_zeta.cpp ; diff --git a/test/test_cbrt.cpp b/test/test_cbrt.cpp index 0cc07a32a..a0fdea428 100644 --- a/test/test_cbrt.cpp +++ b/test/test_cbrt.cpp @@ -3,12 +3,6 @@ // Distributed under the Boost Software License, Version 1.0. // https://www.boost.org/LICENSE_1_0.txt -#include -#include -#include -#include -#include - #include #if defined(__clang__) @@ -21,9 +15,10 @@ #include -#include - #include +#include +#include +#include #include namespace local diff --git a/test/test_zeta.cpp b/test/test_zeta.cpp new file mode 100644 index 000000000..17d47799e --- /dev/null +++ b/test/test_zeta.cpp @@ -0,0 +1,61 @@ +// Copyright 2024 Matt Borland +// Copyright 2024 Christoper Kormanyos +// Distributed under the Boost Software License, Version 1.0. +// https://www.boost.org/LICENSE_1_0.txt + +#include + +#if defined(__clang__) +# pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wfloat-equal" +#elif defined(__GNUC__) +# pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wfloat-equal" +#endif + +#include + +#include + +template +auto test_zeta() -> void +{ + using decimal_type = DecimalType; + using float_type = FloatType; + + // TODO(ckormanyos) Make actual tests whe the implementation is ready. + BOOST_TEST_EQ(riemann_zeta(decimal_type { 15 , -1 }), decimal_type { 0 }); + + { + std::mt19937_64 gen; + + std::uniform_real_distribution dist(static_cast(1.1L), static_cast(100.1L)); + + for(auto i = static_cast(UINT8_C(0)); i < static_cast(UINT8_C(10)); ++i) + { + static_cast(i); + + const decimal_type inf { std::numeric_limits::infinity() * static_cast(dist(gen)) }; + const decimal_type nan { std::numeric_limits::quiet_NaN() * static_cast(dist(gen)) }; + const decimal_type zero { decimal_type { 0 } * static_cast(dist(gen)) }; + + BOOST_TEST_EQ(riemann_zeta(inf), decimal_type { 1 } ); + BOOST_TEST (isinf(riemann_zeta(-inf)) && signbit(riemann_zeta(-inf))); + BOOST_TEST (isnan(riemann_zeta(nan))); + BOOST_TEST (isnan(riemann_zeta(-nan))); + + const decimal_type minus_half { -5, -1 }; + + BOOST_TEST_EQ(riemann_zeta(zero), minus_half); + } + } +} + +int main() +{ + using decimal_type = ::boost::decimal::decimal32; + + test_zeta(); + + return boost::report_errors(); +} From a5febf14084b536a63b142cd8002ca3a3c9b0b20 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sat, 25 May 2024 18:24:25 +0200 Subject: [PATCH 02/11] Preliminary Riemann-zeta imp/tests --- .../boost/decimal/detail/cmath/frexp10.hpp | 5 +- .../detail/cmath/impl/riemann_zeta_impl.hpp | 189 ++++++++++++++++++ .../decimal/detail/cmath/riemann_zeta.hpp | 147 +++++++++++++- test/test_zeta.cpp | 164 ++++++++++++++- 4 files changed, 487 insertions(+), 18 deletions(-) create mode 100644 include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp diff --git a/include/boost/decimal/detail/cmath/frexp10.hpp b/include/boost/decimal/detail/cmath/frexp10.hpp index 393494a59..85f7f4468 100644 --- a/include/boost/decimal/detail/cmath/frexp10.hpp +++ b/include/boost/decimal/detail/cmath/frexp10.hpp @@ -21,8 +21,9 @@ namespace boost { namespace decimal { // Returns the normalized significand and exponent to be cohort agnostic -// Returns num in the range [1'000'000, 9'999'999] -// +// Returns num in the range +// [1e06, 1e06 - 1] for decimal32 +// [1e15, 1e15 - 1] for decimal64 // If the conversion can not be performed returns UINT32_MAX and exp = 0 BOOST_DECIMAL_EXPORT template constexpr auto frexp10(T num, int* expptr) noexcept -> typename T::significand_type diff --git a/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp b/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp new file mode 100644 index 000000000..8f6028f1a --- /dev/null +++ b/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp @@ -0,0 +1,189 @@ +// 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_IMPL_RIEMANN_ZETA_IMPL_HPP +#define BOOST_DECIMAL_DETAIL_CMATH_IMPL_RIEMANN_ZETA_IMPL_HPP + +#include +#include + +#ifndef BOOST_DECIMAL_BUILD_MODULE +#include +#include +#include +#endif + +namespace boost { +namespace decimal { +namespace detail { + +namespace riemann_zeta_detail { + +template +struct prime_table_imp +{ + using decimal_type = T; + + using prime_table_t = std::array; + + // Table[Prime[n], {n, 1, 36, 1}] + static constexpr prime_table_t primes = + {{ + // Table[Prime[n], {n, 1, 26, 1}] + decimal_type { 2 }, decimal_type { 3 }, decimal_type { 5 }, decimal_type { 7 }, + decimal_type { 11 }, decimal_type { 13 }, decimal_type { 17 }, decimal_type { 19 }, + decimal_type { 23 }, decimal_type { 29 }, decimal_type { 31 }, decimal_type { 37 }, + decimal_type { 41 }, decimal_type { 43 }, decimal_type { 47 }, decimal_type { 53 }, + decimal_type { 59 }, decimal_type { 61 }, decimal_type { 67 }, decimal_type { 71 }, + decimal_type { 73 }, decimal_type { 79 }, decimal_type { 83 }, decimal_type { 89 }, + decimal_type { 97 }, decimal_type { 101 }, decimal_type { 103 }, decimal_type { 107 }, + decimal_type { 109 }, decimal_type { 113 }, decimal_type { 127 }, decimal_type { 131 }, + decimal_type { 137 }, decimal_type { 139 }, decimal_type { 149 }, decimal_type { 151 } + }}; +}; + +template +struct riemann_zeta_table_imp +{ +private: + using d32_coeffs_t = std::array; + using d32_fast_coeffs_t = std::array; + using d64_coeffs_t = std::array; + using d128_coeffs_t = std::array; + +public: + static constexpr d32_coeffs_t d32_coeffs = + {{ + // N[Series[Zeta[x], {x, 1, 4}], 19] + + +::boost::decimal::decimal32 { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma + +::boost::decimal::decimal32 { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) + -::boost::decimal::decimal32 { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 + -::boost::decimal::decimal32 { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 + +::boost::decimal::decimal32 { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 + }}; + + static constexpr d32_fast_coeffs_t d32_fast_coeffs = + {{ + // N[Series[Zeta[x], {x, 1, 4}], 19] + + +::boost::decimal::decimal32_fast { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma + +::boost::decimal::decimal32_fast { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) + -::boost::decimal::decimal32_fast { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 + -::boost::decimal::decimal32_fast { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 + +::boost::decimal::decimal32_fast { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 + }}; + + static constexpr d64_coeffs_t d64_coeffs = + {{ + // N[Series[Zeta[x], {x, 1, 6}], 19] + + +::boost::decimal::decimal64 { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma + +::boost::decimal::decimal64 { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) + -::boost::decimal::decimal64 { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 + -::boost::decimal::decimal64 { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 + +::boost::decimal::decimal64 { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 + -::boost::decimal::decimal64 { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 + -::boost::decimal::decimal64 { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 + }}; + + static constexpr d128_coeffs_t d128_coeffs = + {{ + +::boost::decimal::decimal128 { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma + +::boost::decimal::decimal128 { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) + -::boost::decimal::decimal128 { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 + -::boost::decimal::decimal128 { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 + +::boost::decimal::decimal128 { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 + -::boost::decimal::decimal128 { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 + -::boost::decimal::decimal128 { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 + }}; +}; + +#if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) + +template +constexpr typename riemann_zeta_table_imp::d32_coeffs_t riemann_zeta_table_imp::d32_coeffs; + +template +constexpr typename riemann_zeta_table_imp::d64_coeffs_t riemann_zeta_table_imp::d64_coeffs; + +template +constexpr typename riemann_zeta_table_imp::d128_coeffs_t riemann_zeta_table_imp::d128_coeffs; + +template +constexpr typename riemann_zeta_table_imp::d32_fast_coeffs_t riemann_zeta_table_imp::d32_fast_coeffs; + +template +constexpr typename prime_table_imp::prime_table_t prime_table_imp::primes; + +#endif + +} //namespace lgamma_detail + +using riemann_zeta_table = riemann_zeta_detail::riemann_zeta_table_imp; + +template +constexpr auto riemann_zeta_series_expansion(T x) noexcept; + +template <> +constexpr auto riemann_zeta_series_expansion(decimal32 x) noexcept +{ + return taylor_series_result(x, riemann_zeta_table::d32_coeffs); +} + +template <> +constexpr auto riemann_zeta_series_expansion(decimal32_fast x) noexcept +{ + return taylor_series_result(x, riemann_zeta_table::d32_fast_coeffs); +} + +template <> +constexpr auto riemann_zeta_series_expansion(decimal64 x) noexcept +{ + return taylor_series_result(x, riemann_zeta_table::d64_coeffs); +} + +template <> +constexpr auto riemann_zeta_series_expansion(decimal128 x) noexcept +{ + return taylor_series_result(x, riemann_zeta_table::d128_coeffs); +} + +template +using prime_table_t = riemann_zeta_detail::prime_table_imp::prime_table_t; + +template +using prime_table = riemann_zeta_detail::prime_table_imp; + +template +constexpr auto riemann_zeta_decimal_order(T x) noexcept + BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) +{ + int n { }; + + const T fr10 = frexp10(x, &n); + + constexpr int order_bias + { + std::numeric_limits::digits10 < 10 ? 6 + : std::numeric_limits::digits10 < 20 ? 15 + : 33 + }; + + return n + order_bias; +} + +template +constexpr auto riemann_zeta_factorial(int nf) noexcept + BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) +{ + return { nf <= 1 ? T { 1 } : riemann_zeta_factorial(nf - 1) * nf }; +} + +} //namespace detail +} //namespace decimal +} //namespace boost + +#endif //BOOST_DECIMAL_DETAIL_CMATH_IMPL_RIEMANN_ZETA_IMPL_HPP diff --git a/include/boost/decimal/detail/cmath/riemann_zeta.hpp b/include/boost/decimal/detail/cmath/riemann_zeta.hpp index 8ed02b891..8849ccd82 100644 --- a/include/boost/decimal/detail/cmath/riemann_zeta.hpp +++ b/include/boost/decimal/detail/cmath/riemann_zeta.hpp @@ -6,11 +6,15 @@ #ifndef BOOST_DECIMAL_DETAIL_CMATH_RIEMANN_ZETA_HPP #define BOOST_DECIMAL_DETAIL_CMATH_RIEMANN_ZETA_HPP -#include -#include -#include -#include +#include // NOLINT(llvm-include-order) +#include +#include +#include +#include +#include #include +#include +#include #ifndef BOOST_DECIMAL_BUILD_MODULE #include @@ -30,6 +34,8 @@ constexpr auto riemann_zeta_impl(T x) noexcept constexpr T one { 1, 0 }; + const bool is_neg { signbit(x) }; + T result { }; if (fpc == FP_ZERO) @@ -40,7 +46,7 @@ constexpr auto riemann_zeta_impl(T x) noexcept { if (fpc == FP_INFINITE) { - result = (signbit(x) ? -std::numeric_limits::infinity() : one); + result = (is_neg ? -std::numeric_limits::infinity() : one); } else { @@ -49,8 +55,135 @@ constexpr auto riemann_zeta_impl(T x) noexcept } else { - // TODO(ckormanyos) Implement the Riemann-zeta function. - result = T { 0 }; + if (is_neg) + { + // Handle Riemann-zeta reflection. + const T two_pi_term = pow(numbers::pi_v * 2, x) / numbers::pi_v; + const T chi = (two_pi_term * sin((numbers::pi_v * x) / 2)) * tgamma(one - x); + + result = chi * riemann_zeta(one - x); + } + else + { + constexpr int asymp_cutoff + { + std::numeric_limits::digits10 < 10 ? T { 2, 1 } // 20 + : std::numeric_limits::digits10 < 20 ? T { 5, 1 } // 50 + : T { 15, 1 } // 150 + }; + + if(x > asymp_cutoff) + { + result = one; + } + else if(x > T { 99, -2 } && x < T { 101, -2 }) + { + if(x > one || x < one) + { + // Use a Taylor series near the discontinuity at x=1. + + const T dx { x - one }; + + result = one / dx + detail::riemann_zeta_series_expansion(dx); + } + else + { + result = std::numeric_limits::quiet_NaN(); + } + } + else + { + // Test the conditions for the expansion of the product of primes. Set up a + // test for the product of primes. The expansion in the product of primes can + // be used if the number of prime-power terms remains reasonably small. + // This test for being reasonably "small" is checked in relation + // to the precision of the type and the size of primes available in the table. + + using prime_table_type = detail::prime_table_t; + + constexpr std::size_t n_primes = std::tuple_size::value; + + constexpr T lg10_max_prime { log10(detail::prime_table::primes[n_primes - 1U]) }; + + if((x * lg10_max_prime) > std::numeric_limits::digits10) + { + // Perform the product of primes. + result = one; + + for(std::size_t p = static_cast(UINT8_C(0)); p < n_primes; ++p) + { + const T prime_p_pow_s = pow(detail::prime_table::primes[p], x); + + const T pps_term { prime_p_pow_s / (prime_p_pow_s - one) }; + + if((pps_term - one) < std::numeric_limits::epsilon()) + { + break; + } + + result *= pps_term; + } + } + else + { + // Use the accelerated alternating converging series for Zeta as shown in: + // http://numbers.computation.free.fr/Constants/Miscellaneous/zetaevaluations.html + // taken from P. Borwein, "An Efficient Algorithm for the Riemann Zeta Function", + // January 1995. + + // Compute the coefficients dk in a loop and calculate the zeta function sum + // within the same loop on the fly. + + // Set up the factorials and powers for the calculation of the coefficients dk. + // Note that j = n at this stage in the calculation. Also note that the value of + // dn is equal to the value of d0 at the end of the loop. + + // Use N = (digits * 1.45) + {|imag(s)| * 1.1} + constexpr int nd { static_cast(std::numeric_limits::digits10 * 1.5F) }; + + bool neg_term = ((nd % 2) == 0); + + T n_plus_j_minus_one_fact = riemann_zeta_factorial((nd + nd) - 1); + T four_pow_j = pow(T { 4 }, nd); + T n_minus_j_fact = one; + T two_j_fact = n_plus_j_minus_one_fact * (2 * nd); + + T dn = (n_plus_j_minus_one_fact * four_pow_j) / (n_minus_j_fact * two_j_fact); + + T jps = pow(T { nd }, x); + + result = ((!neg_term) ? dn : -dn) / jps; + + for(auto j = nd - 1; j >= 0; --j) + { + const bool j_is_zero = (j == 0); + + const int two_jp1_two_j = ((2 * j) + 1) * (2 * ((!j_is_zero) ? j : 1)); + + n_plus_j_minus_one_fact /= (nd + j); + four_pow_j /= 4; + n_minus_j_fact *= (nd - j); + two_j_fact /= two_jp1_two_j; + + dn += ((n_plus_j_minus_one_fact * four_pow_j) / (n_minus_j_fact * two_j_fact)); + + if(!j_is_zero) + { + // Increment the zeta function sum. + jps = pow(T { j }, x); + + neg_term = (!neg_term); + + result += ((!neg_term) ? dn : -dn) / jps; + } + } + + const T two_pow_one_minus_s { pow(T { 2 }, one - x) }; + + result /= (dn * (one - two_pow_one_minus_s)); + } + } + } } return result; diff --git a/test/test_zeta.cpp b/test/test_zeta.cpp index 17d47799e..6da368177 100644 --- a/test/test_zeta.cpp +++ b/test/test_zeta.cpp @@ -3,8 +3,6 @@ // Distributed under the Boost Software License, Version 1.0. // https://www.boost.org/LICENSE_1_0.txt -#include - #if defined(__clang__) # pragma clang diagnostic push # pragma clang diagnostic ignored "-Wfloat-equal" @@ -13,23 +11,140 @@ # pragma GCC diagnostic ignored "-Wfloat-equal" #endif +#include + #include +#include +#include +#include #include +namespace local +{ + template + auto time_point() noexcept -> IntegralTimePointType + { + using local_integral_time_point_type = IntegralTimePointType; + using local_clock_type = ClockType; + + const auto current_now = + static_cast + ( + std::chrono::duration_cast + ( + local_clock_type::now().time_since_epoch() + ).count() + ); + + return static_cast(current_now); + } + + template + auto is_close_fraction(const NumericType& a, + const NumericType& b, + const NumericType& tol) noexcept -> bool + { + using std::fabs; + + auto result_is_ok = bool { }; + + NumericType delta { }; + + if(b == static_cast(0)) + { + delta = fabs(a - b); // LCOV_EXCL_LINE + + result_is_ok = (delta < tol); // LCOV_EXCL_LINE + } + else + { + delta = fabs(1 - (a / b)); + + result_is_ok = (delta < tol); + } + + // LCOV_EXCL_START + if (!result_is_ok) + { + std::cerr << std::setprecision(std::numeric_limits::digits10) << "a: " << a + << "\nb: " << b + << "\ndelta: " << delta + << "\ntol: " << tol << std::endl; + } + // LCOV_EXCL_STOP + + return result_is_ok; + } + + template + auto test_riemann_zeta(const int tol_factor) -> bool + { + using decimal_type = DecimalType; + using float_type = FloatType; + + std::random_device rd; + std::mt19937_64 gen(rd()); + + gen.seed(time_point()); + + auto dis_r = + std::uniform_real_distribution + { + static_cast(1.001L), + static_cast(12.3) + }; + + bool result_is_ok { true }; + + auto trials = static_cast(UINT8_C(0)); + + #if !defined(BOOST_DECIMAL_REDUCE_TEST_DEPTH) + constexpr auto count = (sizeof(decimal_type) == static_cast(UINT8_C(4))) ? static_cast(UINT32_C(0x80)) : static_cast(UINT32_C(0x20)); + #else + constexpr auto count = (sizeof(decimal_type) == static_cast(UINT8_C(4))) ? static_cast(UINT32_C(0x10)) : static_cast(UINT32_C(0x4)); + #endif + + for( ; trials < count; ++trials) + { + const auto x_flt = dis_r(gen); + const auto x_dec = static_cast(x_flt); + + const auto val_flt = boost::math::zeta(x_flt); + const auto val_dec = riemann_zeta(x_dec); + + const auto result_log_is_ok = is_close_fraction(val_flt, static_cast(val_dec), std::numeric_limits::epsilon() * static_cast(tol_factor)); + + result_is_ok = (result_log_is_ok && result_is_ok); + + if(!result_log_is_ok) + { + // LCOV_EXCL_START + std::cerr << "x_flt : " << std::scientific << std::setprecision(std::numeric_limits::digits10) << x_flt << std::endl; + std::cerr << "val_flt: " << std::scientific << std::setprecision(std::numeric_limits::digits10) << val_flt << std::endl; + std::cerr << "val_dec: " << std::scientific << std::setprecision(std::numeric_limits::digits10) << val_dec << std::endl; + + break; + // LCOV_EXCL_STOP + } + } + + BOOST_TEST(result_is_ok); + + return result_is_ok; + } + template -auto test_zeta() -> void +auto test_riemann_zeta_edge() -> void { using decimal_type = DecimalType; using float_type = FloatType; - // TODO(ckormanyos) Make actual tests whe the implementation is ready. - BOOST_TEST_EQ(riemann_zeta(decimal_type { 15 , -1 }), decimal_type { 0 }); - { std::mt19937_64 gen; - std::uniform_real_distribution dist(static_cast(1.1L), static_cast(100.1L)); + std::uniform_real_distribution dist(static_cast(1.1L), static_cast(101.1L)); for(auto i = static_cast(UINT8_C(0)); i < static_cast(UINT8_C(10)); ++i) { @@ -43,6 +158,7 @@ auto test_zeta() -> void BOOST_TEST (isinf(riemann_zeta(-inf)) && signbit(riemann_zeta(-inf))); BOOST_TEST (isnan(riemann_zeta(nan))); BOOST_TEST (isnan(riemann_zeta(-nan))); + BOOST_TEST (isnan(riemann_zeta(decimal_type { 1 }))); const decimal_type minus_half { -5, -1 }; @@ -51,11 +167,41 @@ auto test_zeta() -> void } } +} // namespace local + int main() { - using decimal_type = ::boost::decimal::decimal32; + bool result_is_ok { true }; + + { + using decimal_type = ::boost::decimal::decimal32; + + const bool result_rz32_is_ok = local::test_riemann_zeta(768); + + result_is_ok = (result_rz32_is_ok && result_is_ok); + + BOOST_TEST(result_is_ok); + } + + { + using decimal_type = ::boost::decimal::decimal64; + + const bool result_rz64_is_ok = local::test_riemann_zeta(2048); + + result_is_ok = (result_rz64_is_ok && result_is_ok); + + BOOST_TEST(result_is_ok); + } + + { + using decimal_type = ::boost::decimal::decimal32; + + local::test_riemann_zeta_edge(); + + BOOST_TEST(result_is_ok); + } - test_zeta(); + static_cast(result_is_ok); return boost::report_errors(); } From dfa3c71b30c6f69c1cac3fb761ec00283b85b95d Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sat, 25 May 2024 18:27:31 +0200 Subject: [PATCH 03/11] Correct a trivial typo on typename --- include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp b/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp index 8f6028f1a..6780a0784 100644 --- a/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp @@ -152,7 +152,7 @@ constexpr auto riemann_zeta_series_expansion(decimal128 x) noexcept } template -using prime_table_t = riemann_zeta_detail::prime_table_imp::prime_table_t; +using prime_table_t = typename riemann_zeta_detail::prime_table_imp::prime_table_t; template using prime_table = riemann_zeta_detail::prime_table_imp; From d88b8fd8c323e131bd266b7d996738017036d8ad Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sat, 25 May 2024 18:42:20 +0200 Subject: [PATCH 04/11] Silence a ton of Boost warnings --- test/test_zeta.cpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/test/test_zeta.cpp b/test/test_zeta.cpp index 6da368177..a3596b5c0 100644 --- a/test/test_zeta.cpp +++ b/test/test_zeta.cpp @@ -3,11 +3,24 @@ // Distributed under the Boost Software License, Version 1.0. // https://www.boost.org/LICENSE_1_0.txt +// Propogates up from boost.math +#define _SILENCE_CXX23_DENORM_DEPRECATION_WARNING + +#include + #if defined(__clang__) # pragma clang diagnostic push +# pragma clang diagnostic ignored "-Wold-style-cast" +# pragma clang diagnostic ignored "-Wundef" +# pragma clang diagnostic ignored "-Wconversion" +# pragma clang diagnostic ignored "-Wsign-conversion" # pragma clang diagnostic ignored "-Wfloat-equal" #elif defined(__GNUC__) # pragma GCC diagnostic push +# pragma GCC diagnostic ignored "-Wold-style-cast" +# pragma GCC diagnostic ignored "-Wundef" +# pragma GCC diagnostic ignored "-Wconversion" +# pragma GCC diagnostic ignored "-Wsign-conversion" # pragma GCC diagnostic ignored "-Wfloat-equal" #endif From db5cccc208e96d071ccb50411d59a1eea2031c0a Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sat, 25 May 2024 20:14:07 +0200 Subject: [PATCH 05/11] Up Riemann-zeta tols versus Boost.Math --- test/test_zeta.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_zeta.cpp b/test/test_zeta.cpp index a3596b5c0..d2bad880c 100644 --- a/test/test_zeta.cpp +++ b/test/test_zeta.cpp @@ -189,7 +189,7 @@ int main() { using decimal_type = ::boost::decimal::decimal32; - const bool result_rz32_is_ok = local::test_riemann_zeta(768); + const bool result_rz32_is_ok = local::test_riemann_zeta(2048); result_is_ok = (result_rz32_is_ok && result_is_ok); @@ -199,7 +199,7 @@ int main() { using decimal_type = ::boost::decimal::decimal64; - const bool result_rz64_is_ok = local::test_riemann_zeta(2048); + const bool result_rz64_is_ok = local::test_riemann_zeta(8192); result_is_ok = (result_rz64_is_ok && result_is_ok); From 87f30ca84bb84004bd8381aa5cefc8968d7764d9 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sat, 25 May 2024 21:14:32 +0200 Subject: [PATCH 06/11] Reduce a test range --- test/test_zeta.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_zeta.cpp b/test/test_zeta.cpp index d2bad880c..c89447e95 100644 --- a/test/test_zeta.cpp +++ b/test/test_zeta.cpp @@ -105,7 +105,7 @@ namespace local auto dis_r = std::uniform_real_distribution { - static_cast(1.001L), + static_cast(1.005L), static_cast(12.3) }; From 014f0295117fd8624f646b645262adab4ca591e5 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 26 May 2024 09:03:22 +0200 Subject: [PATCH 07/11] Finish Riemann-zeta --- .../detail/cmath/impl/riemann_zeta_impl.hpp | 67 +++--- .../decimal/detail/cmath/riemann_zeta.hpp | 33 +-- test/test_zeta.cpp | 194 +++++++++++++++--- 3 files changed, 232 insertions(+), 62 deletions(-) diff --git a/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp b/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp index 6780a0784..dd4b9a451 100644 --- a/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp @@ -28,10 +28,10 @@ struct prime_table_imp using prime_table_t = std::array; - // Table[Prime[n], {n, 1, 36, 1}] static constexpr prime_table_t primes = {{ - // Table[Prime[n], {n, 1, 26, 1}] + // Table[Prime[n], {n, 1, 36, 1}] + decimal_type { 2 }, decimal_type { 3 }, decimal_type { 5 }, decimal_type { 7 }, decimal_type { 11 }, decimal_type { 13 }, decimal_type { 17 }, decimal_type { 19 }, decimal_type { 23 }, decimal_type { 29 }, decimal_type { 31 }, decimal_type { 37 }, @@ -48,56 +48,73 @@ template struct riemann_zeta_table_imp { private: - using d32_coeffs_t = std::array; - using d32_fast_coeffs_t = std::array; - using d64_coeffs_t = std::array; - using d128_coeffs_t = std::array; + using d32_coeffs_t = std::array; + using d32_fast_coeffs_t = std::array; + using d64_coeffs_t = std::array; + using d128_coeffs_t = std::array; public: static constexpr d32_coeffs_t d32_coeffs = {{ - // N[Series[Zeta[x], {x, 1, 4}], 19] + // N[Series[Zeta[x], {x, 1, 6}], 19] +::boost::decimal::decimal32 { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma +::boost::decimal::decimal32 { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) -::boost::decimal::decimal32 { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 -::boost::decimal::decimal32 { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 +::boost::decimal::decimal32 { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 + -::boost::decimal::decimal32 { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 + -::boost::decimal::decimal32 { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 }}; static constexpr d32_fast_coeffs_t d32_fast_coeffs = {{ - // N[Series[Zeta[x], {x, 1, 4}], 19] + // N[Series[Zeta[x], {x, 1, 6}], 19] +::boost::decimal::decimal32_fast { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma +::boost::decimal::decimal32_fast { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) -::boost::decimal::decimal32_fast { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 -::boost::decimal::decimal32_fast { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 +::boost::decimal::decimal32_fast { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 + -::boost::decimal::decimal32_fast { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 + -::boost::decimal::decimal32_fast { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 }}; static constexpr d64_coeffs_t d64_coeffs = {{ - // N[Series[Zeta[x], {x, 1, 6}], 19] - - +::boost::decimal::decimal64 { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma - +::boost::decimal::decimal64 { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) - -::boost::decimal::decimal64 { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 - -::boost::decimal::decimal64 { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 - +::boost::decimal::decimal64 { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 - -::boost::decimal::decimal64 { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 - -::boost::decimal::decimal64 { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 + // N[Series[Zeta[x], {x, 1, 9}], 19] + + +::boost::decimal::decimal64 { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma + +::boost::decimal::decimal64 { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) + -::boost::decimal::decimal64 { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 + -::boost::decimal::decimal64 { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 + +::boost::decimal::decimal64 { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 + -::boost::decimal::decimal64 { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 + -::boost::decimal::decimal64 { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 + +::boost::decimal::decimal64 { UINT64_C(1046209458447918742), - 19 - 6 }, // * (x - 1)^7 + -::boost::decimal::decimal64 { UINT64_C(8733218100273797361), - 19 - 8 }, // * (x - 1)^8 + +::boost::decimal::decimal64 { UINT64_C(9478277782762358956), - 19 - 10 }, // * (x - 1)^9 }}; static constexpr d128_coeffs_t d128_coeffs = {{ - +::boost::decimal::decimal128 { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma - +::boost::decimal::decimal128 { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) - -::boost::decimal::decimal128 { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 - -::boost::decimal::decimal128 { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 - +::boost::decimal::decimal128 { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 - -::boost::decimal::decimal128 { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 - -::boost::decimal::decimal128 { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 + // N[Series[Zeta[x], {x, 1, 14}], 36] + + +::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(312909238939453), UINT64_C(7916302232898517972) }, -34 }, + +::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(394735489323855), UINT64_C(10282954930524890450) }, -35 }, + -::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(262657820647143), UINT64_C(7801536535536173172) }, -36 }, + -::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(185564311701532), UINT64_C(15687007158497646588) }, -37 }, + +::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(525244016002584), UINT64_C(12277750447068982866) }, -38 }, + -::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(358384752584293), UINT64_C(18370286456371002882) }, -39 }, + -::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(179773779887752), UINT64_C(17772011513518515048) }, -40 }, + +::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(56715128386205), UINT64_C(15292499466693711883) }, -40 }, + -::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(473428701855329), UINT64_C(926484760170384186) }, -42 }, + +::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(513818468174601), UINT64_C(18105240268308765734) }, -44 }, + +::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(306743667337648), UINT64_C(15567754919026551912) }, -44 }, + -::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(366931412745108), UINT64_C(2220247416524400302) }, -45 }, + +::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(189307984255553), UINT64_C(8448217616480074192) }, -46 }, + +::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(239089604329878), UINT64_C(14831803080673374292) }, -48 }, + -::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(130092671757244), UINT64_C(16458215134170057406) }, -48 }, }}; }; @@ -179,7 +196,7 @@ template constexpr auto riemann_zeta_factorial(int nf) noexcept BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T) { - return { nf <= 1 ? T { 1 } : riemann_zeta_factorial(nf - 1) * nf }; + return { (nf <= 1) ? T { 1 } : riemann_zeta_factorial(nf - 1) * nf }; } } //namespace detail diff --git a/include/boost/decimal/detail/cmath/riemann_zeta.hpp b/include/boost/decimal/detail/cmath/riemann_zeta.hpp index 8849ccd82..14b1b8318 100644 --- a/include/boost/decimal/detail/cmath/riemann_zeta.hpp +++ b/include/boost/decimal/detail/cmath/riemann_zeta.hpp @@ -58,10 +58,12 @@ constexpr auto riemann_zeta_impl(T x) noexcept if (is_neg) { // Handle Riemann-zeta reflection. - const T two_pi_term = pow(numbers::pi_v * 2, x) / numbers::pi_v; - const T chi = (two_pi_term * sin((numbers::pi_v * x) / 2)) * tgamma(one - x); - result = chi * riemann_zeta(one - x); + const T dmx { one - x }; + const T two_pi_term { pow(numbers::pi_v * 2, x) / numbers::pi_v }; + const T chi { (two_pi_term * sin((numbers::pi_v * x) / 2)) * tgamma(dmx) }; + + result = chi * riemann_zeta(dmx); } else { @@ -74,20 +76,24 @@ constexpr auto riemann_zeta_impl(T x) noexcept if(x > asymp_cutoff) { + // Check the argument is large and then simply return 1. + result = one; } - else if(x > T { 99, -2 } && x < T { 101, -2 }) + else if((x > T { 99, -2 }) && (x < T { 101, -2 })) { - if(x > one || x < one) + if((x > one) || (x < one)) { // Use a Taylor series near the discontinuity at x=1. const T dx { x - one }; - result = one / dx + detail::riemann_zeta_series_expansion(dx); + result = (one / dx) + detail::riemann_zeta_series_expansion(dx); } else { + // The argument is exaclty one. The result is complex-infinity. + result = std::numeric_limits::quiet_NaN(); } } @@ -103,7 +109,7 @@ constexpr auto riemann_zeta_impl(T x) noexcept constexpr std::size_t n_primes = std::tuple_size::value; - constexpr T lg10_max_prime { log10(detail::prime_table::primes[n_primes - 1U]) }; + const T lg10_max_prime { log10(detail::prime_table::primes[n_primes - 1U]) }; if((x * lg10_max_prime) > std::numeric_limits::digits10) { @@ -138,12 +144,15 @@ constexpr auto riemann_zeta_impl(T x) noexcept // Note that j = n at this stage in the calculation. Also note that the value of // dn is equal to the value of d0 at the end of the loop. - // Use N = (digits * 1.45) + {|imag(s)| * 1.1} + // Use nd = (digits * 1.45) + {|imag(s)| * 1.1} + // Here we have, however, only real valued arguments so this + // is streamlined a tad. Also instead of 1.45, we simply use 1.5. + constexpr int nd { static_cast(std::numeric_limits::digits10 * 1.5F) }; - bool neg_term = ((nd % 2) == 0); + bool neg_term { (nd % 2) == 0 }; - T n_plus_j_minus_one_fact = riemann_zeta_factorial((nd + nd) - 1); + T n_plus_j_minus_one_fact = riemann_zeta_factorial(2 * nd - 1); T four_pow_j = pow(T { 4 }, nd); T n_minus_j_fact = one; T two_j_fact = n_plus_j_minus_one_fact * (2 * nd); @@ -156,9 +165,9 @@ constexpr auto riemann_zeta_impl(T x) noexcept for(auto j = nd - 1; j >= 0; --j) { - const bool j_is_zero = (j == 0); + const bool j_is_zero { j == 0 }; - const int two_jp1_two_j = ((2 * j) + 1) * (2 * ((!j_is_zero) ? j : 1)); + const int two_jp1_two_j { ((2 * j) + 1) * (2 * ((!j_is_zero) ? j : 1)) }; n_plus_j_minus_one_fact /= (nd + j); four_pow_j /= 4; diff --git a/test/test_zeta.cpp b/test/test_zeta.cpp index c89447e95..1a226e274 100644 --- a/test/test_zeta.cpp +++ b/test/test_zeta.cpp @@ -92,7 +92,7 @@ namespace local } template - auto test_riemann_zeta(const int tol_factor) -> bool + auto test_riemann_zeta(const int tol_factor, const long double range_lo, const long double range_hi) -> bool { using decimal_type = DecimalType; using float_type = FloatType; @@ -102,11 +102,11 @@ namespace local gen.seed(time_point()); - auto dis_r = + auto dis = std::uniform_real_distribution { - static_cast(1.005L), - static_cast(12.3) + static_cast(range_lo), + static_cast(range_hi) }; bool result_is_ok { true }; @@ -121,7 +121,7 @@ namespace local for( ; trials < count; ++trials) { - const auto x_flt = dis_r(gen); + const auto x_flt = dis(gen); const auto x_dec = static_cast(x_flt); const auto val_flt = boost::math::zeta(x_flt); @@ -180,41 +180,185 @@ auto test_riemann_zeta_edge() -> void } } +auto test_riemann_zeta_128_lo(const int tol_factor) -> bool +{ + using decimal_type = boost::decimal::decimal128; + + using str_ctrl_array_type = std::array; + + const str_ctrl_array_type ctrl_strings = + {{ + // Table[N[Zeta[1 + n/1000], 36], {n, 5, 7, 1}] + "200.577579622956683652084654605524346", + "167.244319052140751595350397994287573", + "143.434867995431699170218293588670480", + }}; + + std::array::value> rz_values { }; + std::array::value> ctrl_values { }; + + int nx { 5 }; + + bool result_is_ok { true }; + + const decimal_type my_tol { std::numeric_limits::epsilon() * static_cast(tol_factor) }; + + for(auto i = static_cast(UINT8_C(0)); i < std::tuple_size::value; ++i) + { + const decimal_type x_arg = + decimal_type + { + decimal_type { 1 } + + decimal_type { nx, -3 } + }; + + ++nx; + + rz_values[i] = riemann_zeta(x_arg); + + static_cast + ( + from_chars(ctrl_strings[i], ctrl_strings[i] + std::strlen(ctrl_strings[i]), ctrl_values[i]) + ); + + const auto result_rz_is_ok = is_close_fraction(rz_values[i], ctrl_values[i], my_tol); + + result_is_ok = (result_rz_is_ok && result_is_ok); + } + + return result_is_ok; +} + +auto test_riemann_zeta_128_hi(const int tol_factor) -> bool +{ + using decimal_type = boost::decimal::decimal128; + + using str_ctrl_array_type = std::array; + + const str_ctrl_array_type ctrl_strings = + {{ + // Table[N[Zeta[n + n/10], 36], {n, 1, 9, 1}] + "10.5844484649508098263864007917355230", + "1.49054325650689350825344649551165452", + "1.15194479472077368855082683374115056", + "1.05928172597983541766404502818685201", + "1.02520457995468569459240582819540529", + "1.01116101415427096427312532266653516", + "1.00504987929596499812178165124883599", + "1.00231277790982194674469422849347780", + "1.00106679698357801585766465214764188", + }}; + + std::array::value> rz_values { }; + std::array::value> ctrl_values { }; + + int nx { 1 }; + + bool result_is_ok { true }; + + const decimal_type my_tol { std::numeric_limits::epsilon() * static_cast(tol_factor) }; + + for(auto i = static_cast(UINT8_C(0)); i < std::tuple_size::value; ++i) + { + const decimal_type x_arg = + decimal_type + { + decimal_type { nx } + + decimal_type { nx, -1 } + }; + + ++nx; + + rz_values[i] = riemann_zeta(x_arg); + + static_cast + ( + from_chars(ctrl_strings[i], ctrl_strings[i] + std::strlen(ctrl_strings[i]), ctrl_values[i]) + ); + + const auto result_rz_is_ok = is_close_fraction(rz_values[i], ctrl_values[i], my_tol); + + result_is_ok = (result_rz_is_ok && result_is_ok); + } + + return result_is_ok; +} + } // namespace local int main() { - bool result_is_ok { true }; + bool result_is_ok { true }; - { - using decimal_type = ::boost::decimal::decimal32; + { + using decimal_type = ::boost::decimal::decimal32; - const bool result_rz32_is_ok = local::test_riemann_zeta(2048); + const bool result_rz32_is_ok = local::test_riemann_zeta(128, 1.1L, 5.6L); - result_is_ok = (result_rz32_is_ok && result_is_ok); + result_is_ok = (result_rz32_is_ok && result_is_ok); - BOOST_TEST(result_is_ok); - } + BOOST_TEST(result_is_ok); + } - { - using decimal_type = ::boost::decimal::decimal64; + { + using decimal_type = ::boost::decimal::decimal32; - const bool result_rz64_is_ok = local::test_riemann_zeta(8192); + const bool result_rz32_is_ok = local::test_riemann_zeta(1024, 1.005L, 1.1L); - result_is_ok = (result_rz64_is_ok && result_is_ok); + result_is_ok = (result_rz32_is_ok && result_is_ok); - BOOST_TEST(result_is_ok); - } + BOOST_TEST(result_is_ok); + } - { - using decimal_type = ::boost::decimal::decimal32; + { + using decimal_type = ::boost::decimal::decimal32; - local::test_riemann_zeta_edge(); + const bool result_rz32_is_ok = local::test_riemann_zeta(512, -3.6L, -2.3L); - BOOST_TEST(result_is_ok); - } + result_is_ok = (result_rz32_is_ok && result_is_ok); + + BOOST_TEST(result_is_ok); + } + + { + using decimal_type = ::boost::decimal::decimal64; + + const bool result_rz64_is_ok = local::test_riemann_zeta(256, 1.1L, 12.3L); + + result_is_ok = (result_rz64_is_ok && result_is_ok); + + BOOST_TEST(result_is_ok); + } + + { + using decimal_type = ::boost::decimal::decimal64; + + const bool result_rz64_is_ok = local::test_riemann_zeta(1024, 1.005L, 1.1L); + + result_is_ok = (result_rz64_is_ok && result_is_ok); + + BOOST_TEST(result_is_ok); + } + + { + using decimal_type = ::boost::decimal::decimal32; + + local::test_riemann_zeta_edge(); + + BOOST_TEST(result_is_ok); + } + + { + const auto result_rz_128_lo_is_ok = local::test_riemann_zeta_128_lo(4096); + const auto result_rz_128_hi_is_ok = local::test_riemann_zeta_128_hi(4096); + + BOOST_TEST(result_rz_128_lo_is_ok); + BOOST_TEST(result_rz_128_hi_is_ok); + + result_is_ok = (result_rz_128_lo_is_ok && result_rz_128_hi_is_ok && result_is_ok); + } - static_cast(result_is_ok); + result_is_ok = ((boost::report_errors() == 0) && result_is_ok); - return boost::report_errors(); + return (result_is_ok ? 0 : -1); } From 528a6e46cad5ca7bb295390c6a5937d37dbf1622 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 26 May 2024 10:01:08 +0200 Subject: [PATCH 08/11] Increase depth of Riemann-zeta tests --- test/test_zeta.cpp | 67 +++++++++++++++++++++++++++++++++------------- 1 file changed, 48 insertions(+), 19 deletions(-) diff --git a/test/test_zeta.cpp b/test/test_zeta.cpp index 1a226e274..8918fccc0 100644 --- a/test/test_zeta.cpp +++ b/test/test_zeta.cpp @@ -33,6 +33,10 @@ #include #include +template auto my_zero() -> DecimalType&; +template auto my_nan () -> DecimalType&; +template auto my_inf () -> DecimalType&; + namespace local { template -auto test_riemann_zeta_edge() -> void +auto test_riemann_zeta_edge() -> bool { using decimal_type = DecimalType; using float_type = FloatType; + bool result_is_ok { true }; + { std::mt19937_64 gen; + gen.seed(time_point()); + std::uniform_real_distribution dist(static_cast(1.1L), static_cast(101.1L)); for(auto i = static_cast(UINT8_C(0)); i < static_cast(UINT8_C(10)); ++i) { static_cast(i); - const decimal_type inf { std::numeric_limits::infinity() * static_cast(dist(gen)) }; - const decimal_type nan { std::numeric_limits::quiet_NaN() * static_cast(dist(gen)) }; - const decimal_type zero { decimal_type { 0 } * static_cast(dist(gen)) }; + { + const decimal_type inf { ::my_inf () * decimal_type(dist(gen)) }; + + bool result_inf_is_ok { }; - BOOST_TEST_EQ(riemann_zeta(inf), decimal_type { 1 } ); - BOOST_TEST (isinf(riemann_zeta(-inf)) && signbit(riemann_zeta(-inf))); - BOOST_TEST (isnan(riemann_zeta(nan))); - BOOST_TEST (isnan(riemann_zeta(-nan))); - BOOST_TEST (isnan(riemann_zeta(decimal_type { 1 }))); + result_inf_is_ok = (riemann_zeta(inf) == decimal_type { 1 } ); result_is_ok = (result_inf_is_ok && result_is_ok); BOOST_TEST(result_is_ok); + result_inf_is_ok = ( isinf(riemann_zeta(-inf)) + && signbit(riemann_zeta(-inf))); result_is_ok = (result_inf_is_ok && result_is_ok); BOOST_TEST(result_is_ok); + } - const decimal_type minus_half { -5, -1 }; + { + const decimal_type nan { ::my_nan () * decimal_type(dist(gen)) }; - BOOST_TEST_EQ(riemann_zeta(zero), minus_half); + bool result_nan_is_ok { }; + + result_nan_is_ok = (isnan(riemann_zeta(nan))); result_is_ok = (result_nan_is_ok && result_is_ok); BOOST_TEST(result_is_ok); + result_nan_is_ok = (isnan(riemann_zeta(-nan))); result_is_ok = (result_nan_is_ok && result_is_ok); BOOST_TEST(result_is_ok); + result_nan_is_ok = (isnan(riemann_zeta(decimal_type { 1 }))); result_is_ok = (result_nan_is_ok && result_is_ok); BOOST_TEST(result_is_ok); + } + + { + const decimal_type zero { ::my_zero() * decimal_type(dist(gen)) }; + + const decimal_type minus_half { -5, -1 }; + + const bool result_zero_is_ok = (riemann_zeta(zero) == minus_half); result_is_ok = (result_zero_is_ok && result_is_ok); BOOST_TEST(result_is_ok); + } } } + + return result_is_ok; } auto test_riemann_zeta_128_lo(const int tol_factor) -> bool @@ -290,6 +314,15 @@ int main() { bool result_is_ok { true }; + { + using decimal_type = ::boost::decimal::decimal32; + + const bool result_edge_is_ok = local::test_riemann_zeta_edge(); + + result_is_ok = (result_edge_is_ok && result_is_ok); + + BOOST_TEST(result_is_ok); + } { using decimal_type = ::boost::decimal::decimal32; @@ -340,14 +373,6 @@ int main() BOOST_TEST(result_is_ok); } - { - using decimal_type = ::boost::decimal::decimal32; - - local::test_riemann_zeta_edge(); - - BOOST_TEST(result_is_ok); - } - { const auto result_rz_128_lo_is_ok = local::test_riemann_zeta_128_lo(4096); const auto result_rz_128_hi_is_ok = local::test_riemann_zeta_128_hi(4096); @@ -362,3 +387,7 @@ int main() return (result_is_ok ? 0 : -1); } + +template auto my_zero() -> DecimalType& { using decimal_type = DecimalType; static decimal_type val_zero { 0 }; return val_zero; } +template auto my_nan () -> DecimalType& { using decimal_type = DecimalType; static decimal_type val_nan { std::numeric_limits::quiet_NaN() }; return val_nan; } +template auto my_inf () -> DecimalType& { using decimal_type = DecimalType; static decimal_type val_inf { std::numeric_limits::infinity() }; return val_inf; } From d94bdb1808772e1b84a8511b7eb8af332ef15e7d Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 26 May 2024 12:57:27 +0200 Subject: [PATCH 09/11] Use Pade approx for decimal32/fast32 --- .../detail/cmath/impl/riemann_zeta_impl.hpp | 86 +++++++++---------- .../decimal/detail/cmath/riemann_zeta.hpp | 6 +- test/test_zeta.cpp | 15 +++- 3 files changed, 56 insertions(+), 51 deletions(-) diff --git a/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp b/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp index dd4b9a451..8decfbf4f 100644 --- a/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp @@ -48,38 +48,10 @@ template struct riemann_zeta_table_imp { private: - using d32_coeffs_t = std::array; - using d32_fast_coeffs_t = std::array; - using d64_coeffs_t = std::array; - using d128_coeffs_t = std::array; + using d64_coeffs_t = std::array; + using d128_coeffs_t = std::array; public: - static constexpr d32_coeffs_t d32_coeffs = - {{ - // N[Series[Zeta[x], {x, 1, 6}], 19] - - +::boost::decimal::decimal32 { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma - +::boost::decimal::decimal32 { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) - -::boost::decimal::decimal32 { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 - -::boost::decimal::decimal32 { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 - +::boost::decimal::decimal32 { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 - -::boost::decimal::decimal32 { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 - -::boost::decimal::decimal32 { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 - }}; - - static constexpr d32_fast_coeffs_t d32_fast_coeffs = - {{ - // N[Series[Zeta[x], {x, 1, 6}], 19] - - +::boost::decimal::decimal32_fast { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma - +::boost::decimal::decimal32_fast { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) - -::boost::decimal::decimal32_fast { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 - -::boost::decimal::decimal32_fast { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 - +::boost::decimal::decimal32_fast { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 - -::boost::decimal::decimal32_fast { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 - -::boost::decimal::decimal32_fast { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 - }}; - static constexpr d64_coeffs_t d64_coeffs = {{ // N[Series[Zeta[x], {x, 1, 9}], 19] @@ -120,18 +92,12 @@ struct riemann_zeta_table_imp #if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) -template -constexpr typename riemann_zeta_table_imp::d32_coeffs_t riemann_zeta_table_imp::d32_coeffs; - template constexpr typename riemann_zeta_table_imp::d64_coeffs_t riemann_zeta_table_imp::d64_coeffs; template constexpr typename riemann_zeta_table_imp::d128_coeffs_t riemann_zeta_table_imp::d128_coeffs; -template -constexpr typename riemann_zeta_table_imp::d32_fast_coeffs_t riemann_zeta_table_imp::d32_fast_coeffs; - template constexpr typename prime_table_imp::prime_table_t prime_table_imp::primes; @@ -142,30 +108,60 @@ constexpr typename prime_table_imp::prime_table_t prime_table_imp::p using riemann_zeta_table = riemann_zeta_detail::riemann_zeta_table_imp; template -constexpr auto riemann_zeta_series_expansion(T x) noexcept; +constexpr auto riemann_zeta_series_or_pade_expansion(T x) noexcept; template <> -constexpr auto riemann_zeta_series_expansion(decimal32 x) noexcept +constexpr auto riemann_zeta_series_or_pade_expansion(decimal32 x) noexcept { - return taylor_series_result(x, riemann_zeta_table::d32_coeffs); + const decimal32 top = + (decimal32 { UINT64_C(7025346442393055904), -19 + 1 } + + x * (decimal32 { UINT64_C(6331631438687936980), -19 + 1 } + + x * decimal32 { UINT64_C(1671529107642800378), -19 + 1 })); + + const decimal32 bot = + (decimal32 { UINT64_C(1402850698872379326), -19 + 2, true } + + x * (decimal32 { UINT64_C(1302850698872379326), -19 + 2 } + + x * decimal32 { 1 })); + + return top / bot; } template <> -constexpr auto riemann_zeta_series_expansion(decimal32_fast x) noexcept +constexpr auto riemann_zeta_series_or_pade_expansion(decimal32_fast x) noexcept { - return taylor_series_result(x, riemann_zeta_table::d32_fast_coeffs); + const decimal32_fast top = + (decimal32_fast { UINT64_C(7025346442393055904), -19 + 1 } + + x * (decimal32_fast { UINT64_C(6331631438687936980), -19 + 1 } + + x * decimal32_fast { UINT64_C(1671529107642800378), -19 + 1 })); + + const decimal32_fast bot = + (decimal32_fast { UINT64_C(1402850698872379326), -19 + 2, true } + + x * (decimal32_fast { UINT64_C(1302850698872379326), -19 + 2 } + + x * decimal32_fast { 1 })); + + return top / bot; } template <> -constexpr auto riemann_zeta_series_expansion(decimal64 x) noexcept +constexpr auto riemann_zeta_series_or_pade_expansion(decimal64 x) noexcept { - return taylor_series_result(x, riemann_zeta_table::d64_coeffs); + // TODO(ckormanyos) Consider using a Pade approximation for 64-bit. + + constexpr decimal64 one { 1 }; + + const decimal64 dx { x - one }; + + return one / dx + taylor_series_result(dx, riemann_zeta_table::d64_coeffs); } template <> -constexpr auto riemann_zeta_series_expansion(decimal128 x) noexcept +constexpr auto riemann_zeta_series_or_pade_expansion(decimal128 x) noexcept { - return taylor_series_result(x, riemann_zeta_table::d128_coeffs); + constexpr decimal128 one { 1 }; + + const decimal128 dx { x - one }; + + return one / dx + taylor_series_result(dx, riemann_zeta_table::d128_coeffs); } template diff --git a/include/boost/decimal/detail/cmath/riemann_zeta.hpp b/include/boost/decimal/detail/cmath/riemann_zeta.hpp index 14b1b8318..5bf2646ee 100644 --- a/include/boost/decimal/detail/cmath/riemann_zeta.hpp +++ b/include/boost/decimal/detail/cmath/riemann_zeta.hpp @@ -84,11 +84,9 @@ constexpr auto riemann_zeta_impl(T x) noexcept { if((x > one) || (x < one)) { - // Use a Taylor series near the discontinuity at x=1. + // Use a Taylor series (or Pade approximation) near the discontinuity at x=1. - const T dx { x - one }; - - result = (one / dx) + detail::riemann_zeta_series_expansion(dx); + result = detail::riemann_zeta_series_or_pade_expansion(x); } else { diff --git a/test/test_zeta.cpp b/test/test_zeta.cpp index 8918fccc0..88975039b 100644 --- a/test/test_zeta.cpp +++ b/test/test_zeta.cpp @@ -323,6 +323,7 @@ int main() BOOST_TEST(result_is_ok); } + { using decimal_type = ::boost::decimal::decimal32; @@ -336,7 +337,17 @@ int main() { using decimal_type = ::boost::decimal::decimal32; - const bool result_rz32_is_ok = local::test_riemann_zeta(1024, 1.005L, 1.1L); + const bool result_rz32_is_ok = local::test_riemann_zeta(512, 1.05L, 1.1L); + + result_is_ok = (result_rz32_is_ok && result_is_ok); + + BOOST_TEST(result_is_ok); + } + + { + using decimal_type = ::boost::decimal::decimal32_fast; + + const bool result_rz32_is_ok = local::test_riemann_zeta(512, 1.05L, 1.1L); result_is_ok = (result_rz32_is_ok && result_is_ok); @@ -366,7 +377,7 @@ int main() { using decimal_type = ::boost::decimal::decimal64; - const bool result_rz64_is_ok = local::test_riemann_zeta(1024, 1.005L, 1.1L); + const bool result_rz64_is_ok = local::test_riemann_zeta(4096, 1.005L, 1.1L); result_is_ok = (result_rz64_is_ok && result_is_ok); From b01a4967aa63dd50e3cc14920c26531578c37308 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 26 May 2024 15:36:50 +0200 Subject: [PATCH 10/11] Tune up Riemann-zeta and its tests --- .../detail/cmath/impl/riemann_zeta_impl.hpp | 125 ++++++++++++++---- .../decimal/detail/cmath/riemann_zeta.hpp | 2 +- test/test_zeta.cpp | 6 +- 3 files changed, 106 insertions(+), 27 deletions(-) diff --git a/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp b/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp index 8decfbf4f..8c23c1749 100644 --- a/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp +++ b/include/boost/decimal/detail/cmath/impl/riemann_zeta_impl.hpp @@ -48,10 +48,38 @@ template struct riemann_zeta_table_imp { private: - using d64_coeffs_t = std::array; - using d128_coeffs_t = std::array; + using d32_coeffs_t = std::array; + using d32_fast_coeffs_t = std::array; + using d64_coeffs_t = std::array; + using d128_coeffs_t = std::array; public: + static constexpr d32_coeffs_t d32_coeffs = + {{ + // N[Series[Zeta[x], {x, 1, 6}], 19] + + +::boost::decimal::decimal32 { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma + +::boost::decimal::decimal32 { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) + -::boost::decimal::decimal32 { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 + -::boost::decimal::decimal32 { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 + +::boost::decimal::decimal32 { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 + -::boost::decimal::decimal32 { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 + -::boost::decimal::decimal32 { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 + }}; + + static constexpr d32_fast_coeffs_t d32_fast_coeffs = + {{ + // N[Series[Zeta[x], {x, 1, 6}], 19] + + +::boost::decimal::decimal32_fast { UINT64_C(5772156649015328606), - 19 - 0 }, // EulerGamma + +::boost::decimal::decimal32_fast { UINT64_C(7281584548367672486), - 19 - 1 }, // * (x - 1) + -::boost::decimal::decimal32_fast { UINT64_C(4845181596436159242), - 19 - 2 }, // * (x - 1)^2 + -::boost::decimal::decimal32_fast { UINT64_C(3423057367172243110), - 19 - 3 }, // * (x - 1)^3 + +::boost::decimal::decimal32_fast { UINT64_C(9689041939447083573), - 19 - 4 }, // * (x - 1)^4 + -::boost::decimal::decimal32_fast { UINT64_C(6611031810842189181), - 19 - 5 }, // * (x - 1)^5 + -::boost::decimal::decimal32_fast { UINT64_C(3316240908752772359), - 19 - 6 }, // * (x - 1)^6 + }}; + static constexpr d64_coeffs_t d64_coeffs = {{ // N[Series[Zeta[x], {x, 1, 9}], 19] @@ -92,6 +120,12 @@ struct riemann_zeta_table_imp #if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900) +template +constexpr typename riemann_zeta_table_imp::d32_coeffs_t riemann_zeta_table_imp::d32_coeffs; + +template +constexpr typename riemann_zeta_table_imp::d32_fast_coeffs_t riemann_zeta_table_imp::d32_fast_coeffs; + template constexpr typename riemann_zeta_table_imp::d64_coeffs_t riemann_zeta_table_imp::d64_coeffs; @@ -113,45 +147,90 @@ constexpr auto riemann_zeta_series_or_pade_expansion(T x) noexcept; template <> constexpr auto riemann_zeta_series_or_pade_expansion(decimal32 x) noexcept { - const decimal32 top = - (decimal32 { UINT64_C(7025346442393055904), -19 + 1 } - + x * (decimal32 { UINT64_C(6331631438687936980), -19 + 1 } - + x * decimal32 { UINT64_C(1671529107642800378), -19 + 1 })); + constexpr decimal32 one { 1 }; - const decimal32 bot = - (decimal32 { UINT64_C(1402850698872379326), -19 + 2, true } - + x * (decimal32 { UINT64_C(1302850698872379326), -19 + 2 } - + x * decimal32 { 1 })); + const decimal32 dx { x - one }; - return top / bot; + if (fabs(dx) < decimal32 { 5, - 2 }) + { + return one / dx + taylor_series_result(dx, riemann_zeta_table::d32_coeffs); + } + else + { + const decimal32 top = + (decimal32 { UINT64_C(7025346442393055904), -19 + 1 } + + x * (decimal32 { UINT64_C(6331631438687936980), -19 + 1 } + + x * decimal32 { UINT64_C(1671529107642800378), -19 + 1 })); + + const decimal32 bot = + (decimal32 { UINT64_C(1402850698872379326), -19 + 2, true } + + x * (decimal32 { UINT64_C(1302850698872379326), -19 + 2 } + + x * decimal32 { 1 })); + + return top / bot; + } } template <> constexpr auto riemann_zeta_series_or_pade_expansion(decimal32_fast x) noexcept { - const decimal32_fast top = - (decimal32_fast { UINT64_C(7025346442393055904), -19 + 1 } - + x * (decimal32_fast { UINT64_C(6331631438687936980), -19 + 1 } - + x * decimal32_fast { UINT64_C(1671529107642800378), -19 + 1 })); + constexpr decimal32_fast one { 1 }; - const decimal32_fast bot = - (decimal32_fast { UINT64_C(1402850698872379326), -19 + 2, true } - + x * (decimal32_fast { UINT64_C(1302850698872379326), -19 + 2 } - + x * decimal32_fast { 1 })); + const decimal32_fast dx { x - one }; - return top / bot; + if (fabs(dx) < decimal32_fast { 5, -2 }) + { + return one / dx + taylor_series_result(dx, riemann_zeta_table::d32_fast_coeffs); + } + else + { + const decimal32_fast top = + (decimal32_fast { UINT64_C(7025346442393055904), -19 + 1 } + + x * (decimal32_fast { UINT64_C(6331631438687936980), -19 + 1 } + + x * decimal32_fast { UINT64_C(1671529107642800378), -19 + 1 })); + + const decimal32_fast bot = + (decimal32_fast { UINT64_C(1402850698872379326), -19 + 2, true } + + x * (decimal32_fast { UINT64_C(1302850698872379326), -19 + 2 } + + x * decimal32_fast { 1 })); + + return top / bot; + } } template <> constexpr auto riemann_zeta_series_or_pade_expansion(decimal64 x) noexcept { - // TODO(ckormanyos) Consider using a Pade approximation for 64-bit. - constexpr decimal64 one { 1 }; const decimal64 dx { x - one }; - return one / dx + taylor_series_result(dx, riemann_zeta_table::d64_coeffs); + if (fabs(dx) < decimal64 { 5, -2 }) + { + return one / dx + taylor_series_result(dx, riemann_zeta_table::d64_coeffs); + } + else + { + constexpr decimal64 c0 { UINT64_C(4124764818173475125), - 19 + 5 }; + constexpr decimal64 c1 { UINT64_C(4582078064035558510), - 19 + 5 }; + constexpr decimal64 c2 { UINT64_C(1806662427082674333), - 19 + 5 }; + constexpr decimal64 c3 { UINT64_C(3281232347201801441), - 19 + 4 }; + constexpr decimal64 c4 { UINT64_C(3092253262304078300), - 19 + 3 }; + constexpr decimal64 c5 { UINT64_C(1985384224421766402), - 19 + 2 }; + constexpr decimal64 c6 { UINT64_C(1016070109033501213), - 19 + 1 }; + + constexpr decimal64 d0 { UINT64_C(8249529636338921254), - 19 + 5, true }; + constexpr decimal64 d1 { UINT64_C(5997465199121809585), - 19 + 5 }; + constexpr decimal64 d2 { UINT64_C(1915568444415559307), - 19 + 5 }; + constexpr decimal64 d3 { UINT64_C(3021354370625514285), - 19 + 4 }; + constexpr decimal64 d4 { UINT64_C(3227310996533313801), - 19 + 3 }; + constexpr decimal64 d5 { UINT64_C(1987445773667795184), - 19 + 2 }; + + const decimal64 top { c0 + x * (c1 + x * (c2 + x * (c3 + x * (c4 + x * (c5+ x * c6))))) }; + const decimal64 bot { d0 + x * (d1 + x * (d2 + x * (d3 + x * (d4 + x * (d5 + x))))) }; + + return top / bot; + } } template <> diff --git a/include/boost/decimal/detail/cmath/riemann_zeta.hpp b/include/boost/decimal/detail/cmath/riemann_zeta.hpp index 5bf2646ee..91ef98f06 100644 --- a/include/boost/decimal/detail/cmath/riemann_zeta.hpp +++ b/include/boost/decimal/detail/cmath/riemann_zeta.hpp @@ -80,7 +80,7 @@ constexpr auto riemann_zeta_impl(T x) noexcept result = one; } - else if((x > T { 99, -2 }) && (x < T { 101, -2 })) + else if((x > T { 9, -1 }) && (x < T { 11, -1 })) { if((x > one) || (x < one)) { diff --git a/test/test_zeta.cpp b/test/test_zeta.cpp index 88975039b..db2c63158 100644 --- a/test/test_zeta.cpp +++ b/test/test_zeta.cpp @@ -337,7 +337,7 @@ int main() { using decimal_type = ::boost::decimal::decimal32; - const bool result_rz32_is_ok = local::test_riemann_zeta(512, 1.05L, 1.1L); + const bool result_rz32_is_ok = local::test_riemann_zeta(1024, 1.01L, 1.1L); result_is_ok = (result_rz32_is_ok && result_is_ok); @@ -347,7 +347,7 @@ int main() { using decimal_type = ::boost::decimal::decimal32_fast; - const bool result_rz32_is_ok = local::test_riemann_zeta(512, 1.05L, 1.1L); + const bool result_rz32_is_ok = local::test_riemann_zeta(1024, 1.01L, 1.1L); result_is_ok = (result_rz32_is_ok && result_is_ok); @@ -377,7 +377,7 @@ int main() { using decimal_type = ::boost::decimal::decimal64; - const bool result_rz64_is_ok = local::test_riemann_zeta(4096, 1.005L, 1.1L); + const bool result_rz64_is_ok = local::test_riemann_zeta(1024, 1.01L, 1.1L); result_is_ok = (result_rz64_is_ok && result_is_ok); From df60ff6cd10164c5d57cc426862a4a5896024997 Mon Sep 17 00:00:00 2001 From: Christopher Kormanyos Date: Sun, 26 May 2024 16:56:33 +0200 Subject: [PATCH 11/11] Hit some unrelated missed cover lines --- .../decimal/detail/cmath/riemann_zeta.hpp | 9 +--- include/boost/decimal/detail/cmath/tgamma.hpp | 10 +++-- test/test_tgamma.cpp | 44 ++++++++++++++----- 3 files changed, 40 insertions(+), 23 deletions(-) diff --git a/include/boost/decimal/detail/cmath/riemann_zeta.hpp b/include/boost/decimal/detail/cmath/riemann_zeta.hpp index 91ef98f06..2db228965 100644 --- a/include/boost/decimal/detail/cmath/riemann_zeta.hpp +++ b/include/boost/decimal/detail/cmath/riemann_zeta.hpp @@ -44,14 +44,7 @@ constexpr auto riemann_zeta_impl(T x) noexcept } else if (fpc != FP_NORMAL) { - if (fpc == FP_INFINITE) - { - result = (is_neg ? -std::numeric_limits::infinity() : one); - } - else - { - result = x; - } + result = ((fpc == FP_INFINITE) ? (is_neg ? -std::numeric_limits::infinity() : one) : x); } else { diff --git a/include/boost/decimal/detail/cmath/tgamma.hpp b/include/boost/decimal/detail/cmath/tgamma.hpp index 7fa9d29b4..e1cca95de 100644 --- a/include/boost/decimal/detail/cmath/tgamma.hpp +++ b/include/boost/decimal/detail/cmath/tgamma.hpp @@ -33,31 +33,33 @@ constexpr auto tgamma_impl(T x) noexcept const auto is_pure_int = (nx == x); + const bool is_neg = signbit(x); + const auto fpc = fpclassify(x); if (fpc != FP_NORMAL) { if (fpc == FP_ZERO) { - result = (signbit(x) ? -std::numeric_limits::infinity() : std::numeric_limits::infinity()); + result = (is_neg ? -std::numeric_limits::infinity() : std::numeric_limits::infinity()); } else if(fpc == FP_INFINITE) { - result = (signbit(x) ? std::numeric_limits::quiet_NaN() : std::numeric_limits::infinity()); + result = (is_neg ? std::numeric_limits::quiet_NaN() : std::numeric_limits::infinity()); } else { result = x; } } - else if ((nx < 0) && is_pure_int && ((nx & 1) != 0)) + else if (is_neg && is_pure_int) { // Pure negative integer argument. result = std::numeric_limits::quiet_NaN(); } else { - if (signbit(x)) + if (is_neg) { // Reflection for negative argument. const auto ga = tgamma(-x); diff --git a/test/test_tgamma.cpp b/test/test_tgamma.cpp index 6cb0edd68..4601063ba 100644 --- a/test/test_tgamma.cpp +++ b/test/test_tgamma.cpp @@ -21,8 +21,8 @@ #include -template auto my_zero() -> DecimalType& { using decimal_type = DecimalType; static decimal_type val_zero { 0, 0 }; return val_zero; } -template auto my_one () -> DecimalType& { using decimal_type = DecimalType; static decimal_type val_one { 1, 0 }; return val_one; } +template auto my_zero() -> DecimalType&; +template auto my_one () -> DecimalType&; namespace local { @@ -304,7 +304,7 @@ namespace local { static_cast(i); - const auto val_zero_neg = tgamma(-::my_zero()); + const auto val_zero_neg = tgamma(-(::my_zero() * static_cast(dist(gen)))); const auto result_val_zero_neg_is_ok = (isinf(val_zero_neg) && signbit(val_zero_neg)); @@ -315,13 +315,13 @@ namespace local for(auto i = static_cast(UINT8_C(0)); i < static_cast(UINT8_C(6)); ++i) { - const auto n_neg = static_cast(-static_cast(i) - 1); + decimal_type dnx { ::my_zero() * static_cast(dist(gen)) }; - const auto val_neg_int = tgamma( decimal_type { n_neg, 0 } ); + dnx -= static_cast(i + 1); - const auto is_odd = ((n_neg & 1) != 0); + const auto val_neg_int = tgamma(dnx); - const auto result_val_neg_int_is_ok = (is_odd ? isnan(val_neg_int) : (fpclassify(val_neg_int) == FP_NORMAL)); + const auto result_val_neg_int_is_ok = isnan(val_neg_int); BOOST_TEST(result_val_neg_int_is_ok); @@ -444,11 +444,11 @@ namespace local { using decimal_type = boost::decimal::decimal128; - using str_ctrl_array_type = std::array; + using str_ctrl_array_type = std::array; const str_ctrl_array_type ctrl_strings = {{ - // Table[N[Gamma[n + n/10 + n/100 + n/1000], 36], {n, 1, 131, 10}] + // Table[N[Gamma[n + n/10 + n/100 + n/1000], 36], {n, 1, 191, 10}] "0.947008281162266001895790481785841941", "6.86303089001025022525468906807854872E7", "3.15793281780505944512262743601561476E21", @@ -462,7 +462,13 @@ namespace local "4.76763037027821868276349648015607359E180", "4.62395515046183569847627307320617350E203", "1.22680267570425015175034111397510637E227", - "8.18925182002285090986591692926519438E250" + "8.18925182002285090986591692926519438E250", + "1.28134405415265103961333749220602490E275", + "4.42253704896092684478952587741331734E299", + "3.19451354412535995695298989136255493E324", + "4.61171076932972412633999770033353340E349", + "1.27758574231803927960543278875893523E375", + "6.55079490827236721494047351435261992E400", }}; std::array::value> tg_values { }; @@ -530,6 +536,17 @@ auto main() -> int result_is_ok = (result_tgamma_is_ok && result_is_ok); } + { + using decimal_type = boost::decimal::decimal32_fast; + using float_type = float; + + const auto result_tgamma_is_ok = local::test_tgamma(768, 2.1L, 23.4L); + + BOOST_TEST(result_tgamma_is_ok); + + result_is_ok = (result_tgamma_is_ok && result_is_ok); + } + { using decimal_type = boost::decimal::decimal64; using float_type = double; @@ -589,7 +606,9 @@ auto main() -> int { const auto result_tgamma128_lo_is_ok = local::test_tgamma_128_lo(4096); - const auto result_tgamma128_hi_is_ok = local::test_tgamma_128_hi(4096); + + // TODO(ckormanyos) Can we get better asymptotic accuracy with more coefficients at 128-bit? + const auto result_tgamma128_hi_is_ok = local::test_tgamma_128_hi(0x10000); BOOST_TEST(result_tgamma128_lo_is_ok); BOOST_TEST(result_tgamma128_hi_is_ok); @@ -601,3 +620,6 @@ auto main() -> int return (result_is_ok ? 0 : -1); } + +template auto my_zero() -> DecimalType& { using decimal_type = DecimalType; static decimal_type val_zero { 0 }; return val_zero; } +template auto my_one () -> DecimalType& { using decimal_type = DecimalType; static decimal_type val_one { 1 }; return val_one; }