630 lines · plain
1///////////////////////////////////////////////////////////////////////////////2// Copyright 2013 John Maddock3// Distributed under the Boost4// Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_BERNOULLI_DETAIL_HPP8#define BOOST_MATH_BERNOULLI_DETAIL_HPP9 10#include <boost/math/tools/atomic.hpp>11#include <boost/math/tools/toms748_solve.hpp>12#include <boost/math/tools/cxx03_warn.hpp>13#include <boost/math/tools/throw_exception.hpp>14#include <boost/math/tools/config.hpp>15#include <boost/math/special_functions/fpclassify.hpp>16#include <vector>17#include <type_traits>18 19#if defined(BOOST_MATH_HAS_THREADS) && !defined(BOOST_NO_CXX11_HDR_MUTEX) && !defined(BOOST_MATH_NO_ATOMIC_INT)20#include <mutex>21#else22# define BOOST_MATH_BERNOULLI_NOTHREADS23#endif24 25namespace boost{ namespace math{ namespace detail{26//27// Asymptotic expansion for B2n due to28// Luschny LogB3 formula (http://www.luschny.de/math/primes/bernincl.html)29//30template <class T, class Policy>31T b2n_asymptotic(int n)32{33 BOOST_MATH_STD_USING34 const auto nx = static_cast<T>(n);35 const T nx2(nx * nx);36 37 const T approximate_log_of_bernoulli_bn =38 ((boost::math::constants::half<T>() + nx) * log(nx))39 + ((boost::math::constants::half<T>() - nx) * log(boost::math::constants::pi<T>()))40 + (((T(3) / 2) - nx) * boost::math::constants::ln_two<T>())41 + ((nx * (T(2) - (nx2 * 7) * (1 + ((nx2 * 30) * ((nx2 * 12) - 1))))) / (((nx2 * nx2) * nx2) * 2520));42 return ((n / 2) & 1 ? 1 : -1) * (approximate_log_of_bernoulli_bn > tools::log_max_value<T>()43 ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, nx, Policy())44 : static_cast<T>(exp(approximate_log_of_bernoulli_bn)));45}46 47template <class T, class Policy>48T t2n_asymptotic(int n)49{50 BOOST_MATH_STD_USING51 // Just get B2n and convert to a Tangent number:52 T t2n = fabs(b2n_asymptotic<T, Policy>(2 * n)) / (2 * n);53 T p2 = ldexp(T(1), n);54 if(tools::max_value<T>() / p2 < t2n)55 {56 return policies::raise_overflow_error<T>("boost::math::tangent_t2n<%1%>(std::size_t)", nullptr, T(n), Policy());57 }58 t2n *= p2;59 p2 -= 1;60 if(tools::max_value<T>() / p2 < t2n)61 {62 return policies::raise_overflow_error<T>("boost::math::tangent_t2n<%1%>(std::size_t)", nullptr, Policy());63 }64 t2n *= p2;65 return t2n;66}67//68// We need to know the approximate value of /n/ which will69// cause bernoulli_b2n<T>(n) to return infinity - this allows70// us to elude a great deal of runtime checking for values below71// n, and only perform the full overflow checks when we know that we're72// getting close to the point where our calculations will overflow.73// We use Luschny's LogB3 formula (http://www.luschny.de/math/primes/bernincl.html)74// to find the limit, and since we're dealing with the log of the Bernoulli numbers75// we need only perform the calculation at double precision and not with T76// (which may be a multiprecision type). The limit returned is within 1 of the true77// limit for all the types tested. Note that although the code below is basically78// the same as b2n_asymptotic above, it has been recast as a continuous real-valued79// function as this makes the root finding go smoother/faster. It also omits the80// sign of the Bernoulli number.81//82struct max_bernoulli_root_functor83{84 explicit max_bernoulli_root_functor(unsigned long long t) : target(static_cast<double>(t)) {}85 double operator()(double n) const86 {87 BOOST_MATH_STD_USING88 89 // Luschny LogB3(n) formula.90 91 const double nx2(n * n);92 93 const double approximate_log_of_bernoulli_bn94 = ((boost::math::constants::half<double>() + n) * log(n))95 + ((boost::math::constants::half<double>() - n) * log(boost::math::constants::pi<double>()))96 + (((static_cast<double>(3) / 2) - n) * boost::math::constants::ln_two<double>())97 + ((n * (2 - (nx2 * 7) * (1 + ((nx2 * 30) * ((nx2 * 12) - 1))))) / (((nx2 * nx2) * nx2) * 2520));98 99 return approximate_log_of_bernoulli_bn - target;100 }101private:102 double target;103};104 105template <class T, class Policy>106inline std::size_t find_bernoulli_overflow_limit(const std::false_type&)107{108 // Set a limit on how large the result can ever be:109 static const auto max_result = static_cast<double>((std::numeric_limits<std::size_t>::max)() - 1000u);110 111 unsigned long long t = static_cast<unsigned long long>(lltrunc(boost::math::tools::log_max_value<T>()));112 max_bernoulli_root_functor fun(t);113 boost::math::tools::equal_floor tol;114 std::uintmax_t max_iter = boost::math::policies::get_max_root_iterations<Policy>();115 double result = boost::math::tools::toms748_solve(fun, sqrt(static_cast<double>(t)), static_cast<double>(t), tol, max_iter).first / 2;116 if (result > max_result)117 {118 result = max_result;119 }120 121 return static_cast<std::size_t>(result);122}123 124template <class T, class Policy>125inline std::size_t find_bernoulli_overflow_limit(const std::true_type&)126{127 return max_bernoulli_index<bernoulli_imp_variant<T>::value>::value;128}129 130template <class T, class Policy>131std::size_t b2n_overflow_limit()132{133 // This routine is called at program startup if it's called at all:134 // that guarantees safe initialization of the static variable.135 using tag_type = std::integral_constant<bool, (bernoulli_imp_variant<T>::value >= 1) && (bernoulli_imp_variant<T>::value <= 3)>;136 static const std::size_t lim = find_bernoulli_overflow_limit<T, Policy>(tag_type());137 return lim;138}139 140//141// The tangent numbers grow larger much more rapidly than the Bernoulli numbers do....142// so to compute the Bernoulli numbers from the tangent numbers, we need to avoid spurious143// overflow in the calculation, we can do this by scaling all the tangent number by some scale factor:144//145template <class T, typename std::enable_if<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), bool>::type = true>146inline T tangent_scale_factor()147{148 BOOST_MATH_STD_USING149 return ldexp(T(1), std::numeric_limits<T>::min_exponent + 5);150}151 152template <class T, typename std::enable_if<!std::numeric_limits<T>::is_specialized || !(std::numeric_limits<T>::radix == 2), bool>::type = true>153inline T tangent_scale_factor()154{155 return tools::min_value<T>() * 16;156}157 158//159// We need something to act as a cache for our calculated Bernoulli numbers. In order to160// ensure both fast access and thread safety, we need a stable table which may be extended161// in size, but which never reallocates: that way values already calculated may be accessed162// concurrently with another thread extending the table with new values.163//164// Very very simple vector class that will never allocate more than once, we could use165// boost::container::static_vector here, but that allocates on the stack, which may well166// cause issues for the amount of memory we want in the extreme case...167//168template <class T>169struct fixed_vector : private std::allocator<T>170{171 using size_type = unsigned;172 using iterator = T*;173 using const_iterator = const T*;174 fixed_vector() : m_used(0)175 {176 std::size_t overflow_limit = 5 + b2n_overflow_limit<T, policies::policy<> >();177 m_capacity = static_cast<unsigned>((std::min)(overflow_limit, static_cast<std::size_t>(100000u)));178 m_data = this->allocate(m_capacity);179 }180 ~fixed_vector()181 {182 using allocator_type = std::allocator<T>;183 using allocator_traits = std::allocator_traits<allocator_type>;184 allocator_type& alloc = *this;185 for(unsigned i = 0; i < m_used; ++i)186 {187 allocator_traits::destroy(alloc, &m_data[i]);188 }189 allocator_traits::deallocate(alloc, m_data, m_capacity);190 }191 T& operator[](unsigned n) { BOOST_MATH_ASSERT(n < m_used); return m_data[n]; }192 const T& operator[](unsigned n)const { BOOST_MATH_ASSERT(n < m_used); return m_data[n]; }193 unsigned size()const { return m_used; }194 unsigned size() { return m_used; }195 bool resize(unsigned n, const T& val)196 {197 if(n > m_capacity)198 {199#ifndef BOOST_MATH_NO_EXCEPTIONS200 BOOST_MATH_THROW_EXCEPTION(std::runtime_error("Exhausted storage for Bernoulli numbers."));201#else202 return false;203#endif204 }205 for(unsigned i = m_used; i < n; ++i)206 new (m_data + i) T(val);207 m_used = n;208 return true;209 }210 bool resize(unsigned n) { return resize(n, T()); }211 T* begin() { return m_data; }212 T* end() { return m_data + m_used; }213 T* begin()const { return m_data; }214 T* end()const { return m_data + m_used; }215 unsigned capacity()const { return m_capacity; }216 void clear() { m_used = 0; }217private:218 T* m_data;219 unsigned m_used {};220 unsigned m_capacity;221};222 223template <class T, class Policy>224class bernoulli_numbers_cache225{226public:227 bernoulli_numbers_cache() : m_overflow_limit((std::numeric_limits<std::size_t>::max)())228 , m_counter(0)229 , m_current_precision(boost::math::tools::digits<T>())230 {}231 232 using container_type = fixed_vector<T>;233 234 bool tangent(std::size_t m)235 {236 static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1;237 238 if (!tn.resize(static_cast<typename container_type::size_type>(m), T(0U)))239 {240 return false;241 }242 243 BOOST_MATH_INSTRUMENT_VARIABLE(min_overflow_index);244 245 std::size_t prev_size = m_intermediates.size();246 m_intermediates.resize(m, T(0U));247 248 if(prev_size == 0)249 {250 m_intermediates[1] = tangent_scale_factor<T>() /*T(1U)*/;251 tn[0U] = T(0U);252 tn[1U] = tangent_scale_factor<T>()/* T(1U)*/;253 BOOST_MATH_INSTRUMENT_VARIABLE(tn[0]);254 BOOST_MATH_INSTRUMENT_VARIABLE(tn[1]);255 }256 257 for(std::size_t i = std::max<size_t>(2, prev_size); i < m; i++)258 {259 bool overflow_check = false;260 if(i >= min_overflow_index && (boost::math::tools::max_value<T>() / (i-1) < m_intermediates[1]) )261 {262 std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>());263 break;264 }265 m_intermediates[1] = m_intermediates[1] * (i-1);266 for(std::size_t j = 2; j <= i; j++)267 {268 overflow_check =269 (i >= min_overflow_index) && (270 (boost::math::tools::max_value<T>() / (i - j) < m_intermediates[j])271 || (boost::math::tools::max_value<T>() / (i - j + 2) < m_intermediates[j-1])272 || (boost::math::tools::max_value<T>() - m_intermediates[j] * (i - j) < m_intermediates[j-1] * (i - j + 2))273 || ((boost::math::isinf)(m_intermediates[j]))274 );275 276 if(overflow_check)277 {278 std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>());279 break;280 }281 m_intermediates[j] = m_intermediates[j] * (i - j) + m_intermediates[j-1] * (i - j + 2);282 }283 if(overflow_check)284 break; // already filled the tn...285 tn[static_cast<typename container_type::size_type>(i)] = m_intermediates[i];286 BOOST_MATH_INSTRUMENT_VARIABLE(i);287 BOOST_MATH_INSTRUMENT_VARIABLE(tn[static_cast<typename container_type::size_type>(i)]);288 }289 return true;290 }291 292 bool tangent_numbers_series(const std::size_t m)293 {294 BOOST_MATH_STD_USING295 static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1;296 297 typename container_type::size_type old_size = bn.size();298 299 if (!tangent(m))300 return false;301 if (!bn.resize(static_cast<typename container_type::size_type>(m)))302 return false;303 304 if(!old_size)305 {306 bn[0] = 1;307 old_size = 1;308 }309 310 T power_two(ldexp(T(1), static_cast<int>(2 * old_size)));311 312 for(std::size_t i = old_size; i < m; i++)313 {314 T b(static_cast<T>(i * 2));315 //316 // Not only do we need to take care to avoid spurious over/under flow in317 // the calculation, but we also need to avoid overflow altogether in case318 // we're calculating with a type where "bad things" happen in that case:319 //320 b = b / (power_two * tangent_scale_factor<T>());321 b /= (power_two - 1);322 bool overflow_check = (i >= min_overflow_index) && (tools::max_value<T>() / tn[static_cast<typename container_type::size_type>(i)] < b);323 if(overflow_check)324 {325 m_overflow_limit = i;326 while(i < m)327 {328 b = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : tools::max_value<T>();329 bn[static_cast<typename container_type::size_type>(i)] = ((i % 2U) ? b : T(-b));330 ++i;331 }332 break;333 }334 else335 {336 b *= tn[static_cast<typename container_type::size_type>(i)];337 }338 339 power_two = ldexp(power_two, 2);340 341 const bool b_neg = i % 2 == 0;342 343 bn[static_cast<typename container_type::size_type>(i)] = ((!b_neg) ? b : T(-b));344 }345 return true;346 }347 348 template <class OutputIterator>349 OutputIterator copy_bernoulli_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol)350 {351 //352 // There are basically 3 thread safety options:353 //354 // 1) There are no threads (BOOST_MATH_HAS_THREADS is not defined).355 // 2) There are threads, but we do not have a true atomic integer type,356 // in this case we just use a mutex to guard against race conditions.357 // 3) There are threads, and we have an atomic integer: in this case we can358 // use the double-checked locking pattern to avoid thread synchronisation359 // when accessing values already in the cache.360 //361 // First off handle the common case for overflow and/or asymptotic expansion:362 //363 if(start + n > bn.capacity())364 {365 if(start < bn.capacity())366 {367 out = copy_bernoulli_numbers(out, start, bn.capacity() - start, pol);368 n -= bn.capacity() - start;369 start = static_cast<std::size_t>(bn.capacity());370 }371 if(start < b2n_overflow_limit<T, Policy>() + 2u)372 {373 for(; n; ++start, --n)374 {375 *out = b2n_asymptotic<T, Policy>(static_cast<int>(start * 2U));376 ++out;377 }378 }379 for(; n; ++start, --n)380 {381 *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(start), pol);382 ++out;383 }384 return out;385 }386 387 #if defined(BOOST_MATH_HAS_THREADS) && defined(BOOST_MATH_BERNOULLI_NOTHREADS) && !defined(BOOST_MATH_BERNOULLI_UNTHREADED)388 // Add a static_assert on instantiation if we have threads, but no C++11 threading support.389 static_assert(sizeof(T) == 1, "Unsupported configuration: your platform appears to have either no atomic integers, or no std::mutex. If you are happy with thread-unsafe code, then you may define BOOST_MATH_BERNOULLI_UNTHREADED to suppress this error.");390 #elif defined(BOOST_MATH_BERNOULLI_NOTHREADS)391 //392 // Single threaded code, very simple:393 //394 if(m_current_precision < boost::math::tools::digits<T>())395 {396 bn.clear();397 tn.clear();398 m_intermediates.clear();399 m_current_precision = boost::math::tools::digits<T>();400 }401 if(start + n >= bn.size())402 {403 std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));404 if (!tangent_numbers_series(new_size))405 {406 return std::fill_n(out, n, policies::raise_evaluation_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", "Unable to allocate Bernoulli numbers cache for %1% values", T(start + n), pol));407 }408 }409 410 for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)411 {412 *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol) : bn[i];413 ++out;414 }415 #else416 //417 // Double-checked locking pattern, lets us access cached already cached values418 // without locking:419 //420 // Get the counter and see if we need to calculate more constants:421 //422 if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)423 || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))424 {425 std::lock_guard<std::mutex> l(m_mutex);426 427 if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)428 || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))429 {430 if(static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>())431 {432 bn.clear();433 tn.clear();434 m_intermediates.clear();435 m_counter.store(0, std::memory_order_release);436 m_current_precision = boost::math::tools::digits<T>();437 }438 if(start + n >= bn.size())439 {440 std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));441 if (!tangent_numbers_series(new_size))442 return std::fill_n(out, n, policies::raise_evaluation_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", "Unable to allocate Bernoulli numbers cache for %1% values", T(new_size), pol));443 }444 m_counter.store(static_cast<atomic_integer_type>(bn.size()), std::memory_order_release);445 }446 }447 448 for(std::size_t i = (std::max)(static_cast<std::size_t>(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)449 {450 *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol) : bn[static_cast<typename container_type::size_type>(i)];451 ++out;452 }453 454 #endif // BOOST_MATH_HAS_THREADS455 return out;456 }457 458 template <class OutputIterator>459 OutputIterator copy_tangent_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol)460 {461 //462 // There are basically 3 thread safety options:463 //464 // 1) There are no threads (BOOST_MATH_HAS_THREADS is not defined).465 // 2) There are threads, but we do not have a true atomic integer type,466 // in this case we just use a mutex to guard against race conditions.467 // 3) There are threads, and we have an atomic integer: in this case we can468 // use the double-checked locking pattern to avoid thread synchronisation469 // when accessing values already in the cache.470 //471 //472 // First off handle the common case for overflow and/or asymptotic expansion:473 //474 if(start + n > bn.capacity())475 {476 if(start < bn.capacity())477 {478 out = copy_tangent_numbers(out, start, bn.capacity() - start, pol);479 n -= bn.capacity() - start;480 start = static_cast<std::size_t>(bn.capacity());481 }482 if(start < b2n_overflow_limit<T, Policy>() + 2u)483 {484 for(; n; ++start, --n)485 {486 *out = t2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start));487 ++out;488 }489 }490 for(; n; ++start, --n)491 {492 *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(start), pol);493 ++out;494 }495 return out;496 }497 498 #if defined(BOOST_MATH_BERNOULLI_NOTHREADS)499 //500 // Single threaded code, very simple:501 //502 if(m_current_precision < boost::math::tools::digits<T>())503 {504 bn.clear();505 tn.clear();506 m_intermediates.clear();507 m_current_precision = boost::math::tools::digits<T>();508 }509 if(start + n >= bn.size())510 {511 std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));512 if (!tangent_numbers_series(new_size))513 return std::fill_n(out, n, policies::raise_evaluation_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", "Unable to allocate Bernoulli numbers cache for %1% values", T(start + n), pol));514 }515 516 for(std::size_t i = start; i < start + n; ++i)517 {518 if(i >= m_overflow_limit)519 *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol);520 else521 {522 if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])523 *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol);524 else525 *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();526 }527 ++out;528 }529 #elif defined(BOOST_MATH_NO_ATOMIC_INT)530 static_assert(sizeof(T) == 1, "Unsupported configuration: your platform appears to have no atomic integers. If you are happy with thread-unsafe code, then you may define BOOST_MATH_BERNOULLI_UNTHREADED to suppress this error.");531 #else532 //533 // Double-checked locking pattern, lets us access cached already cached values534 // without locking:535 //536 // Get the counter and see if we need to calculate more constants:537 //538 if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)539 || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))540 {541 std::lock_guard<std::mutex> l(m_mutex);542 543 if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)544 || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))545 {546 if(static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>())547 {548 bn.clear();549 tn.clear();550 m_intermediates.clear();551 m_counter.store(0, std::memory_order_release);552 m_current_precision = boost::math::tools::digits<T>();553 }554 if(start + n >= bn.size())555 {556 std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));557 if (!tangent_numbers_series(new_size))558 return std::fill_n(out, n, policies::raise_evaluation_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", "Unable to allocate Bernoulli numbers cache for %1% values", T(start + n), pol));559 }560 m_counter.store(static_cast<atomic_integer_type>(bn.size()), std::memory_order_release);561 }562 }563 564 for(std::size_t i = start; i < start + n; ++i)565 {566 if(i >= m_overflow_limit)567 *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol);568 else569 {570 if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])571 *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol);572 else573 *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();574 }575 ++out;576 }577 578 #endif // BOOST_MATH_HAS_THREADS579 return out;580 }581 582private:583 //584 // The caches for Bernoulli and tangent numbers, once allocated,585 // these must NEVER EVER reallocate as it breaks our thread586 // safety guarantees:587 //588 fixed_vector<T> bn, tn;589 std::vector<T> m_intermediates;590 // The value at which we know overflow has already occurred for the Bn:591 std::size_t m_overflow_limit;592 593 #if !defined(BOOST_MATH_BERNOULLI_NOTHREADS)594 std::mutex m_mutex;595 atomic_counter_type m_counter, m_current_precision;596 #else597 int m_counter;598 int m_current_precision;599 #endif // BOOST_MATH_HAS_THREADS600};601 602template <class T, class Policy>603inline 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()604{605 //606 // When numeric_limits<>::digits is zero, the type has either not specialized numeric_limits at all607 // or it's precision can vary at runtime. So make the cache thread_local so that each thread can608 // have it's own precision if required:609 //610 static611#ifndef BOOST_MATH_NO_THREAD_LOCAL_WITH_NON_TRIVIAL_TYPES612 BOOST_MATH_THREAD_LOCAL613#endif614 bernoulli_numbers_cache<T, Policy> data;615 return data;616}617template <class T, class Policy>618inline 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()619{620 //621 // Note that we rely on C++11 thread-safe initialization here:622 //623 static bernoulli_numbers_cache<T, Policy> data;624 return data;625}626 627}}}628 629#endif // BOOST_MATH_BERNOULLI_DETAIL_HPP630