Chris@102: ///////////////////////////////////////////////////////////////////////////////
Chris@102: //  Copyright 2013 John Maddock
Chris@102: //  Distributed under the Boost
Chris@102: //  Software License, Version 1.0. (See accompanying file
Chris@102: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@102: 
Chris@102: #ifndef BOOST_MATH_BERNOULLI_DETAIL_HPP
Chris@102: #define BOOST_MATH_BERNOULLI_DETAIL_HPP
Chris@102: 
Chris@102: #include <boost/config.hpp>
Chris@102: #include <boost/detail/lightweight_mutex.hpp>
Chris@102: #include <boost/utility/enable_if.hpp>
Chris@102: #include <boost/math/tools/toms748_solve.hpp>
Chris@102: 
Chris@102: #ifdef BOOST_HAS_THREADS
Chris@102: 
Chris@102: #ifndef BOOST_NO_CXX11_HDR_ATOMIC
Chris@102: #  include <atomic>
Chris@102: #  define BOOST_MATH_ATOMIC_NS std
Chris@102: #if ATOMIC_INT_LOCK_FREE == 2
Chris@102: typedef std::atomic<int> atomic_counter_type;
Chris@102: typedef int atomic_integer_type;
Chris@102: #elif ATOMIC_SHORT_LOCK_FREE == 2
Chris@102: typedef std::atomic<short> atomic_counter_type;
Chris@102: typedef short atomic_integer_type;
Chris@102: #elif ATOMIC_LONG_LOCK_FREE == 2
Chris@102: typedef std::atomic<long> atomic_counter_type;
Chris@102: typedef long atomic_integer_type;
Chris@102: #elif ATOMIC_LLONG_LOCK_FREE == 2
Chris@102: typedef std::atomic<long long> atomic_counter_type;
Chris@102: typedef long long atomic_integer_type;
Chris@102: #else
Chris@102: #  define BOOST_MATH_NO_ATOMIC_INT
Chris@102: #endif
Chris@102: 
Chris@102: #else // BOOST_NO_CXX11_HDR_ATOMIC
Chris@102: //
Chris@102: // We need Boost.Atomic, but on any platform that supports auto-linking we do
Chris@102: // not need to link against a separate library:
Chris@102: //
Chris@102: #define BOOST_ATOMIC_NO_LIB
Chris@102: #include <boost/atomic.hpp>
Chris@102: #  define BOOST_MATH_ATOMIC_NS boost
Chris@102: 
Chris@102: namespace boost{ namespace math{ namespace detail{
Chris@102: 
Chris@102: //
Chris@102: // We need a type to use as an atomic counter:
Chris@102: //
Chris@102: #if BOOST_ATOMIC_INT_LOCK_FREE == 2
Chris@102: typedef boost::atomic<int> atomic_counter_type;
Chris@102: typedef int atomic_integer_type;
Chris@102: #elif BOOST_ATOMIC_SHORT_LOCK_FREE == 2
Chris@102: typedef boost::atomic<short> atomic_counter_type;
Chris@102: typedef short atomic_integer_type;
Chris@102: #elif BOOST_ATOMIC_LONG_LOCK_FREE == 2
Chris@102: typedef boost::atomic<long> atomic_counter_type;
Chris@102: typedef long atomic_integer_type;
Chris@102: #elif BOOST_ATOMIC_LLONG_LOCK_FREE == 2
Chris@102: typedef boost::atomic<long long> atomic_counter_type;
Chris@102: typedef long long atomic_integer_type;
Chris@102: #else
Chris@102: #  define BOOST_MATH_NO_ATOMIC_INT
Chris@102: #endif
Chris@102: 
Chris@102: }}} // namespaces
Chris@102: 
Chris@102: #endif  // BOOST_NO_CXX11_HDR_ATOMIC
Chris@102: 
Chris@102: #endif // BOOST_HAS_THREADS
Chris@102: 
Chris@102: namespace boost{ namespace math{ namespace detail{
Chris@102: //
Chris@102: // Asymptotic expansion for B2n due to
Chris@102: // Luschny LogB3 formula (http://www.luschny.de/math/primes/bernincl.html)
Chris@102: //
Chris@102: template <class T, class Policy>
Chris@102: T b2n_asymptotic(int n)
Chris@102: {
Chris@102:    BOOST_MATH_STD_USING
Chris@102:    const T nx = static_cast<T>(n);
Chris@102:    const T nx2(nx * nx);
Chris@102: 
Chris@102:    const T approximate_log_of_bernoulli_bn = 
Chris@102:         ((boost::math::constants::half<T>() + nx) * log(nx))
Chris@102:         + ((boost::math::constants::half<T>() - nx) * log(boost::math::constants::pi<T>()))
Chris@102:         + (((T(3) / 2) - nx) * boost::math::constants::ln_two<T>())
Chris@102:         + ((nx * (T(2) - (nx2 * 7) * (1 + ((nx2 * 30) * ((nx2 * 12) - 1))))) / (((nx2 * nx2) * nx2) * 2520));
Chris@102:    return ((n / 2) & 1 ? 1 : -1) * (approximate_log_of_bernoulli_bn > tools::log_max_value<T>() 
Chris@102:       ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, nx, Policy())
Chris@102:       : static_cast<T>(exp(approximate_log_of_bernoulli_bn)));
Chris@102: }
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: T t2n_asymptotic(int n)
Chris@102: {
Chris@102:    BOOST_MATH_STD_USING
Chris@102:    // Just get B2n and convert to a Tangent number:
Chris@102:    T t2n = fabs(b2n_asymptotic<T, Policy>(2 * n)) / (2 * n);
Chris@102:    T p2 = ldexp(T(1), n);
Chris@102:    if(tools::max_value<T>() / p2 < t2n)
Chris@102:       return policies::raise_overflow_error<T>("boost::math::tangent_t2n<%1%>(std::size_t)", 0, T(n), Policy());
Chris@102:    t2n *= p2;
Chris@102:    p2 -= 1;
Chris@102:    if(tools::max_value<T>() / p2 < t2n)
Chris@102:       return policies::raise_overflow_error<T>("boost::math::tangent_t2n<%1%>(std::size_t)", 0, Policy());
Chris@102:    t2n *= p2;
Chris@102:    return t2n;
Chris@102: }
Chris@102: //
Chris@102: // We need to know the approximate value of /n/ which will
Chris@102: // cause bernoulli_b2n<T>(n) to return infinity - this allows
Chris@102: // us to elude a great deal of runtime checking for values below
Chris@102: // n, and only perform the full overflow checks when we know that we're
Chris@102: // getting close to the point where our calculations will overflow.
Chris@102: // We use Luschny's LogB3 formula (http://www.luschny.de/math/primes/bernincl.html) 
Chris@102: // to find the limit, and since we're dealing with the log of the Bernoulli numbers
Chris@102: // we need only perform the calculation at double precision and not with T
Chris@102: // (which may be a multiprecision type).  The limit returned is within 1 of the true
Chris@102: // limit for all the types tested.  Note that although the code below is basically
Chris@102: // the same as b2n_asymptotic above, it has been recast as a continuous real-valued 
Chris@102: // function as this makes the root finding go smoother/faster.  It also omits the
Chris@102: // sign of the Bernoulli number.
Chris@102: //
Chris@102: struct max_bernoulli_root_functor
Chris@102: {
Chris@102:    max_bernoulli_root_functor(long long t) : target(static_cast<double>(t)) {}
Chris@102:    double operator()(double n)
Chris@102:    {
Chris@102:       BOOST_MATH_STD_USING
Chris@102: 
Chris@102:       // Luschny LogB3(n) formula.
Chris@102: 
Chris@102:       const double nx2(n * n);
Chris@102: 
Chris@102:       const double approximate_log_of_bernoulli_bn
Chris@102:          =   ((boost::math::constants::half<double>() + n) * log(n))
Chris@102:            + ((boost::math::constants::half<double>() - n) * log(boost::math::constants::pi<double>()))
Chris@102:            + (((double(3) / 2) - n) * boost::math::constants::ln_two<double>())
Chris@102:            + ((n * (2 - (nx2 * 7) * (1 + ((nx2 * 30) * ((nx2 * 12) - 1))))) / (((nx2 * nx2) * nx2) * 2520));
Chris@102: 
Chris@102:       return approximate_log_of_bernoulli_bn - target;
Chris@102:    }
Chris@102: private:
Chris@102:    double target;
Chris@102: };
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: inline std::size_t find_bernoulli_overflow_limit(const mpl::false_&)
Chris@102: {
Chris@102:    long long t = lltrunc(boost::math::tools::log_max_value<T>());
Chris@102:    max_bernoulli_root_functor fun(t);
Chris@102:    boost::math::tools::equal_floor tol;
Chris@102:    boost::uintmax_t max_iter = boost::math::policies::get_max_root_iterations<Policy>();
Chris@102:    return static_cast<std::size_t>(boost::math::tools::toms748_solve(fun, sqrt(double(t)), double(t), tol, max_iter).first) / 2;
Chris@102: }
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: inline std::size_t find_bernoulli_overflow_limit(const mpl::true_&)
Chris@102: {
Chris@102:    return max_bernoulli_index<bernoulli_imp_variant<T>::value>::value;
Chris@102: }
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: std::size_t b2n_overflow_limit()
Chris@102: {
Chris@102:    // This routine is called at program startup if it's called at all:
Chris@102:    // that guarantees safe initialization of the static variable.
Chris@102:    typedef mpl::bool_<(bernoulli_imp_variant<T>::value >= 1) && (bernoulli_imp_variant<T>::value <= 3)> tag_type;
Chris@102:    static const std::size_t lim = find_bernoulli_overflow_limit<T, Policy>(tag_type());
Chris@102:    return lim;
Chris@102: }
Chris@102: 
Chris@102: //
Chris@102: // The tangent numbers grow larger much more rapidly than the Bernoulli numbers do....
Chris@102: // so to compute the Bernoulli numbers from the tangent numbers, we need to avoid spurious
Chris@102: // overflow in the calculation, we can do this by scaling all the tangent number by some scale factor:
Chris@102: //
Chris@102: template <class T>
Chris@102: inline typename enable_if_c<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), T>::type tangent_scale_factor()
Chris@102: {
Chris@102:    BOOST_MATH_STD_USING
Chris@102:    return ldexp(T(1), std::numeric_limits<T>::min_exponent + 5);
Chris@102: }
Chris@102: template <class T>
Chris@102: inline typename disable_if_c<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), T>::type tangent_scale_factor()
Chris@102: {
Chris@102:    return tools::min_value<T>() * 16;
Chris@102: }
Chris@102: //
Chris@102: // Initializer: ensure all our constants are initialized prior to the first call of main:
Chris@102: //
Chris@102: template <class T, class Policy>
Chris@102: struct bernoulli_initializer
Chris@102: {
Chris@102:    struct init
Chris@102:    {
Chris@102:       init()
Chris@102:       {
Chris@102:          //
Chris@102:          // We call twice, once to initialize our static table, and once to
Chris@102:          // initialize our dymanic table:
Chris@102:          //
Chris@102:          boost::math::bernoulli_b2n<T>(2, Policy());
Chris@102:          try{
Chris@102:             boost::math::bernoulli_b2n<T>(max_bernoulli_b2n<T>::value + 1, Policy());
Chris@102:          } catch(const std::overflow_error&){}
Chris@102:          boost::math::tangent_t2n<T>(2, Policy());
Chris@102:       }
Chris@102:       void force_instantiate()const{}
Chris@102:    };
Chris@102:    static const init initializer;
Chris@102:    static void force_instantiate()
Chris@102:    {
Chris@102:       initializer.force_instantiate();
Chris@102:    }
Chris@102: };
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: const typename bernoulli_initializer<T, Policy>::init bernoulli_initializer<T, Policy>::initializer;
Chris@102: 
Chris@102: //
Chris@102: // We need something to act as a cache for our calculated Bernoulli numbers.  In order to
Chris@102: // ensure both fast access and thread safety, we need a stable table which may be extended
Chris@102: // in size, but which never reallocates: that way values already calculated may be accessed
Chris@102: // concurrently with another thread extending the table with new values.
Chris@102: //
Chris@102: // Very very simple vector class that will never allocate more than once, we could use
Chris@102: // boost::container::static_vector here, but that allocates on the stack, which may well
Chris@102: // cause issues for the amount of memory we want in the extreme case...
Chris@102: //
Chris@102: template <class T>
Chris@102: struct fixed_vector : private std::allocator<T>
Chris@102: {
Chris@102:    typedef unsigned size_type;
Chris@102:    typedef T* iterator;
Chris@102:    typedef const T* const_iterator;
Chris@102:    fixed_vector() : m_used(0)
Chris@102:    { 
Chris@102:       std::size_t overflow_limit = 5 + b2n_overflow_limit<T, policies::policy<> >();
Chris@102:       m_capacity = static_cast<unsigned>((std::min)(overflow_limit, static_cast<std::size_t>(100000u)));
Chris@102:       m_data = this->allocate(m_capacity); 
Chris@102:    }
Chris@102:    ~fixed_vector()
Chris@102:    {
Chris@102:       for(unsigned i = 0; i < m_used; ++i)
Chris@102:          this->destroy(&m_data[i]);
Chris@102:       this->deallocate(m_data, m_capacity);
Chris@102:    }
Chris@102:    T& operator[](unsigned n) { BOOST_ASSERT(n < m_used); return m_data[n]; }
Chris@102:    const T& operator[](unsigned n)const { BOOST_ASSERT(n < m_used); return m_data[n]; }
Chris@102:    unsigned size()const { return m_used; }
Chris@102:    unsigned size() { return m_used; }
Chris@102:    void resize(unsigned n, const T& val)
Chris@102:    {
Chris@102:       if(n > m_capacity)
Chris@102:          throw std::runtime_error("Exhausted storage for Bernoulli numbers.");
Chris@102:       for(unsigned i = m_used; i < n; ++i)
Chris@102:          new (m_data + i) T(val);
Chris@102:       m_used = n;
Chris@102:    }
Chris@102:    void resize(unsigned n) { resize(n, T()); }
Chris@102:    T* begin() { return m_data; }
Chris@102:    T* end() { return m_data + m_used; }
Chris@102:    T* begin()const { return m_data; }
Chris@102:    T* end()const { return m_data + m_used; }
Chris@102:    unsigned capacity()const { return m_capacity; }
Chris@102: private:
Chris@102:    T* m_data;
Chris@102:    unsigned m_used, m_capacity;
Chris@102: };
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: class bernoulli_numbers_cache
Chris@102: {
Chris@102: public:
Chris@102:    bernoulli_numbers_cache() : m_overflow_limit((std::numeric_limits<std::size_t>::max)())
Chris@102: #if defined(BOOST_HAS_THREADS) && !defined(BOOST_MATH_NO_ATOMIC_INT)
Chris@102:       , m_counter(0)
Chris@102: #endif
Chris@102:    {}
Chris@102: 
Chris@102:    typedef fixed_vector<T> container_type;
Chris@102: 
Chris@102:    void tangent(std::size_t m)
Chris@102:    {
Chris@102:       static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1;
Chris@102:       tn.resize(static_cast<typename container_type::size_type>(m), T(0U));
Chris@102: 
Chris@102:       BOOST_MATH_INSTRUMENT_VARIABLE(min_overflow_index);
Chris@102: 
Chris@102:       std::size_t prev_size = m_intermediates.size();
Chris@102:       m_intermediates.resize(m, T(0U));
Chris@102: 
Chris@102:       if(prev_size == 0)
Chris@102:       {
Chris@102:          m_intermediates[1] = tangent_scale_factor<T>() /*T(1U)*/;
Chris@102:          tn[0U] = T(0U);
Chris@102:          tn[1U] = tangent_scale_factor<T>()/* T(1U)*/;
Chris@102:          BOOST_MATH_INSTRUMENT_VARIABLE(tn[0]);
Chris@102:          BOOST_MATH_INSTRUMENT_VARIABLE(tn[1]);
Chris@102:       }
Chris@102: 
Chris@102:       for(std::size_t i = std::max<size_t>(2, prev_size); i < m; i++)
Chris@102:       {
Chris@102:          bool overflow_check = false;
Chris@102:          if(i >= min_overflow_index && (boost::math::tools::max_value<T>() / (i-1) < m_intermediates[1]) )
Chris@102:          {
Chris@102:             std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>());
Chris@102:             break;
Chris@102:          }
Chris@102:          m_intermediates[1] = m_intermediates[1] * (i-1);
Chris@102:          for(std::size_t j = 2; j <= i; j++)
Chris@102:          {
Chris@102:             overflow_check =
Chris@102:                   (i >= min_overflow_index) && (
Chris@102:                   (boost::math::tools::max_value<T>() / (i - j) < m_intermediates[j])
Chris@102:                   || (boost::math::tools::max_value<T>() / (i - j + 2) < m_intermediates[j-1])
Chris@102:                   || (boost::math::tools::max_value<T>() - m_intermediates[j] * (i - j) < m_intermediates[j-1] * (i - j + 2))
Chris@102:                   || ((boost::math::isinf)(m_intermediates[j]))
Chris@102:                 );
Chris@102: 
Chris@102:             if(overflow_check)
Chris@102:             {
Chris@102:                std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>());
Chris@102:                break;
Chris@102:             }
Chris@102:             m_intermediates[j] = m_intermediates[j] * (i - j) + m_intermediates[j-1] * (i - j + 2);
Chris@102:          }
Chris@102:          if(overflow_check)
Chris@102:             break; // already filled the tn...
Chris@102:          tn[static_cast<typename container_type::size_type>(i)] = m_intermediates[i];
Chris@102:          BOOST_MATH_INSTRUMENT_VARIABLE(i);
Chris@102:          BOOST_MATH_INSTRUMENT_VARIABLE(tn[static_cast<typename container_type::size_type>(i)]);
Chris@102:       }
Chris@102:    }
Chris@102: 
Chris@102:    void tangent_numbers_series(const std::size_t m)
Chris@102:    {
Chris@102:       BOOST_MATH_STD_USING
Chris@102:       static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1;
Chris@102: 
Chris@102:       typename container_type::size_type old_size = bn.size();
Chris@102: 
Chris@102:       tangent(m);
Chris@102:       bn.resize(static_cast<typename container_type::size_type>(m));
Chris@102: 
Chris@102:       if(!old_size)
Chris@102:       {
Chris@102:          bn[0] = 1;
Chris@102:          old_size = 1;
Chris@102:       }
Chris@102: 
Chris@102:       T power_two(ldexp(T(1), static_cast<int>(2 * old_size)));
Chris@102: 
Chris@102:       for(std::size_t i = old_size; i < m; i++)
Chris@102:       {
Chris@102:          T b(static_cast<T>(i * 2));
Chris@102:          //
Chris@102:          // Not only do we need to take care to avoid spurious over/under flow in
Chris@102:          // the calculation, but we also need to avoid overflow altogether in case
Chris@102:          // we're calculating with a type where "bad things" happen in that case:
Chris@102:          //
Chris@102:          b  = b / (power_two * tangent_scale_factor<T>());
Chris@102:          b /= (power_two - 1);
Chris@102:          bool overflow_check = (i >= min_overflow_index) && (tools::max_value<T>() / tn[static_cast<typename container_type::size_type>(i)] < b);
Chris@102:          if(overflow_check)
Chris@102:          {
Chris@102:             m_overflow_limit = i;
Chris@102:             while(i < m)
Chris@102:             {
Chris@102:                b = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : tools::max_value<T>();
Chris@102:                bn[static_cast<typename container_type::size_type>(i)] = ((i % 2U) ? b : T(-b));
Chris@102:                ++i;
Chris@102:             }
Chris@102:             break;
Chris@102:          }
Chris@102:          else
Chris@102:          {
Chris@102:             b *= tn[static_cast<typename container_type::size_type>(i)];
Chris@102:          }
Chris@102: 
Chris@102:          power_two = ldexp(power_two, 2);
Chris@102: 
Chris@102:          const bool b_neg = i % 2 == 0;
Chris@102: 
Chris@102:          bn[static_cast<typename container_type::size_type>(i)] = ((!b_neg) ? b : T(-b));
Chris@102:       }
Chris@102:    }
Chris@102: 
Chris@102:    template <class OutputIterator>
Chris@102:    OutputIterator copy_bernoulli_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol)
Chris@102:    {
Chris@102:       //
Chris@102:       // There are basically 3 thread safety options:
Chris@102:       //
Chris@102:       // 1) There are no threads (BOOST_HAS_THREADS is not defined).
Chris@102:       // 2) There are threads, but we do not have a true atomic integer type, 
Chris@102:       //    in this case we just use a mutex to guard against race conditions.
Chris@102:       // 3) There are threads, and we have an atomic integer: in this case we can
Chris@102:       //    use the double-checked locking pattern to avoid thread synchronisation
Chris@102:       //    when accessing values already in the cache.
Chris@102:       //
Chris@102:       // First off handle the common case for overflow and/or asymptotic expansion:
Chris@102:       //
Chris@102:       if(start + n > bn.capacity())
Chris@102:       {
Chris@102:          if(start < bn.capacity())
Chris@102:          {
Chris@102:             out = copy_bernoulli_numbers(out, start, bn.capacity() - start, pol);
Chris@102:             n -= bn.capacity() - start;
Chris@102:             start = static_cast<std::size_t>(bn.capacity());
Chris@102:          }
Chris@102:          if(start < b2n_overflow_limit<T, Policy>() + 2u)
Chris@102:          {
Chris@102:             for(; n; ++start, --n)
Chris@102:             {
Chris@102:                *out = b2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start * 2U));
Chris@102:                ++out;
Chris@102:             }
Chris@102:          }
Chris@102:          for(; n; ++start, --n)
Chris@102:          {
Chris@102:             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(start), pol);
Chris@102:             ++out;
Chris@102:          }
Chris@102:          return out;
Chris@102:       }
Chris@102:    #if !defined(BOOST_HAS_THREADS)
Chris@102:       //
Chris@102:       // Single threaded code, very simple:
Chris@102:       //
Chris@102:       if(start + n >= bn.size())
Chris@102:       {
Chris@102:          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()));
Chris@102:          tangent_numbers_series(new_size);
Chris@102:       }
Chris@102: 
Chris@102:       for(std::size_t i = (std::max)(max_bernoulli_b2n<T>::value + 1, start); i < start + n; ++i)
Chris@102:       {
Chris@102:          *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
Chris@102:          ++out;
Chris@102:       }
Chris@102:    #elif defined(BOOST_MATH_NO_ATOMIC_INT)
Chris@102:       //
Chris@102:       // We need to grab a mutex every time we get here, for both readers and writers:
Chris@102:       //
Chris@102:       boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
Chris@102:       if(start + n >= bn.size())
Chris@102:       {
Chris@102:          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()));
Chris@102:          tangent_numbers_series(new_size);
Chris@102:       }
Chris@102: 
Chris@102:       for(std::size_t i = (std::max)(max_bernoulli_b2n<T>::value + 1, start); i < start + n; ++i)
Chris@102:       {
Chris@102:          *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
Chris@102:          ++out;
Chris@102:       }
Chris@102: 
Chris@102:    #else
Chris@102:       //
Chris@102:       // Double-checked locking pattern, lets us access cached already cached values
Chris@102:       // without locking:
Chris@102:       //
Chris@102:       // Get the counter and see if we need to calculate more constants:
Chris@102:       //
Chris@102:       if(static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
Chris@102:       {
Chris@102:          boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
Chris@102: 
Chris@102:          if(static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
Chris@102:          {
Chris@102:             if(start + n >= bn.size())
Chris@102:             {
Chris@102:                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()));
Chris@102:                tangent_numbers_series(new_size);
Chris@102:             }
Chris@102:             m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release);
Chris@102:          }
Chris@102:       }
Chris@102: 
Chris@102:       for(std::size_t i = (std::max)(static_cast<std::size_t>(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
Chris@102:       {
Chris@102:          *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[static_cast<typename container_type::size_type>(i)];
Chris@102:          ++out;
Chris@102:       }
Chris@102: 
Chris@102:    #endif
Chris@102:       return out;
Chris@102:    }
Chris@102: 
Chris@102:    template <class OutputIterator>
Chris@102:    OutputIterator copy_tangent_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol)
Chris@102:    {
Chris@102:       //
Chris@102:       // There are basically 3 thread safety options:
Chris@102:       //
Chris@102:       // 1) There are no threads (BOOST_HAS_THREADS is not defined).
Chris@102:       // 2) There are threads, but we do not have a true atomic integer type, 
Chris@102:       //    in this case we just use a mutex to guard against race conditions.
Chris@102:       // 3) There are threads, and we have an atomic integer: in this case we can
Chris@102:       //    use the double-checked locking pattern to avoid thread synchronisation
Chris@102:       //    when accessing values already in the cache.
Chris@102:       //
Chris@102:       //
Chris@102:       // First off handle the common case for overflow and/or asymptotic expansion:
Chris@102:       //
Chris@102:       if(start + n > bn.capacity())
Chris@102:       {
Chris@102:          if(start < bn.capacity())
Chris@102:          {
Chris@102:             out = copy_tangent_numbers(out, start, bn.capacity() - start, pol);
Chris@102:             n -= bn.capacity() - start;
Chris@102:             start = static_cast<std::size_t>(bn.capacity());
Chris@102:          }
Chris@102:          if(start < b2n_overflow_limit<T, Policy>() + 2u)
Chris@102:          {
Chris@102:             for(; n; ++start, --n)
Chris@102:             {
Chris@102:                *out = t2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start));
Chris@102:                ++out;
Chris@102:             }
Chris@102:          }
Chris@102:          for(; n; ++start, --n)
Chris@102:          {
Chris@102:             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(start), pol);
Chris@102:             ++out;
Chris@102:          }
Chris@102:          return out;
Chris@102:       }
Chris@102:    #if !defined(BOOST_HAS_THREADS)
Chris@102:       //
Chris@102:       // Single threaded code, very simple:
Chris@102:       //
Chris@102:       if(start + n >= bn.size())
Chris@102:       {
Chris@102:          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()));
Chris@102:          tangent_numbers_series(new_size);
Chris@102:       }
Chris@102: 
Chris@102:       for(std::size_t i = start; i < start + n; ++i)
Chris@102:       {
Chris@102:          if(i >= m_overflow_limit)
Chris@102:             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
Chris@102:          else
Chris@102:          {
Chris@102:             if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
Chris@102:                *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
Chris@102:             else
Chris@102:                *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
Chris@102:          }
Chris@102:          ++out;
Chris@102:       }
Chris@102:    #elif defined(BOOST_MATH_NO_ATOMIC_INT)
Chris@102:       //
Chris@102:       // We need to grab a mutex every time we get here, for both readers and writers:
Chris@102:       //
Chris@102:       boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
Chris@102:       if(start + n >= bn.size())
Chris@102:       {
Chris@102:          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()));
Chris@102:          tangent_numbers_series(new_size);
Chris@102:       }
Chris@102: 
Chris@102:       for(std::size_t i = start; i < start + n; ++i)
Chris@102:       {
Chris@102:          if(i >= m_overflow_limit)
Chris@102:             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
Chris@102:          else
Chris@102:          {
Chris@102:             if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
Chris@102:                *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
Chris@102:             else
Chris@102:                *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
Chris@102:          }
Chris@102:          ++out;
Chris@102:       }
Chris@102: 
Chris@102:    #else
Chris@102:       //
Chris@102:       // Double-checked locking pattern, lets us access cached already cached values
Chris@102:       // without locking:
Chris@102:       //
Chris@102:       // Get the counter and see if we need to calculate more constants:
Chris@102:       //
Chris@102:       if(static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
Chris@102:       {
Chris@102:          boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
Chris@102: 
Chris@102:          if(static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
Chris@102:          {
Chris@102:             if(start + n >= bn.size())
Chris@102:             {
Chris@102:                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()));
Chris@102:                tangent_numbers_series(new_size);
Chris@102:             }
Chris@102:             m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release);
Chris@102:          }
Chris@102:       }
Chris@102: 
Chris@102:       for(std::size_t i = start; i < start + n; ++i)
Chris@102:       {
Chris@102:          if(i >= m_overflow_limit)
Chris@102:             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
Chris@102:          else
Chris@102:          {
Chris@102:             if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
Chris@102:                *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
Chris@102:             else
Chris@102:                *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
Chris@102:          }
Chris@102:          ++out;
Chris@102:       }
Chris@102: 
Chris@102:    #endif
Chris@102:       return out;
Chris@102:    }
Chris@102: 
Chris@102: private:
Chris@102:    //
Chris@102:    // The caches for Bernoulli and tangent numbers, once allocated,
Chris@102:    // these must NEVER EVER reallocate as it breaks our thread
Chris@102:    // safety guarentees:
Chris@102:    //
Chris@102:    fixed_vector<T> bn, tn;
Chris@102:    std::vector<T> m_intermediates;
Chris@102:    // The value at which we know overflow has already occured for the Bn:
Chris@102:    std::size_t m_overflow_limit;
Chris@102: #if !defined(BOOST_HAS_THREADS)
Chris@102: #elif defined(BOOST_MATH_NO_ATOMIC_INT)
Chris@102:    boost::detail::lightweight_mutex m_mutex;
Chris@102: #else
Chris@102:    boost::detail::lightweight_mutex m_mutex;
Chris@102:    atomic_counter_type m_counter;
Chris@102: #endif
Chris@102: };
Chris@102: 
Chris@102: template <class T, class Policy>
Chris@102: inline bernoulli_numbers_cache<T, Policy>& get_bernoulli_numbers_cache()
Chris@102: {
Chris@102:    //
Chris@102:    // Force this function to be called at program startup so all the static variables
Chris@102:    // get initailzed then (thread safety).
Chris@102:    //
Chris@102:    bernoulli_initializer<T, Policy>::force_instantiate();
Chris@102:    static bernoulli_numbers_cache<T, Policy> data;
Chris@102:    return data;
Chris@102: }
Chris@102: 
Chris@102: }}}
Chris@102: 
Chris@102: #endif // BOOST_MATH_BERNOULLI_DETAIL_HPP