special_functions/detail/bernoulli_details.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //  Copyright 2013 John Maddock
0003 //  Distributed under the Boost
0004 //  Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_BERNOULLI_DETAIL_HPP
0008 #define BOOST_MATH_BERNOULLI_DETAIL_HPP
0009
0010 #include <boost/math/tools/atomic.hpp>
0011 #include <boost/math/tools/toms748_solve.hpp>
0012 #include <boost/math/tools/cxx03_warn.hpp>
0013 #include <boost/math/tools/throw_exception.hpp>
0014 #include <boost/math/tools/config.hpp>
0015 #include <boost/math/special_functions/fpclassify.hpp>
0016 #include <vector>
0017 #include <type_traits>
0018
0019 #if defined(BOOST_HAS_THREADS) && !defined(BOOST_NO_CXX11_HDR_MUTEX) && !defined(BOOST_MATH_NO_ATOMIC_INT)
0020 #include <mutex>
0021 #else
0022 #  define BOOST_MATH_BERNOULLI_NOTHREADS
0023 #endif
0024
0025 namespace boost{ namespace math{ namespace detail{
0026 //
0027 // Asymptotic expansion for B2n due to
0028 // Luschny LogB3 formula (http://www.luschny.de/math/primes/bernincl.html)
0029 //
0030 template <class T, class Policy>
0031 T b2n_asymptotic(int n)
0032 {
0033    BOOST_MATH_STD_USING
0034    const auto nx = static_cast<T>(n);
0035    const T nx2(nx * nx);
0036
0037    const T approximate_log_of_bernoulli_bn =
0038         ((boost::math::constants::half<T>() + nx) * log(nx))
0039         + ((boost::math::constants::half<T>() - nx) * log(boost::math::constants::pi<T>()))
0040         + (((T(3) / 2) - nx) * boost::math::constants::ln_two<T>())
0041         + ((nx * (T(2) - (nx2 * 7) * (1 + ((nx2 * 30) * ((nx2 * 12) - 1))))) / (((nx2 * nx2) * nx2) * 2520));
0042    return ((n / 2) & 1 ? 1 : -1) * (approximate_log_of_bernoulli_bn > tools::log_max_value<T>()
0043       ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, nx, Policy())
0044       : static_cast<T>(exp(approximate_log_of_bernoulli_bn)));
0045 }
0046
0047 template <class T, class Policy>
0048 T t2n_asymptotic(int n)
0049 {
0050    BOOST_MATH_STD_USING
0051    // Just get B2n and convert to a Tangent number:
0052    T t2n = fabs(b2n_asymptotic<T, Policy>(2 * n)) / (2 * n);
0053    T p2 = ldexp(T(1), n);
0054    if(tools::max_value<T>() / p2 < t2n)
0055    {
0056       return policies::raise_overflow_error<T>("boost::math::tangent_t2n<%1%>(std::size_t)", nullptr, T(n), Policy());
0057    }
0058    t2n *= p2;
0059    p2 -= 1;
0060    if(tools::max_value<T>() / p2 < t2n)
0061    {
0062       return policies::raise_overflow_error<T>("boost::math::tangent_t2n<%1%>(std::size_t)", nullptr, Policy());
0063    }
0064    t2n *= p2;
0065    return t2n;
0066 }
0067 //
0068 // We need to know the approximate value of /n/ which will
0069 // cause bernoulli_b2n<T>(n) to return infinity - this allows
0070 // us to elude a great deal of runtime checking for values below
0071 // n, and only perform the full overflow checks when we know that we're
0072 // getting close to the point where our calculations will overflow.
0073 // We use Luschny's LogB3 formula (http://www.luschny.de/math/primes/bernincl.html)
0074 // to find the limit, and since we're dealing with the log of the Bernoulli numbers
0075 // we need only perform the calculation at double precision and not with T
0076 // (which may be a multiprecision type).  The limit returned is within 1 of the true
0077 // limit for all the types tested.  Note that although the code below is basically
0078 // the same as b2n_asymptotic above, it has been recast as a continuous real-valued
0079 // function as this makes the root finding go smoother/faster.  It also omits the
0080 // sign of the Bernoulli number.
0081 //
0082 struct max_bernoulli_root_functor
0083 {
0084    explicit max_bernoulli_root_functor(unsigned long long t) : target(static_cast<double>(t)) {}
0085    double operator()(double n) const
0086    {
0087       BOOST_MATH_STD_USING
0088
0089       // Luschny LogB3(n) formula.
0090
0091       const double nx2(n * n);
0092
0093       const double approximate_log_of_bernoulli_bn
0094          =   ((boost::math::constants::half<double>() + n) * log(n))
0095            + ((boost::math::constants::half<double>() - n) * log(boost::math::constants::pi<double>()))
0096            + (((static_cast<double>(3) / 2) - n) * boost::math::constants::ln_two<double>())
0097            + ((n * (2 - (nx2 * 7) * (1 + ((nx2 * 30) * ((nx2 * 12) - 1))))) / (((nx2 * nx2) * nx2) * 2520));
0098
0099       return approximate_log_of_bernoulli_bn - target;
0100    }
0101 private:
0102    double target;
0103 };
0104
0105 template <class T, class Policy>
0106 inline std::size_t find_bernoulli_overflow_limit(const std::false_type&)
0107 {
0108    // Set a limit on how large the result can ever be:
0109    static const auto max_result = static_cast<double>((std::numeric_limits<std::size_t>::max)() - 1000u);
0110
0111    unsigned long long t = lltrunc(boost::math::tools::log_max_value<T>());
0112    max_bernoulli_root_functor fun(t);
0113    boost::math::tools::equal_floor tol;
0114    std::uintmax_t max_iter = boost::math::policies::get_max_root_iterations<Policy>();
0115    double result = boost::math::tools::toms748_solve(fun, sqrt(static_cast<double>(t)), static_cast<double>(t), tol, max_iter).first / 2;
0116    if (result > max_result)
0117    {
0118       result = max_result;
0119    }
0120
0121    return static_cast<std::size_t>(result);
0122 }
0123
0124 template <class T, class Policy>
0125 inline std::size_t find_bernoulli_overflow_limit(const std::true_type&)
0126 {
0127    return max_bernoulli_index<bernoulli_imp_variant<T>::value>::value;
0128 }
0129
0130 template <class T, class Policy>
0131 std::size_t b2n_overflow_limit()
0132 {
0133    // This routine is called at program startup if it's called at all:
0134    // that guarantees safe initialization of the static variable.
0135    using tag_type = std::integral_constant<bool, (bernoulli_imp_variant<T>::value >= 1) && (bernoulli_imp_variant<T>::value <= 3)>;
0136    static const std::size_t lim = find_bernoulli_overflow_limit<T, Policy>(tag_type());
0137    return lim;
0138 }
0139
0140 //
0141 // The tangent numbers grow larger much more rapidly than the Bernoulli numbers do....
0142 // so to compute the Bernoulli numbers from the tangent numbers, we need to avoid spurious
0143 // overflow in the calculation, we can do this by scaling all the tangent number by some scale factor:
0144 //
0145 template <class T, typename std::enable_if<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), bool>::type = true>
0146 inline T tangent_scale_factor()
0147 {
0148    BOOST_MATH_STD_USING
0149    return ldexp(T(1), std::numeric_limits<T>::min_exponent + 5);
0150 }
0151
0152 template <class T, typename std::enable_if<!std::numeric_limits<T>::is_specialized || !(std::numeric_limits<T>::radix == 2), bool>::type = true>
0153 inline T tangent_scale_factor()
0154 {
0155    return tools::min_value<T>() * 16;
0156 }
0157
0158 //
0159 // We need something to act as a cache for our calculated Bernoulli numbers.  In order to
0160 // ensure both fast access and thread safety, we need a stable table which may be extended
0161 // in size, but which never reallocates: that way values already calculated may be accessed
0162 // concurrently with another thread extending the table with new values.
0163 //
0164 // Very very simple vector class that will never allocate more than once, we could use
0165 // boost::container::static_vector here, but that allocates on the stack, which may well
0166 // cause issues for the amount of memory we want in the extreme case...
0167 //
0168 template <class T>
0169 struct fixed_vector : private std::allocator<T>
0170 {
0171    using size_type = unsigned;
0172    using iterator = T*;
0173    using const_iterator = const T*;
0174    fixed_vector() : m_used(0)
0175    {
0176       std::size_t overflow_limit = 5 + b2n_overflow_limit<T, policies::policy<> >();
0177       m_capacity = static_cast<unsigned>((std::min)(overflow_limit, static_cast<std::size_t>(100000u)));
0178       m_data = this->allocate(m_capacity);
0179    }
0180    ~fixed_vector()
0181    {
0182       using allocator_type = std::allocator<T>;
0183       using allocator_traits = std::allocator_traits<allocator_type>;
0184       allocator_type& alloc = *this;
0185       for(unsigned i = 0; i < m_used; ++i)
0186       {
0187          allocator_traits::destroy(alloc, &m_data[i]);
0188       }
0189       allocator_traits::deallocate(alloc, m_data, m_capacity);
0190    }
0191    T& operator[](unsigned n) { BOOST_MATH_ASSERT(n < m_used); return m_data[n]; }
0192    const T& operator[](unsigned n)const { BOOST_MATH_ASSERT(n < m_used); return m_data[n]; }
0193    unsigned size()const { return m_used; }
0194    unsigned size() { return m_used; }
0195    bool resize(unsigned n, const T& val)
0196    {
0197       if(n > m_capacity)
0198       {
0199 #ifndef BOOST_NO_EXCEPTIONS
0200          BOOST_MATH_THROW_EXCEPTION(std::runtime_error("Exhausted storage for Bernoulli numbers."));
0201 #else
0202          return false;
0203 #endif
0204       }
0205       for(unsigned i = m_used; i < n; ++i)
0206          new (m_data + i) T(val);
0207       m_used = n;
0208       return true;
0209    }
0210    bool resize(unsigned n) { return resize(n, T()); }
0211    T* begin() { return m_data; }
0212    T* end() { return m_data + m_used; }
0213    T* begin()const { return m_data; }
0214    T* end()const { return m_data + m_used; }
0215    unsigned capacity()const { return m_capacity; }
0216    void clear() { m_used = 0; }
0217 private:
0218    T* m_data;
0219    unsigned m_used {};
0220    unsigned m_capacity;
0221 };
0222
0223 template <class T, class Policy>
0224 class bernoulli_numbers_cache
0225 {
0226 public:
0227    bernoulli_numbers_cache() : m_overflow_limit((std::numeric_limits<std::size_t>::max)())
0228       , m_counter(0)
0229       , m_current_precision(boost::math::tools::digits<T>())
0230    {}
0231
0232    using container_type = fixed_vector<T>;
0233
0234    bool tangent(std::size_t m)
0235    {
0236       static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1;
0237
0238       if (!tn.resize(static_cast<typename container_type::size_type>(m), T(0U)))
0239       {
0240          return false;
0241       }
0242
0243       BOOST_MATH_INSTRUMENT_VARIABLE(min_overflow_index);
0244
0245       std::size_t prev_size = m_intermediates.size();
0246       m_intermediates.resize(m, T(0U));
0247
0248       if(prev_size == 0)
0249       {
0250          m_intermediates[1] = tangent_scale_factor<T>() /*T(1U)*/;
0251          tn[0U] = T(0U);
0252          tn[1U] = tangent_scale_factor<T>()/* T(1U)*/;
0253          BOOST_MATH_INSTRUMENT_VARIABLE(tn[0]);
0254          BOOST_MATH_INSTRUMENT_VARIABLE(tn[1]);
0255       }
0256
0257       for(std::size_t i = std::max<size_t>(2, prev_size); i < m; i++)
0258       {
0259          bool overflow_check = false;
0260          if(i >= min_overflow_index && (boost::math::tools::max_value<T>() / (i-1) < m_intermediates[1]) )
0261          {
0262             std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>());
0263             break;
0264          }
0265          m_intermediates[1] = m_intermediates[1] * (i-1);
0266          for(std::size_t j = 2; j <= i; j++)
0267          {
0268             overflow_check =
0269                   (i >= min_overflow_index) && (
0270                   (boost::math::tools::max_value<T>() / (i - j) < m_intermediates[j])
0271                   || (boost::math::tools::max_value<T>() / (i - j + 2) < m_intermediates[j-1])
0272                   || (boost::math::tools::max_value<T>() - m_intermediates[j] * (i - j) < m_intermediates[j-1] * (i - j + 2))
0273                   || ((boost::math::isinf)(m_intermediates[j]))
0274                 );
0275
0276             if(overflow_check)
0277             {
0278                std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>());
0279                break;
0280             }
0281             m_intermediates[j] = m_intermediates[j] * (i - j) + m_intermediates[j-1] * (i - j + 2);
0282          }
0283          if(overflow_check)
0284             break; // already filled the tn...
0285          tn[static_cast<typename container_type::size_type>(i)] = m_intermediates[i];
0286          BOOST_MATH_INSTRUMENT_VARIABLE(i);
0287          BOOST_MATH_INSTRUMENT_VARIABLE(tn[static_cast<typename container_type::size_type>(i)]);
0288       }
0289       return true;
0290    }
0291
0292    bool tangent_numbers_series(const std::size_t m)
0293    {
0294       BOOST_MATH_STD_USING
0295       static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1;
0296
0297       typename container_type::size_type old_size = bn.size();
0298
0299       if (!tangent(m))
0300          return false;
0301       if (!bn.resize(static_cast<typename container_type::size_type>(m)))
0302          return false;
0303
0304       if(!old_size)
0305       {
0306          bn[0] = 1;
0307          old_size = 1;
0308       }
0309
0310       T power_two(ldexp(T(1), static_cast<int>(2 * old_size)));
0311
0312       for(std::size_t i = old_size; i < m; i++)
0313       {
0314          T b(static_cast<T>(i * 2));
0315          //
0316          // Not only do we need to take care to avoid spurious over/under flow in
0317          // the calculation, but we also need to avoid overflow altogether in case
0318          // we're calculating with a type where "bad things" happen in that case:
0319          //
0320          b  = b / (power_two * tangent_scale_factor<T>());
0321          b /= (power_two - 1);
0322          bool overflow_check = (i >= min_overflow_index) && (tools::max_value<T>() / tn[static_cast<typename container_type::size_type>(i)] < b);
0323          if(overflow_check)
0324          {
0325             m_overflow_limit = i;
0326             while(i < m)
0327             {
0328                b = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : tools::max_value<T>();
0329                bn[static_cast<typename container_type::size_type>(i)] = ((i % 2U) ? b : T(-b));
0330                ++i;
0331             }
0332             break;
0333          }
0334          else
0335          {
0336             b *= tn[static_cast<typename container_type::size_type>(i)];
0337          }
0338
0339          power_two = ldexp(power_two, 2);
0340
0341          const bool b_neg = i % 2 == 0;
0342
0343          bn[static_cast<typename container_type::size_type>(i)] = ((!b_neg) ? b : T(-b));
0344       }
0345       return true;
0346    }
0347
0348    template <class OutputIterator>
0349    OutputIterator copy_bernoulli_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol)
0350    {
0351       //
0352       // There are basically 3 thread safety options:
0353       //
0354       // 1) There are no threads (BOOST_HAS_THREADS is not defined).
0355       // 2) There are threads, but we do not have a true atomic integer type,
0356       //    in this case we just use a mutex to guard against race conditions.
0357       // 3) There are threads, and we have an atomic integer: in this case we can
0358       //    use the double-checked locking pattern to avoid thread synchronisation
0359       //    when accessing values already in the cache.
0360       //
0361       // First off handle the common case for overflow and/or asymptotic expansion:
0362       //
0363       if(start + n > bn.capacity())
0364       {
0365          if(start < bn.capacity())
0366          {
0367             out = copy_bernoulli_numbers(out, start, bn.capacity() - start, pol);
0368             n -= bn.capacity() - start;
0369             start = static_cast<std::size_t>(bn.capacity());
0370          }
0371          if(start < b2n_overflow_limit<T, Policy>() + 2u)
0372          {
0373             for(; n; ++start, --n)
0374             {
0375                *out = b2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start * 2U));
0376                ++out;
0377             }
0378          }
0379          for(; n; ++start, --n)
0380          {
0381             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(start), pol);
0382             ++out;
0383          }
0384          return out;
0385       }
0386
0387       #if defined(BOOST_HAS_THREADS) && defined(BOOST_MATH_BERNOULLI_NOTHREADS) && !defined(BOOST_MATH_BERNOULLI_UNTHREADED)
0388       // Add a static_assert on instantiation if we have threads, but no C++11 threading support.
0389       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.");
0390       #elif defined(BOOST_MATH_BERNOULLI_NOTHREADS)
0391       //
0392       // Single threaded code, very simple:
0393       //
0394       if(m_current_precision < boost::math::tools::digits<T>())
0395       {
0396          bn.clear();
0397          tn.clear();
0398          m_intermediates.clear();
0399          m_current_precision = boost::math::tools::digits<T>();
0400       }
0401       if(start + n >= bn.size())
0402       {
0403          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()));
0404          if (!tangent_numbers_series(new_size))
0405          {
0406             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));
0407          }
0408       }
0409
0410       for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
0411       {
0412          *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol) : bn[i];
0413          ++out;
0414       }
0415       #else
0416       //
0417       // Double-checked locking pattern, lets us access cached already cached values
0418       // without locking:
0419       //
0420       // Get the counter and see if we need to calculate more constants:
0421       //
0422       if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)
0423          || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))
0424       {
0425          std::lock_guard<std::mutex> l(m_mutex);
0426
0427          if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)
0428             || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))
0429          {
0430             if(static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>())
0431             {
0432                bn.clear();
0433                tn.clear();
0434                m_intermediates.clear();
0435                m_counter.store(0, std::memory_order_release);
0436                m_current_precision = boost::math::tools::digits<T>();
0437             }
0438             if(start + n >= bn.size())
0439             {
0440                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()));
0441                if (!tangent_numbers_series(new_size))
0442                   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));
0443             }
0444             m_counter.store(static_cast<atomic_integer_type>(bn.size()), std::memory_order_release);
0445          }
0446       }
0447
0448       for(std::size_t i = (std::max)(static_cast<std::size_t>(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
0449       {
0450          *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)];
0451          ++out;
0452       }
0453
0454       #endif // BOOST_HAS_THREADS
0455       return out;
0456    }
0457
0458    template <class OutputIterator>
0459    OutputIterator copy_tangent_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol)
0460    {
0461       //
0462       // There are basically 3 thread safety options:
0463       //
0464       // 1) There are no threads (BOOST_HAS_THREADS is not defined).
0465       // 2) There are threads, but we do not have a true atomic integer type,
0466       //    in this case we just use a mutex to guard against race conditions.
0467       // 3) There are threads, and we have an atomic integer: in this case we can
0468       //    use the double-checked locking pattern to avoid thread synchronisation
0469       //    when accessing values already in the cache.
0470       //
0471       //
0472       // First off handle the common case for overflow and/or asymptotic expansion:
0473       //
0474       if(start + n > bn.capacity())
0475       {
0476          if(start < bn.capacity())
0477          {
0478             out = copy_tangent_numbers(out, start, bn.capacity() - start, pol);
0479             n -= bn.capacity() - start;
0480             start = static_cast<std::size_t>(bn.capacity());
0481          }
0482          if(start < b2n_overflow_limit<T, Policy>() + 2u)
0483          {
0484             for(; n; ++start, --n)
0485             {
0486                *out = t2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start));
0487                ++out;
0488             }
0489          }
0490          for(; n; ++start, --n)
0491          {
0492             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(start), pol);
0493             ++out;
0494          }
0495          return out;
0496       }
0497
0498       #if defined(BOOST_MATH_BERNOULLI_NOTHREADS)
0499       //
0500       // Single threaded code, very simple:
0501       //
0502       if(m_current_precision < boost::math::tools::digits<T>())
0503       {
0504          bn.clear();
0505          tn.clear();
0506          m_intermediates.clear();
0507          m_current_precision = boost::math::tools::digits<T>();
0508       }
0509       if(start + n >= bn.size())
0510       {
0511          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()));
0512          if (!tangent_numbers_series(new_size))
0513             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));
0514       }
0515
0516       for(std::size_t i = start; i < start + n; ++i)
0517       {
0518          if(i >= m_overflow_limit)
0519             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol);
0520          else
0521          {
0522             if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
0523                *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol);
0524             else
0525                *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
0526          }
0527          ++out;
0528       }
0529       #elif defined(BOOST_MATH_NO_ATOMIC_INT)
0530       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.");
0531       #else
0532       //
0533       // Double-checked locking pattern, lets us access cached already cached values
0534       // without locking:
0535       //
0536       // Get the counter and see if we need to calculate more constants:
0537       //
0538       if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)
0539          || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))
0540       {
0541          std::lock_guard<std::mutex> l(m_mutex);
0542
0543          if((static_cast<std::size_t>(m_counter.load(std::memory_order_consume)) < start + n)
0544             || (static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>()))
0545          {
0546             if(static_cast<int>(m_current_precision.load(std::memory_order_consume)) < boost::math::tools::digits<T>())
0547             {
0548                bn.clear();
0549                tn.clear();
0550                m_intermediates.clear();
0551                m_counter.store(0, std::memory_order_release);
0552                m_current_precision = boost::math::tools::digits<T>();
0553             }
0554             if(start + n >= bn.size())
0555             {
0556                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()));
0557                if (!tangent_numbers_series(new_size))
0558                   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));
0559             }
0560             m_counter.store(static_cast<atomic_integer_type>(bn.size()), std::memory_order_release);
0561          }
0562       }
0563
0564       for(std::size_t i = start; i < start + n; ++i)
0565       {
0566          if(i >= m_overflow_limit)
0567             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol);
0568          else
0569          {
0570             if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
0571                *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", nullptr, T(i), pol);
0572             else
0573                *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
0574          }
0575          ++out;
0576       }
0577
0578       #endif // BOOST_HAS_THREADS
0579       return out;
0580    }
0581
0582 private:
0583    //
0584    // The caches for Bernoulli and tangent numbers, once allocated,
0585    // these must NEVER EVER reallocate as it breaks our thread
0586    // safety guarantees:
0587    //
0588    fixed_vector<T> bn, tn;
0589    std::vector<T> m_intermediates;
0590    // The value at which we know overflow has already occurred for the Bn:
0591    std::size_t m_overflow_limit;
0592
0593    #if !defined(BOOST_MATH_BERNOULLI_NOTHREADS)
0594    std::mutex m_mutex;
0595    atomic_counter_type m_counter, m_current_precision;
0596    #else
0597    int m_counter;
0598    int m_current_precision;
0599    #endif // BOOST_HAS_THREADS
0600 };
0601
0602 template <class T, class Policy>
0603 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()
0604 {
0605    //
0606    // When numeric_limits<>::digits is zero, the type has either not specialized numeric_limits at all
0607    // or it's precision can vary at runtime.  So make the cache thread_local so that each thread can
0608    // have it's own precision if required:
0609    //
0610    static
0611 #ifndef BOOST_MATH_NO_THREAD_LOCAL_WITH_NON_TRIVIAL_TYPES
0612       BOOST_MATH_THREAD_LOCAL
0613 #endif
0614       bernoulli_numbers_cache<T, Policy> data;
0615    return data;
0616 }
0617 template <class T, class Policy>
0618 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()
0619 {
0620    //
0621    // Note that we rely on C++11 thread-safe initialization here:
0622    //
0623    static bernoulli_numbers_cache<T, Policy> data;
0624    return data;
0625 }
0626
0627 }}}
0628
0629 #endif // BOOST_MATH_BERNOULLI_DETAIL_HPP