mirror of
https://github.com/boostorg/multiprecision.git
synced 2026-01-19 04:22:11 +00:00
410 lines
14 KiB
C++
410 lines
14 KiB
C++
///////////////////////////////////////////////////////////////////
|
|
// Copyright Christopher Kormanyos 2020 - 2022. //
|
|
// Distributed under the Boost Software License, //
|
|
// Version 1.0. (See accompanying file LICENSE_1_0.txt //
|
|
// or copy at http://www.boost.org/LICENSE_1_0.txt) //
|
|
///////////////////////////////////////////////////////////////////
|
|
|
|
#include <array>
|
|
#include <cstdint>
|
|
#include <ctime>
|
|
#include <iomanip>
|
|
#include <iostream>
|
|
#include <memory>
|
|
#include <utility>
|
|
#include <vector>
|
|
|
|
#define EXAMPLE008_BERNOULLI_USE_LOCAL_PI
|
|
|
|
#if !defined(BOOST_MP_STANDALONE)
|
|
#define BOOST_MP_STANDALONE
|
|
#endif
|
|
|
|
#if !defined(EXAMPLE008_BERNOULLI_USE_LOCAL_PI)
|
|
#if !defined(BOOST_MATH_STANDALONE)
|
|
#define BOOST_MATH_STANDALONE
|
|
#endif
|
|
#include <boost/math/constants/constants.hpp>
|
|
#endif
|
|
|
|
#include <boost/multiprecision/cpp_bin_float.hpp>
|
|
|
|
// cd /mnt/c/Users/User/Documents/Ks/PC_Software/Test
|
|
// g++ -Wall -Wextra -Wconversion -Wsign-conversion -Wshadow -Wundef -O3 -std=c++11 -I/mnt/c/boost/modular_boost/boost/libs/multiprecision/include -I/mnt/c/boost/modular_boost/boost/libs/config/include standalone_bernoulli_tgamma.cpp -o standalone_bernoulli_tgamma.exe
|
|
|
|
// D:\MinGW_nuwen\MinGW\bin\g++ -Wall -Wextra -Wconversion -Wsign-conversion -Wshadow -Wundef -O3 -std=c++11 -IC:\boost\modular_boost\boost\libs\multiprecision\include -IC:\boost\modular_boost\boost\libs\config\include test.cpp -o test.exe
|
|
namespace example008_bernoulli
|
|
{
|
|
constexpr std::int32_t wide_decimal_digits10 = INT32_C(1001);
|
|
|
|
using wide_float_backend_type = boost::multiprecision::cpp_bin_float<wide_decimal_digits10, boost::multiprecision::digit_base_10, std::allocator<void>>;
|
|
|
|
using wide_float_type = boost::multiprecision::number<wide_float_backend_type, boost::multiprecision::et_off>;
|
|
|
|
template<typename FloatingPointType>
|
|
auto pi() -> FloatingPointType
|
|
{
|
|
return static_cast<FloatingPointType>(3.1415926535897932384626433832795029L); // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
|
|
}
|
|
|
|
#if defined(EXAMPLE008_BERNOULLI_USE_LOCAL_PI)
|
|
template<typename FloatingPointType>
|
|
auto calc_pi() -> FloatingPointType
|
|
{
|
|
// Use a quadratically convergent Gauss AGM to compute pi.
|
|
|
|
using floating_point_type = FloatingPointType;
|
|
|
|
floating_point_type val_pi;
|
|
|
|
floating_point_type a(1U);
|
|
|
|
// Initialize bB to 0.5.
|
|
floating_point_type bB(0.5F); // NOLINT(readability-identifier-naming,cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
|
|
|
|
// Initialize t to 0.375.
|
|
floating_point_type t(static_cast<floating_point_type>(3U) / 8U); // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
|
|
|
|
floating_point_type s(bB);
|
|
|
|
// This loop is designed for a maximum of several million
|
|
// decimal digits of pi. The index k should reach no higher
|
|
// than about 25 or 30. After about 20 iterations, the precision
|
|
// is about one million decimal digits.
|
|
|
|
const auto digits10_iteration_goal =
|
|
static_cast<std::uint32_t>
|
|
(
|
|
static_cast<std::uint32_t>(std::numeric_limits<floating_point_type>::digits10 / 2)
|
|
+ static_cast<std::uint32_t>(9U)
|
|
);
|
|
|
|
using std::log;
|
|
using std::lround;
|
|
|
|
const auto digits10_scale =
|
|
static_cast<std::uint32_t>
|
|
(
|
|
lround
|
|
(
|
|
static_cast<float>(1000.0F * log(static_cast<float>(std::numeric_limits<floating_point_type>::radix))) / log(10.0F)
|
|
)
|
|
);
|
|
|
|
for(auto k = static_cast<unsigned>(UINT8_C(0));
|
|
k < static_cast<unsigned>(UINT8_C(48));
|
|
++k)
|
|
{
|
|
using std::sqrt;
|
|
|
|
a += sqrt(bB);
|
|
a /= 2U;
|
|
val_pi = a;
|
|
val_pi *= a;
|
|
bB = val_pi;
|
|
bB -= t;
|
|
bB *= 2U;
|
|
|
|
floating_point_type iterate_term(bB);
|
|
|
|
iterate_term -= val_pi;
|
|
iterate_term *= static_cast<unsigned long long>(1ULL << (k + 1U)); // NOLINT(google-runtime-int)
|
|
|
|
s += iterate_term;
|
|
|
|
// Test the number of precise digits from this iteration.
|
|
// If it is there are enough precise digits, then the calculation
|
|
// is finished.
|
|
const auto ib =
|
|
(std::max)
|
|
(
|
|
static_cast<std::int32_t>(0),
|
|
static_cast<std::int32_t>(-ilogb(iterate_term))
|
|
);
|
|
|
|
const auto digits10_of_iteration =
|
|
static_cast<std::uint32_t>
|
|
(
|
|
static_cast<std::uint64_t>(static_cast<std::uint64_t>(ib) * digits10_scale) / UINT32_C(1000)
|
|
);
|
|
|
|
// Estimate the approximate decimal digits of this iteration term.
|
|
// If we have attained about half of the total desired digits
|
|
// with this iteration term, then the calculation is finished
|
|
// because the change from the next iteration will be
|
|
// insignificantly small.
|
|
|
|
if(digits10_of_iteration > digits10_iteration_goal)
|
|
{
|
|
break;
|
|
}
|
|
|
|
t = val_pi;
|
|
t += bB;
|
|
t /= 4U;
|
|
}
|
|
|
|
return (val_pi + bB) / s;
|
|
}
|
|
|
|
template<typename FloatingPointType>
|
|
auto my_pi() -> const FloatingPointType&
|
|
{
|
|
using floating_point_type = FloatingPointType;
|
|
|
|
static const floating_point_type local_pi = calc_pi<floating_point_type>();
|
|
|
|
return local_pi;
|
|
}
|
|
#endif
|
|
|
|
template<>
|
|
auto pi() -> wide_float_type
|
|
{
|
|
#if defined(EXAMPLE008_BERNOULLI_USE_LOCAL_PI)
|
|
return my_pi<wide_float_type>();
|
|
#else
|
|
return boost::math::constants::pi<wide_float_type>();
|
|
#endif
|
|
}
|
|
|
|
auto bernoulli_table() -> std::vector<wide_float_type>&
|
|
{
|
|
static std::vector<wide_float_type>
|
|
bernoulli_table
|
|
(
|
|
static_cast<std::size_t>
|
|
(
|
|
static_cast<float>(std::numeric_limits<wide_float_type>::digits10) * 0.95F
|
|
)
|
|
);
|
|
|
|
return bernoulli_table;
|
|
}
|
|
|
|
template<typename FloatingPointType>
|
|
auto bernoulli_b(FloatingPointType* bn, std::uint32_t n) -> void
|
|
{
|
|
using floating_point_type = FloatingPointType;
|
|
|
|
// See reference "Computing Bernoulli and Tangent Numbers", Richard P. Brent.
|
|
// See also the book Richard P. Brent and Paul Zimmermann, "Modern Computer Arithmetic",
|
|
// Cambridge University Press, 2010, p. 237.
|
|
|
|
const auto m = static_cast<std::uint32_t>(n / 2U);
|
|
|
|
std::vector<floating_point_type> tangent_numbers(m + 1U);
|
|
|
|
tangent_numbers[0U] = 0U;
|
|
tangent_numbers[1U] = 1U;
|
|
|
|
for(std::uint32_t k = 1U; k < m; ++k)
|
|
{
|
|
tangent_numbers[k + 1U] = tangent_numbers[k] * k;
|
|
}
|
|
|
|
for(auto k = static_cast<std::uint32_t>(2U); k <= m; ++k)
|
|
{
|
|
for(auto j = k; j <= m; ++j)
|
|
{
|
|
const std::uint32_t j_minus_k = j - k;
|
|
|
|
tangent_numbers[j] = (tangent_numbers[j - 1] * j_minus_k)
|
|
+ (tangent_numbers[j] * (j_minus_k + 2U));
|
|
}
|
|
}
|
|
|
|
floating_point_type two_pow_two_m(4U);
|
|
|
|
for(std::uint32_t i = 1U; i < static_cast<std::uint32_t>(n / 2U); ++i)
|
|
{
|
|
const auto two_i = static_cast<std::uint32_t>(2U * i);
|
|
|
|
const floating_point_type b = (tangent_numbers[i] * two_i) / (two_pow_two_m * (two_pow_two_m - 1));
|
|
|
|
const bool b_neg = ((two_i % 4U) == 0U);
|
|
|
|
bn[two_i] = ((!b_neg) ? b : -b); // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
|
|
|
|
two_pow_two_m *= 4U;
|
|
}
|
|
|
|
bn[0U] = 1U; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
|
|
bn[1U] = floating_point_type(-1) / 2U; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
|
|
}
|
|
|
|
template<typename FloatingPointType>
|
|
auto tgamma(const FloatingPointType& x) -> FloatingPointType
|
|
{
|
|
using floating_point_type = FloatingPointType;
|
|
|
|
// Check if the argument should be scaled up for the Bernoulli series expansion.
|
|
static const auto min_arg_n =
|
|
static_cast<std::int32_t>
|
|
(
|
|
static_cast<float>(static_cast<float>(std::numeric_limits<floating_point_type>::digits10) * 0.8F)
|
|
);
|
|
|
|
static const floating_point_type min_arg_x = floating_point_type(min_arg_n);
|
|
|
|
const auto n_recur =
|
|
static_cast<std::uint32_t>
|
|
(
|
|
(x < min_arg_x) ? static_cast<std::uint32_t>(static_cast<std::uint32_t>(min_arg_n - static_cast<std::int32_t>(x)) + 1U)
|
|
: static_cast<std::uint32_t>(0U)
|
|
);
|
|
|
|
floating_point_type xx(x);
|
|
|
|
// Scale the argument up and use downward recursion later for the final result.
|
|
if(n_recur != 0U)
|
|
{
|
|
xx += n_recur;
|
|
}
|
|
|
|
floating_point_type one_over_x_pow_two_n_minus_one = 1 / xx;
|
|
const floating_point_type one_over_x2 = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one;
|
|
floating_point_type sum = (one_over_x_pow_two_n_minus_one * bernoulli_table()[2U]) / 2U;
|
|
|
|
floating_point_type tol = std::numeric_limits<floating_point_type>::epsilon();
|
|
|
|
using std::log;
|
|
|
|
if(xx > 8U) // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
|
|
{
|
|
// In the following code sequence, we extract the approximate logarithm
|
|
// of the argument x and use the leading term of Stirling's approximation,
|
|
// which is Log[Gamma[x]] aprox. (x (Log[x] - 1)) in order to scale
|
|
// the tolerance. In order to do this, we find the built-in floating point
|
|
// approximation of (x (Log[x] - 1)).
|
|
|
|
// Limit fx to the range 8 <= fx <= 10^16, where 8 is chosen to
|
|
// ensure that (log(fx) - 1.0F) remains positive and 10^16 is
|
|
// selected arbitrarily, yet ensured to be rather large.
|
|
auto fx_max = (std::max)(static_cast<floating_point_type>(8U), xx); // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
|
|
|
|
auto fx = (std::min)(fx_max, static_cast<floating_point_type>(UINT64_C(10000000000000000)));
|
|
|
|
tol *= static_cast<float>(fx * (log(fx) - 1.0F));
|
|
}
|
|
|
|
// Perform the Bernoulli series expansion.
|
|
for(auto n2 = static_cast<std::uint32_t>(4U); n2 < static_cast<std::uint32_t>(bernoulli_table().size()); n2 += 2U)
|
|
{
|
|
one_over_x_pow_two_n_minus_one *= one_over_x2;
|
|
|
|
const floating_point_type term =
|
|
(one_over_x_pow_two_n_minus_one * bernoulli_table()[n2])
|
|
/ static_cast<std::uint64_t>(static_cast<std::uint64_t>(n2) * static_cast<std::uint32_t>(n2 - 1U));
|
|
|
|
using std::fabs;
|
|
|
|
if((n2 > 6U) && (fabs(term) < tol)) // NOLINT(cppcoreguidelines-avoid-magic-numbers,readability-magic-numbers)
|
|
{
|
|
break;
|
|
}
|
|
|
|
sum += term;
|
|
}
|
|
|
|
using example008_bernoulli::pi;
|
|
using std::exp;
|
|
|
|
static const floating_point_type half = floating_point_type(1U) / 2U;
|
|
static const floating_point_type half_ln_two_pi = log(pi<floating_point_type>() * 2U) / 2U;
|
|
|
|
floating_point_type g = exp(((((xx - half) * log(xx)) - xx) + half_ln_two_pi) + sum);
|
|
|
|
// Rescale the result using downward recursion if necessary.
|
|
for(auto k = static_cast<std::uint32_t>(0U); k < n_recur; ++k)
|
|
{
|
|
g /= --xx;
|
|
}
|
|
|
|
return g;
|
|
}
|
|
} // namespace example008_bernoulli
|
|
|
|
auto example008_bernoulli_tgamma() -> bool
|
|
{
|
|
const std::clock_t start = std::clock();
|
|
|
|
// Initialize the table of Bernoulli numbers.
|
|
example008_bernoulli::bernoulli_b
|
|
(
|
|
example008_bernoulli::bernoulli_table().data(),
|
|
static_cast<std::uint32_t>(example008_bernoulli::bernoulli_table().size())
|
|
);
|
|
|
|
// In this example, we compute values of Gamma[1/2 + n].
|
|
|
|
// We will make use of the relation
|
|
// (2n)!
|
|
// Gamma[1/2 + n] = -------- * Sqrt[Pi].
|
|
// (4^n) n!
|
|
|
|
// Table[Factorial[2 n]/((4^n) Factorial[n]), {n, 0, 17, 1}]
|
|
const std::array<std::pair<std::uint64_t, std::uint32_t>, 18U> ratios =
|
|
{{
|
|
{ UINT64_C( 1), UINT32_C( 1) },
|
|
{ UINT64_C( 1), UINT32_C( 2) },
|
|
{ UINT64_C( 3), UINT32_C( 4) },
|
|
{ UINT64_C( 15), UINT32_C( 8) },
|
|
{ UINT64_C( 105), UINT32_C( 16) },
|
|
{ UINT64_C( 945), UINT32_C( 32) },
|
|
{ UINT64_C( 10395), UINT32_C( 64) },
|
|
{ UINT64_C( 135135), UINT32_C( 128) },
|
|
{ UINT64_C( 2027025), UINT32_C( 256) },
|
|
{ UINT64_C( 34459425), UINT32_C( 512) },
|
|
{ UINT64_C( 654729075), UINT32_C( 1024) },
|
|
{ UINT64_C( 13749310575), UINT32_C( 2048) },
|
|
{ UINT64_C( 316234143225), UINT32_C( 4096) },
|
|
{ UINT64_C( 7905853580625), UINT32_C( 8192) },
|
|
{ UINT64_C( 213458046676875), UINT32_C( 16384) },
|
|
{ UINT64_C( 6190283353629375), UINT32_C( 32768) },
|
|
{ UINT64_C( 191898783962510625), UINT32_C( 65536) },
|
|
{ UINT64_C(6332659870762850625), UINT32_C(131072) }
|
|
}};
|
|
|
|
bool result_is_ok = true;
|
|
|
|
using example008_bernoulli::wide_float_type;
|
|
|
|
const wide_float_type tol (std::numeric_limits<wide_float_type>::epsilon() * UINT32_C(100000));
|
|
const wide_float_type half(0.5F);
|
|
|
|
for(auto i = static_cast<std::size_t>(0U); i < ratios.size(); ++i)
|
|
{
|
|
// Calculate Gamma[1/2 + i]
|
|
|
|
const wide_float_type g = example008_bernoulli::tgamma(half + i);
|
|
|
|
// Calculate the control value.
|
|
|
|
using example008_bernoulli::pi;
|
|
using std::fabs;
|
|
using std::sqrt;
|
|
|
|
const wide_float_type control = (sqrt(pi<wide_float_type>()) * ratios[i].first) / ratios[i].second; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index)
|
|
|
|
const wide_float_type closeness = fabs(1 - (g / control));
|
|
|
|
result_is_ok &= (closeness < tol);
|
|
}
|
|
|
|
const std::clock_t stop = std::clock();
|
|
|
|
std::cout << "Time example008_bernoulli_tgamma(): "
|
|
<< static_cast<float>(stop - start) / static_cast<float>(CLOCKS_PER_SEC)
|
|
<< std::endl;
|
|
|
|
return result_is_ok;
|
|
}
|
|
|
|
int main()
|
|
{
|
|
const bool result_is_ok = example008_bernoulli_tgamma();
|
|
|
|
std::cout << "result_is_ok: " << std::boolalpha << result_is_ok << std::endl;
|
|
}
|