#ifndef BOOST_MATH_BERNOULLI_DETAIL_HPP
#define BOOST_MATH_BERNOULLI_DETAIL_HPP
-#include <boost/config.hpp>
-#include <boost/detail/lightweight_mutex.hpp>
#include <boost/math/tools/atomic.hpp>
-#include <boost/utility/enable_if.hpp>
#include <boost/math/tools/toms748_solve.hpp>
#include <boost/math/tools/cxx03_warn.hpp>
+#include <boost/math/tools/throw_exception.hpp>
+#include <boost/math/tools/config.hpp>
+#include <boost/math/special_functions/fpclassify.hpp>
#include <vector>
+#include <type_traits>
+
+#if defined(BOOST_HAS_THREADS) && !defined(BOOST_NO_CXX11_HDR_MUTEX) && !defined(BOOST_MATH_NO_ATOMIC_INT)
+#include <mutex>
+#else
+# define BOOST_MATH_BERNOULLI_NOTHREADS
+#endif
namespace boost{ namespace math{ namespace detail{
//
//
struct max_bernoulli_root_functor
{
- max_bernoulli_root_functor(ulong_long_type t) : target(static_cast<double>(t)) {}
+ max_bernoulli_root_functor(unsigned long long t) : target(static_cast<double>(t)) {}
double operator()(double n)
{
BOOST_MATH_STD_USING
};
template <class T, class Policy>
-inline std::size_t find_bernoulli_overflow_limit(const boost::false_type&)
+inline std::size_t find_bernoulli_overflow_limit(const std::false_type&)
{
// Set a limit on how large the result can ever be:
static const double max_result = static_cast<double>((std::numeric_limits<std::size_t>::max)() - 1000u);
- ulong_long_type t = lltrunc(boost::math::tools::log_max_value<T>());
+ unsigned long long t = lltrunc(boost::math::tools::log_max_value<T>());
max_bernoulli_root_functor fun(t);
boost::math::tools::equal_floor tol;
- boost::uintmax_t max_iter = boost::math::policies::get_max_root_iterations<Policy>();
+ std::uintmax_t max_iter = boost::math::policies::get_max_root_iterations<Policy>();
double result = boost::math::tools::toms748_solve(fun, sqrt(double(t)), double(t), tol, max_iter).first / 2;
if (result > max_result)
result = max_result;
}
template <class T, class Policy>
-inline std::size_t find_bernoulli_overflow_limit(const boost::true_type&)
+inline std::size_t find_bernoulli_overflow_limit(const std::true_type&)
{
return max_bernoulli_index<bernoulli_imp_variant<T>::value>::value;
}
{
// This routine is called at program startup if it's called at all:
// that guarantees safe initialization of the static variable.
- typedef boost::integral_constant<bool, (bernoulli_imp_variant<T>::value >= 1) && (bernoulli_imp_variant<T>::value <= 3)> tag_type;
+ typedef std::integral_constant<bool, (bernoulli_imp_variant<T>::value >= 1) && (bernoulli_imp_variant<T>::value <= 3)> tag_type;
static const std::size_t lim = find_bernoulli_overflow_limit<T, Policy>(tag_type());
return lim;
}
// so to compute the Bernoulli numbers from the tangent numbers, we need to avoid spurious
// overflow in the calculation, we can do this by scaling all the tangent number by some scale factor:
//
-template <class T>
-inline typename enable_if_c<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), T>::type tangent_scale_factor()
+template <class T, typename std::enable_if<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), bool>::type = true>
+inline T tangent_scale_factor()
{
BOOST_MATH_STD_USING
return ldexp(T(1), std::numeric_limits<T>::min_exponent + 5);
}
-template <class T>
-inline typename disable_if_c<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), T>::type tangent_scale_factor()
+
+template <class T, typename std::enable_if<!std::numeric_limits<T>::is_specialized || !(std::numeric_limits<T>::radix == 2), bool>::type = true>
+inline T tangent_scale_factor()
{
return tools::min_value<T>() * 16;
}
-//
-// Initializer: ensure all our constants are initialized prior to the first call of main:
-//
-template <class T, class Policy>
-struct bernoulli_initializer
-{
- struct init
- {
- init()
- {
- //
- // We call twice, once to initialize our static table, and once to
- // initialize our dymanic table:
- //
- boost::math::bernoulli_b2n<T>(2, Policy());
-#ifndef BOOST_NO_EXCEPTIONS
- try{
-#endif
- boost::math::bernoulli_b2n<T>(max_bernoulli_b2n<T>::value + 1, Policy());
-#ifndef BOOST_NO_EXCEPTIONS
- } catch(const std::overflow_error&){}
-#endif
- boost::math::tangent_t2n<T>(2, Policy());
- }
- void force_instantiate()const{}
- };
- static const init initializer;
- static void force_instantiate()
- {
- initializer.force_instantiate();
- }
-};
-
-template <class T, class Policy>
-const typename bernoulli_initializer<T, Policy>::init bernoulli_initializer<T, Policy>::initializer;
//
// We need something to act as a cache for our calculated Bernoulli numbers. In order to
}
~fixed_vector()
{
-#ifdef BOOST_NO_CXX11_ALLOCATOR
- for(unsigned i = 0; i < m_used; ++i)
- this->destroy(&m_data[i]);
- this->deallocate(m_data, m_capacity);
-#else
typedef std::allocator<T> allocator_type;
typedef std::allocator_traits<allocator_type> allocator_traits;
allocator_type& alloc = *this;
for(unsigned i = 0; i < m_used; ++i)
allocator_traits::destroy(alloc, &m_data[i]);
allocator_traits::deallocate(alloc, m_data, m_capacity);
-#endif
}
- T& operator[](unsigned n) { BOOST_ASSERT(n < m_used); return m_data[n]; }
- const T& operator[](unsigned n)const { BOOST_ASSERT(n < m_used); return m_data[n]; }
+ T& operator[](unsigned n) { BOOST_MATH_ASSERT(n < m_used); return m_data[n]; }
+ const T& operator[](unsigned n)const { BOOST_MATH_ASSERT(n < m_used); return m_data[n]; }
unsigned size()const { return m_used; }
unsigned size() { return m_used; }
void resize(unsigned n, const T& val)
{
if(n > m_capacity)
{
- BOOST_THROW_EXCEPTION(std::runtime_error("Exhausted storage for Bernoulli numbers."));
+ BOOST_MATH_THROW_EXCEPTION(std::runtime_error("Exhausted storage for Bernoulli numbers."));
}
for(unsigned i = m_used; i < n; ++i)
new (m_data + i) T(val);
{
public:
bernoulli_numbers_cache() : m_overflow_limit((std::numeric_limits<std::size_t>::max)())
-#if defined(BOOST_HAS_THREADS) && !defined(BOOST_MATH_NO_ATOMIC_INT)
, m_counter(0)
-#endif
, m_current_precision(boost::math::tools::digits<T>())
{}
}
return out;
}
- #if !defined(BOOST_HAS_THREADS)
- //
- // Single threaded code, very simple:
- //
- if(m_current_precision < boost::math::tools::digits<T>())
- {
- bn.clear();
- tn.clear();
- m_intermediates.clear();
- m_current_precision = boost::math::tools::digits<T>();
- }
- if(start + n >= bn.size())
- {
- 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()));
- tangent_numbers_series(new_size);
- }
- for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
- {
- *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
- ++out;
- }
- #elif defined(BOOST_MATH_NO_ATOMIC_INT)
+ #if defined(BOOST_HAS_THREADS) && defined(BOOST_MATH_BERNOULLI_NOTHREADS) && !defined(BOOST_MATH_BERNOULLI_UNTHREADED)
+ // Add a static_assert on instantiation if we have threads, but no C++11 threading support.
+ 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.");
+ #elif defined(BOOST_MATH_BERNOULLI_NOTHREADS)
//
- // We need to grab a mutex every time we get here, for both readers and writers:
+ // Single threaded code, very simple:
//
- boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
if(m_current_precision < boost::math::tools::digits<T>())
{
bn.clear();
*out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
++out;
}
-
- #else
+ #else
//
// Double-checked locking pattern, lets us access cached already cached values
// without locking:
//
// Get the counter and see if we need to calculate more constants:
//
- if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
- || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()))
+ if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)
+ || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))
{
- boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
+ std::lock_guard<std::mutex> l(m_mutex);
- if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
- || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()))
+ if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)
+ || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))
{
- if(static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())
+ if(static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>())
{
bn.clear();
tn.clear();
m_intermediates.clear();
- m_counter.store(0, BOOST_MATH_ATOMIC_NS::memory_order_release);
+ m_counter.store(0, std::memory_order_release);
m_current_precision = boost::math::tools::digits<T>();
}
if(start + n >= bn.size())
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()));
tangent_numbers_series(new_size);
}
- m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release);
+ m_counter.store(static_cast<atomic_integer_type>(bn.size()), std::memory_order_release);
}
}
++out;
}
- #endif
+ #endif // BOOST_HAS_THREADS
return out;
}
}
return out;
}
- #if !defined(BOOST_HAS_THREADS)
- //
- // Single threaded code, very simple:
- //
- if(m_current_precision < boost::math::tools::digits<T>())
- {
- bn.clear();
- tn.clear();
- m_intermediates.clear();
- m_current_precision = boost::math::tools::digits<T>();
- }
- if(start + n >= bn.size())
- {
- 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()));
- tangent_numbers_series(new_size);
- }
- for(std::size_t i = start; i < start + n; ++i)
- {
- if(i >= m_overflow_limit)
- *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
- else
- {
- if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
- *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
- else
- *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
- }
- ++out;
- }
- #elif defined(BOOST_MATH_NO_ATOMIC_INT)
+ #if defined(BOOST_MATH_BERNOULLI_NOTHREADS)
//
- // We need to grab a mutex every time we get here, for both readers and writers:
+ // Single threaded code, very simple:
//
- boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
if(m_current_precision < boost::math::tools::digits<T>())
{
bn.clear();
}
++out;
}
-
- #else
+ #elif defined(BOOST_MATH_NO_ATOMIC_INT)
+ 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.");
+ #else
//
// Double-checked locking pattern, lets us access cached already cached values
// without locking:
//
// Get the counter and see if we need to calculate more constants:
//
- if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
- || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()))
+ if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)
+ || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))
{
- boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
+ std::lock_guard<std::mutex> l(m_mutex);
- if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
- || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()))
+ if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)
+ || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))
{
- if(static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())
+ if(static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>())
{
bn.clear();
tn.clear();
m_intermediates.clear();
- m_counter.store(0, BOOST_MATH_ATOMIC_NS::memory_order_release);
+ m_counter.store(0, std::memory_order_release);
m_current_precision = boost::math::tools::digits<T>();
}
if(start + n >= bn.size())
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()));
tangent_numbers_series(new_size);
}
- m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release);
+ m_counter.store(static_cast<atomic_integer_type>(bn.size()), std::memory_order_release);
}
}
++out;
}
- #endif
+ #endif // BOOST_HAS_THREADS
return out;
}
std::vector<T> m_intermediates;
// The value at which we know overflow has already occurred for the Bn:
std::size_t m_overflow_limit;
-#if !defined(BOOST_HAS_THREADS)
- int m_current_precision;
-#elif defined(BOOST_MATH_NO_ATOMIC_INT)
- boost::detail::lightweight_mutex m_mutex;
- int m_current_precision;
-#else
- boost::detail::lightweight_mutex m_mutex;
+
+ #if !defined(BOOST_MATH_BERNOULLI_NOTHREADS)
+ std::mutex m_mutex;
atomic_counter_type m_counter, m_current_precision;
-#endif
+ #else
+ int m_counter;
+ int m_current_precision;
+ #endif // BOOST_HAS_THREADS
};
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()
+{
+ //
+ // 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
+ BOOST_MATH_THREAD_LOCAL
+#endif
+ 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()
{
//
- // Force this function to be called at program startup so all the static variables
- // get initialized then (thread safety).
+ // Note that we rely on C++11 thread-safe initialization here:
//
- bernoulli_initializer<T, Policy>::force_instantiate();
static bernoulli_numbers_cache<T, Policy> data;
return data;
}