boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Merge pull request #634 from boostorg/bernoulli_threading_2

Browse Source 
Update for Bernoulli and Lanczos code.
This commit is contained in:
jzmaddock
2021-05-25 08:26:08 +01:00
committed by GitHub
parent 755b81992e 84b8e4ac81
commit 90278fcf9e
5 changed files with 1016 additions and 1314 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
include/boost/math/special_functions/detail/bernoulli_details.hpp
View File

@@ -565,11 +565,12 @@ private:
};
template <class T, class Policy>
inline bernoulli_numbers_cache<T, Policy>& get_bernoulli_numbers_cache()
inline typename std::enable_if<(std::numeric_limits<T>::digits == 0) || (std::numeric_limits<T>::digits >= INT_MAX), bernoulli_numbers_cache<T, Policy>&>::type get_bernoulli_numbers_cache()
{
//
// Force this function to be called at program startup so all the static variables
// get initialized then (thread safety).
// When numeric_limits<>::digits is zero, the type has either not specialized numeric_limits at all
// or it's precision can vary at runtime. So make the cache thread_local so that each thread can
// have it's own precision if required:
//
static
#ifndef BOOST_MATH_NO_THREAD_LOCAL_WITH_NON_TRIVIAL_TYPES
@@ -578,6 +579,15 @@ inline bernoulli_numbers_cache<T, Policy>& get_bernoulli_numbers_cache()
bernoulli_numbers_cache<T, Policy> data;
return data;
}
template <class T, class Policy>
inline typename std::enable_if<std::numeric_limits<T>::digits && (std::numeric_limits<T>::digits < INT_MAX), bernoulli_numbers_cache<T, Policy>&>::type get_bernoulli_numbers_cache()
{
//
// Note that we rely on C++11 thread-safe initialization here:
//
static bernoulli_numbers_cache<T, Policy> data;
return data;
}
}}}

10
include/boost/math/special_functions/detail/lgamma_small.hpp
View File

@@ -482,7 +482,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant<int, 113>&, c
return result;
}
template <class T, class Policy, class Lanczos>
T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant<int, 0>&, const Policy& pol, const Lanczos&)
T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant<int, 0>&, const Policy& pol, const Lanczos& l)
{
//
// No rational approximations are available because either
@@ -511,16 +511,16 @@ T lgamma_small_imp(T z, T zm1, T zm2, const std::integral_constant<int, 0>&, con
{
// special case near 2:
T dz = zm2;
result = dz * log((z + Lanczos::g() - T(0.5)) / boost::math::constants::e<T>());
result += boost::math::log1p(dz / (Lanczos::g() + T(1.5)), pol) * T(1.5);
result = dz * log((z + lanczos_g_near_1_and_2(l) - T(0.5)) / boost::math::constants::e<T>());
result += boost::math::log1p(dz / (lanczos_g_near_1_and_2(l) + T(1.5)), pol) * T(1.5);
result += boost::math::log1p(Lanczos::lanczos_sum_near_2(dz), pol);
}
else
{
// special case near 1:
T dz = zm1;
result = dz * log((z + Lanczos::g() - T(0.5)) / boost::math::constants::e<T>());
result += boost::math::log1p(dz / (Lanczos::g() + T(0.5)), pol) / 2;
result = dz * log((z + lanczos_g_near_1_and_2(l) - T(0.5)) / boost::math::constants::e<T>());
result += boost::math::log1p(dz / (lanczos_g_near_1_and_2(l) + T(0.5)), pol) / 2;
result += boost::math::log1p(Lanczos::lanczos_sum_near_1(dz), pol);
}
return result;

2223
include/boost/math/special_functions/lanczos.hpp
View File

File diff suppressed because it is too large Load Diff

34
test/test_gamma_mp.cpp
View File

@@ -121,20 +121,6 @@ void expected_results()
// Define the max and mean errors expected for
// various compilers and platforms.
//
add_expected_result(
".*", // compiler
".*", // stdlib
".*", // platform
"number<cpp_bin_float<65> >", // test type(s)
".*", // test data group
"lgamma", 9000, 4000); // test function
add_expected_result(
".*", // compiler
".*", // stdlib
".*", // platform
"number<cpp_bin_float<75> >", // test type(s)
".*", // test data group
"lgamma", 60000, 20000); // test function
add_expected_result(
".*", // compiler
".*", // stdlib
@@ -146,16 +132,30 @@ void expected_results()
".*", // compiler
".*", // stdlib
".*", // platform
".*", // test type(s)
"number<cpp_bin_float<[56]5> >", // test type(s)
".*", // test data group
"lgamma", 4800, 2500); // test function
"lgamma", 7000, 3000); // test function
add_expected_result(
".*", // compiler
".*", // stdlib
".*", // platform
"number<cpp_bin_float<75> >", // test type(s)
".*", // test data group
"lgamma", 40000, 15000); // test function
add_expected_result(
".*", // compiler
".*", // stdlib
".*", // platform
".*", // test type(s)
".*", // test data group
"[tl]gamma", 100, 50); // test function
"lgamma", 600, 200); // test function
add_expected_result(
".*", // compiler
".*", // stdlib
".*", // platform
".*", // test type(s)
".*", // test data group
"[tl]gamma", 120, 50); // test function
//
// Finish off by printing out the compiler/stdlib/platform names,
// we do this to make it easier to mark up expected error rates.

47
tools/lanczos_generator.cpp
View File

@@ -25,7 +25,7 @@
// create a lanczos approximation for that type.
//
#define MP_TYPE boost::multiprecision::number<boost::multiprecision::cpp_bin_float<100> >
#define MP_TYPE boost::multiprecision::number<boost::multiprecision::cpp_bin_float<90> >
//
// this is a sort of recursive include, since this file
@@ -4477,6 +4477,38 @@ lanczos_info<T> generate_lanczos(unsigned n, T g)
return result;
}
template <class T>
T lanczos_conditioning_near_1(const lanczos_info<T>& l)
{
using std::abs;
T sum{ 0 };
T abs_sum{ 0 };
for (unsigned i = 1; i <= l.c.size(); ++i)
{
sum += l.c[i - 1] / (i * i);
abs_sum += abs(l.c[i - 1] / (i * i));
}
return abs(abs_sum / sum);
}
template <class T>
T lanczos_conditioning_near_2(const lanczos_info<T>& l)
{
using std::abs;
T sum{ 0 };
T abs_sum{ 0 };
for (unsigned i = 1; i <= l.c.size(); ++i)
{
sum += l.c[i - 1] / (i * i);
abs_sum += abs(l.c[i - 1] / (1 + 2 * i + i * i));
}
return abs(abs_sum / sum);
}
//
// Stopwatch for measuring performance:
//
@@ -4520,7 +4552,9 @@ std::vector<std::vector<T> > const & get_test_data()
data.back().push_back(static_cast<T>(fact));
data.back().push_back(static_cast<T>(log(fact)));
fact = fact * T(k++);
}while(k < 100);
if (k > 500)
break;
}while((std::numeric_limits<T>::max)() / k > fact);
fact = 0.5;
mp_t srpi = sqrt(boost::math::constants::pi<mp_t>());
@@ -4532,7 +4566,7 @@ std::vector<std::vector<T> > const & get_test_data()
data.back().push_back(static_cast<T>(log(fact*srpi)));
fact *= mul;
mul += 1;
}while(mul < 100);
}while(mul < k);
}
return data;
@@ -5019,7 +5053,7 @@ void find_best_lanczos(const char* name, T eps, int max_scan = 100)
lanczos_info<mp_t> best;
best.err = 100; // best case had better be better than this!
std::cout << std::setw(20) << std::right << "N" << std::setw(20) << std::right << "g" << std::setw(20) << std::right << "eps" << std::setw(20) << std::right << "time (ms)\n";
std::cout << std::setw(20) << std::right << "N" << std::setw(20) << std::right << "g" << std::setw(20) << std::right << "eps" << std::setw(20) << std::right << "c at 1" << std::setw(20) << std::right << "c at 2" << std::setw(20) << std::right << "time (ms)\n";
for (int i = 0; i < sizeof(sweet_spots) / sizeof(sweet_spots[0]); ++i)
{
@@ -5031,8 +5065,11 @@ void find_best_lanczos(const char* name, T eps, int max_scan = 100)
{
std::cout << std::setprecision(14) << std::fixed << std::setw(20) << std::right << sweet_spots[i].N
<< std::setw(20) << std::right << sweet_spots[i].g << std::setw(20) << std::right << static_cast<int>(err / eps)
<< std::setw(20) << std::right << (unsigned)lanczos_conditioning_near_1(info) << std::setw(20) << std::right << (unsigned)lanczos_conditioning_near_2(info)
<< std::setw(20) << std::right << (int)(exec_time * 1000) << std::endl;
}
else
std::cout << "Skipping spot with error: " << std::setprecision(5) << err << std::endl;
if (err < best.err)
{
best = info;
@@ -5042,7 +5079,7 @@ void find_best_lanczos(const char* name, T eps, int max_scan = 100)
}
std::cout << std::endl;
if (best.err < 100)
if (best.err / eps < 100)
print_code(best, name, std::numeric_limits<T>::max_digits10, std::numeric_limits<T>::digits);
else
std::cout << "Sorry, no viable approximation was found!!" << std::endl;