boost/math/special_functions/detail/bernoulli_details.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //  Copyright 2013 John Maddock
   3 //  Distributed under the Boost
   4 //  Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_MATH_BERNOULLI_DETAIL_HPP
   8 #define BOOST_MATH_BERNOULLI_DETAIL_HPP
   9
  10 #include <boost/config.hpp>
  11 #include <boost/detail/lightweight_mutex.hpp>
  12 #include <boost/math/tools/atomic.hpp>
  13 #include <boost/utility/enable_if.hpp>
  14 #include <boost/math/tools/toms748_solve.hpp>
  15 #include <vector>
  16
  17 namespace boost{ namespace math{ namespace detail{
  18 //
  19 // Asymptotic expansion for B2n due to
  20 // Luschny LogB3 formula (http://www.luschny.de/math/primes/bernincl.html)
  21 //
  22 template <class T, class Policy>
  23 T b2n_asymptotic(int n)
  24 {
  25    BOOST_MATH_STD_USING
  26    const T nx = static_cast<T>(n);
  27    const T nx2(nx * nx);
  28
  29    const T approximate_log_of_bernoulli_bn =
  30         ((boost::math::constants::half<T>() + nx) * log(nx))
  31         + ((boost::math::constants::half<T>() - nx) * log(boost::math::constants::pi<T>()))
  32         + (((T(3) / 2) - nx) * boost::math::constants::ln_two<T>())
  33         + ((nx * (T(2) - (nx2 * 7) * (1 + ((nx2 * 30) * ((nx2 * 12) - 1))))) / (((nx2 * nx2) * nx2) * 2520));
  34    return ((n / 2) & 1 ? 1 : -1) * (approximate_log_of_bernoulli_bn > tools::log_max_value<T>()
  35       ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, nx, Policy())
  36       : static_cast<T>(exp(approximate_log_of_bernoulli_bn)));
  37 }
  38
  39 template <class T, class Policy>
  40 T t2n_asymptotic(int n)
  41 {
  42    BOOST_MATH_STD_USING
  43    // Just get B2n and convert to a Tangent number:
  44    T t2n = fabs(b2n_asymptotic<T, Policy>(2 * n)) / (2 * n);
  45    T p2 = ldexp(T(1), n);
  46    if(tools::max_value<T>() / p2 < t2n)
  47       return policies::raise_overflow_error<T>("boost::math::tangent_t2n<%1%>(std::size_t)", 0, T(n), Policy());
  48    t2n *= p2;
  49    p2 -= 1;
  50    if(tools::max_value<T>() / p2 < t2n)
  51       return policies::raise_overflow_error<T>("boost::math::tangent_t2n<%1%>(std::size_t)", 0, Policy());
  52    t2n *= p2;
  53    return t2n;
  54 }
  55 //
  56 // We need to know the approximate value of /n/ which will
  57 // cause bernoulli_b2n<T>(n) to return infinity - this allows
  58 // us to elude a great deal of runtime checking for values below
  59 // n, and only perform the full overflow checks when we know that we're
  60 // getting close to the point where our calculations will overflow.
  61 // We use Luschny's LogB3 formula (http://www.luschny.de/math/primes/bernincl.html)
  62 // to find the limit, and since we're dealing with the log of the Bernoulli numbers
  63 // we need only perform the calculation at double precision and not with T
  64 // (which may be a multiprecision type).  The limit returned is within 1 of the true
  65 // limit for all the types tested.  Note that although the code below is basically
  66 // the same as b2n_asymptotic above, it has been recast as a continuous real-valued
  67 // function as this makes the root finding go smoother/faster.  It also omits the
  68 // sign of the Bernoulli number.
  69 //
  70 struct max_bernoulli_root_functor
  71 {
  72    max_bernoulli_root_functor(ulong_long_type t) : target(static_cast<double>(t)) {}
  73    double operator()(double n)
  74    {
  75       BOOST_MATH_STD_USING
  76
  77       // Luschny LogB3(n) formula.
  78
  79       const double nx2(n * n);
  80
  81       const double approximate_log_of_bernoulli_bn
  82          =   ((boost::math::constants::half<double>() + n) * log(n))
  83            + ((boost::math::constants::half<double>() - n) * log(boost::math::constants::pi<double>()))
  84            + (((double(3) / 2) - n) * boost::math::constants::ln_two<double>())
  85            + ((n * (2 - (nx2 * 7) * (1 + ((nx2 * 30) * ((nx2 * 12) - 1))))) / (((nx2 * nx2) * nx2) * 2520));
  86
  87       return approximate_log_of_bernoulli_bn - target;
  88    }
  89 private:
  90    double target;
  91 };
  92
  93 template <class T, class Policy>
  94 inline std::size_t find_bernoulli_overflow_limit(const mpl::false_&)
  95 {
  96    // Set a limit on how large the result can ever be:
  97    static const double max_result = static_cast<double>((std::numeric_limits<std::size_t>::max)() - 1000u);
  98
  99    ulong_long_type t = lltrunc(boost::math::tools::log_max_value<T>());
 100    max_bernoulli_root_functor fun(t);
 101    boost::math::tools::equal_floor tol;
 102    boost::uintmax_t max_iter = boost::math::policies::get_max_root_iterations<Policy>();
 103    double result = boost::math::tools::toms748_solve(fun, sqrt(double(t)), double(t), tol, max_iter).first / 2;
 104    if (result > max_result)
 105       result = max_result;
 106
 107    return static_cast<std::size_t>(result);
 108 }
 109
 110 template <class T, class Policy>
 111 inline std::size_t find_bernoulli_overflow_limit(const mpl::true_&)
 112 {
 113    return max_bernoulli_index<bernoulli_imp_variant<T>::value>::value;
 114 }
 115
 116 template <class T, class Policy>
 117 std::size_t b2n_overflow_limit()
 118 {
 119    // This routine is called at program startup if it's called at all:
 120    // that guarantees safe initialization of the static variable.
 121    typedef mpl::bool_<(bernoulli_imp_variant<T>::value >= 1) && (bernoulli_imp_variant<T>::value <= 3)> tag_type;
 122    static const std::size_t lim = find_bernoulli_overflow_limit<T, Policy>(tag_type());
 123    return lim;
 124 }
 125
 126 //
 127 // The tangent numbers grow larger much more rapidly than the Bernoulli numbers do....
 128 // so to compute the Bernoulli numbers from the tangent numbers, we need to avoid spurious
 129 // overflow in the calculation, we can do this by scaling all the tangent number by some scale factor:
 130 //
 131 template <class T>
 132 inline typename enable_if_c<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), T>::type tangent_scale_factor()
 133 {
 134    BOOST_MATH_STD_USING
 135    return ldexp(T(1), std::numeric_limits<T>::min_exponent + 5);
 136 }
 137 template <class T>
 138 inline typename disable_if_c<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), T>::type tangent_scale_factor()
 139 {
 140    return tools::min_value<T>() * 16;
 141 }
 142 //
 143 // Initializer: ensure all our constants are initialized prior to the first call of main:
 144 //
 145 template <class T, class Policy>
 146 struct bernoulli_initializer
 147 {
 148    struct init
 149    {
 150       init()
 151       {
 152          //
 153          // We call twice, once to initialize our static table, and once to
 154          // initialize our dymanic table:
 155          //
 156          boost::math::bernoulli_b2n<T>(2, Policy());
 157 #ifndef BOOST_NO_EXCEPTIONS
 158          try{
 159 #endif
 160             boost::math::bernoulli_b2n<T>(max_bernoulli_b2n<T>::value + 1, Policy());
 161 #ifndef BOOST_NO_EXCEPTIONS
 162          } catch(const std::overflow_error&){}
 163 #endif
 164          boost::math::tangent_t2n<T>(2, Policy());
 165       }
 166       void force_instantiate()const{}
 167    };
 168    static const init initializer;
 169    static void force_instantiate()
 170    {
 171       initializer.force_instantiate();
 172    }
 173 };
 174
 175 template <class T, class Policy>
 176 const typename bernoulli_initializer<T, Policy>::init bernoulli_initializer<T, Policy>::initializer;
 177
 178 //
 179 // We need something to act as a cache for our calculated Bernoulli numbers.  In order to
 180 // ensure both fast access and thread safety, we need a stable table which may be extended
 181 // in size, but which never reallocates: that way values already calculated may be accessed
 182 // concurrently with another thread extending the table with new values.
 183 //
 184 // Very very simple vector class that will never allocate more than once, we could use
 185 // boost::container::static_vector here, but that allocates on the stack, which may well
 186 // cause issues for the amount of memory we want in the extreme case...
 187 //
 188 template <class T>
 189 struct fixed_vector : private std::allocator<T>
 190 {
 191    typedef unsigned size_type;
 192    typedef T* iterator;
 193    typedef const T* const_iterator;
 194    fixed_vector() : m_used(0)
 195    {
 196       std::size_t overflow_limit = 5 + b2n_overflow_limit<T, policies::policy<> >();
 197       m_capacity = static_cast<unsigned>((std::min)(overflow_limit, static_cast<std::size_t>(100000u)));
 198       m_data = this->allocate(m_capacity);
 199    }
 200    ~fixed_vector()
 201    {
 202 #ifdef BOOST_NO_CXX11_ALLOCATOR
 203       for(unsigned i = 0; i < m_used; ++i)
 204          this->destroy(&m_data[i]);
 205       this->deallocate(m_data, m_capacity);
 206 #else
 207       typedef std::allocator<T> allocator_type;
 208       typedef std::allocator_traits<allocator_type> allocator_traits;
 209       allocator_type& alloc = *this;
 210       for(unsigned i = 0; i < m_used; ++i)
 211          allocator_traits::destroy(alloc, &m_data[i]);
 212       allocator_traits::deallocate(alloc, m_data, m_capacity);
 213 #endif
 214    }
 215    T& operator[](unsigned n) { BOOST_ASSERT(n < m_used); return m_data[n]; }
 216    const T& operator[](unsigned n)const { BOOST_ASSERT(n < m_used); return m_data[n]; }
 217    unsigned size()const { return m_used; }
 218    unsigned size() { return m_used; }
 219    void resize(unsigned n, const T& val)
 220    {
 221       if(n > m_capacity)
 222       {
 223          BOOST_THROW_EXCEPTION(std::runtime_error("Exhausted storage for Bernoulli numbers."));
 224       }
 225       for(unsigned i = m_used; i < n; ++i)
 226          new (m_data + i) T(val);
 227       m_used = n;
 228    }
 229    void resize(unsigned n) { resize(n, T()); }
 230    T* begin() { return m_data; }
 231    T* end() { return m_data + m_used; }
 232    T* begin()const { return m_data; }
 233    T* end()const { return m_data + m_used; }
 234    unsigned capacity()const { return m_capacity; }
 235    void clear() { m_used = 0; }
 236 private:
 237    T* m_data;
 238    unsigned m_used, m_capacity;
 239 };
 240
 241 template <class T, class Policy>
 242 class bernoulli_numbers_cache
 243 {
 244 public:
 245    bernoulli_numbers_cache() : m_overflow_limit((std::numeric_limits<std::size_t>::max)())
 246 #if defined(BOOST_HAS_THREADS) && !defined(BOOST_MATH_NO_ATOMIC_INT)
 247       , m_counter(0)
 248 #endif
 249       , m_current_precision(boost::math::tools::digits<T>())
 250    {}
 251
 252    typedef fixed_vector<T> container_type;
 253
 254    void tangent(std::size_t m)
 255    {
 256       static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1;
 257       tn.resize(static_cast<typename container_type::size_type>(m), T(0U));
 258
 259       BOOST_MATH_INSTRUMENT_VARIABLE(min_overflow_index);
 260
 261       std::size_t prev_size = m_intermediates.size();
 262       m_intermediates.resize(m, T(0U));
 263
 264       if(prev_size == 0)
 265       {
 266          m_intermediates[1] = tangent_scale_factor<T>() /*T(1U)*/;
 267          tn[0U] = T(0U);
 268          tn[1U] = tangent_scale_factor<T>()/* T(1U)*/;
 269          BOOST_MATH_INSTRUMENT_VARIABLE(tn[0]);
 270          BOOST_MATH_INSTRUMENT_VARIABLE(tn[1]);
 271       }
 272
 273       for(std::size_t i = std::max<size_t>(2, prev_size); i < m; i++)
 274       {
 275          bool overflow_check = false;
 276          if(i >= min_overflow_index && (boost::math::tools::max_value<T>() / (i-1) < m_intermediates[1]) )
 277          {
 278             std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>());
 279             break;
 280          }
 281          m_intermediates[1] = m_intermediates[1] * (i-1);
 282          for(std::size_t j = 2; j <= i; j++)
 283          {
 284             overflow_check =
 285                   (i >= min_overflow_index) && (
 286                   (boost::math::tools::max_value<T>() / (i - j) < m_intermediates[j])
 287                   || (boost::math::tools::max_value<T>() / (i - j + 2) < m_intermediates[j-1])
 288                   || (boost::math::tools::max_value<T>() - m_intermediates[j] * (i - j) < m_intermediates[j-1] * (i - j + 2))
 289                   || ((boost::math::isinf)(m_intermediates[j]))
 290                 );
 291
 292             if(overflow_check)
 293             {
 294                std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>());
 295                break;
 296             }
 297             m_intermediates[j] = m_intermediates[j] * (i - j) + m_intermediates[j-1] * (i - j + 2);
 298          }
 299          if(overflow_check)
 300             break; // already filled the tn...
 301          tn[static_cast<typename container_type::size_type>(i)] = m_intermediates[i];
 302          BOOST_MATH_INSTRUMENT_VARIABLE(i);
 303          BOOST_MATH_INSTRUMENT_VARIABLE(tn[static_cast<typename container_type::size_type>(i)]);
 304       }
 305    }
 306
 307    void tangent_numbers_series(const std::size_t m)
 308    {
 309       BOOST_MATH_STD_USING
 310       static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1;
 311
 312       typename container_type::size_type old_size = bn.size();
 313
 314       tangent(m);
 315       bn.resize(static_cast<typename container_type::size_type>(m));
 316
 317       if(!old_size)
 318       {
 319          bn[0] = 1;
 320          old_size = 1;
 321       }
 322
 323       T power_two(ldexp(T(1), static_cast<int>(2 * old_size)));
 324
 325       for(std::size_t i = old_size; i < m; i++)
 326       {
 327          T b(static_cast<T>(i * 2));
 328          //
 329          // Not only do we need to take care to avoid spurious over/under flow in
 330          // the calculation, but we also need to avoid overflow altogether in case
 331          // we're calculating with a type where "bad things" happen in that case:
 332          //
 333          b  = b / (power_two * tangent_scale_factor<T>());
 334          b /= (power_two - 1);
 335          bool overflow_check = (i >= min_overflow_index) && (tools::max_value<T>() / tn[static_cast<typename container_type::size_type>(i)] < b);
 336          if(overflow_check)
 337          {
 338             m_overflow_limit = i;
 339             while(i < m)
 340             {
 341                b = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : tools::max_value<T>();
 342                bn[static_cast<typename container_type::size_type>(i)] = ((i % 2U) ? b : T(-b));
 343                ++i;
 344             }
 345             break;
 346          }
 347          else
 348          {
 349             b *= tn[static_cast<typename container_type::size_type>(i)];
 350          }
 351
 352          power_two = ldexp(power_two, 2);
 353
 354          const bool b_neg = i % 2 == 0;
 355
 356          bn[static_cast<typename container_type::size_type>(i)] = ((!b_neg) ? b : T(-b));
 357       }
 358    }
 359
 360    template <class OutputIterator>
 361    OutputIterator copy_bernoulli_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol)
 362    {
 363       //
 364       // There are basically 3 thread safety options:
 365       //
 366       // 1) There are no threads (BOOST_HAS_THREADS is not defined).
 367       // 2) There are threads, but we do not have a true atomic integer type,
 368       //    in this case we just use a mutex to guard against race conditions.
 369       // 3) There are threads, and we have an atomic integer: in this case we can
 370       //    use the double-checked locking pattern to avoid thread synchronisation
 371       //    when accessing values already in the cache.
 372       //
 373       // First off handle the common case for overflow and/or asymptotic expansion:
 374       //
 375       if(start + n > bn.capacity())
 376       {
 377          if(start < bn.capacity())
 378          {
 379             out = copy_bernoulli_numbers(out, start, bn.capacity() - start, pol);
 380             n -= bn.capacity() - start;
 381             start = static_cast<std::size_t>(bn.capacity());
 382          }
 383          if(start < b2n_overflow_limit<T, Policy>() + 2u)
 384          {
 385             for(; n; ++start, --n)
 386             {
 387                *out = b2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start * 2U));
 388                ++out;
 389             }
 390          }
 391          for(; n; ++start, --n)
 392          {
 393             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(start), pol);
 394             ++out;
 395          }
 396          return out;
 397       }
 398    #if !defined(BOOST_HAS_THREADS)
 399       //
 400       // Single threaded code, very simple:
 401       //
 402       if(m_current_precision < boost::math::tools::digits<T>())
 403       {
 404          bn.clear();
 405          tn.clear();
 406          m_intermediates.clear();
 407          m_current_precision = boost::math::tools::digits<T>();
 408       }
 409       if(start + n >= bn.size())
 410       {
 411          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()));
 412          tangent_numbers_series(new_size);
 413       }
 414
 415       for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
 416       {
 417          *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
 418          ++out;
 419       }
 420    #elif defined(BOOST_MATH_NO_ATOMIC_INT)
 421       //
 422       // We need to grab a mutex every time we get here, for both readers and writers:
 423       //
 424       boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
 425       if(m_current_precision < boost::math::tools::digits<T>())
 426       {
 427          bn.clear();
 428          tn.clear();
 429          m_intermediates.clear();
 430          m_current_precision = boost::math::tools::digits<T>();
 431       }
 432       if(start + n >= bn.size())
 433       {
 434          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()));
 435          tangent_numbers_series(new_size);
 436       }
 437
 438       for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
 439       {
 440          *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
 441          ++out;
 442       }
 443
 444    #else
 445       //
 446       // Double-checked locking pattern, lets us access cached already cached values
 447       // without locking:
 448       //
 449       // Get the counter and see if we need to calculate more constants:
 450       //
 451       if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
 452          || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()))
 453       {
 454          boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
 455
 456          if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
 457             || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()))
 458          {
 459             if(static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())
 460             {
 461                bn.clear();
 462                tn.clear();
 463                m_intermediates.clear();
 464                m_counter.store(0, BOOST_MATH_ATOMIC_NS::memory_order_release);
 465                m_current_precision = boost::math::tools::digits<T>();
 466             }
 467             if(start + n >= bn.size())
 468             {
 469                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()));
 470                tangent_numbers_series(new_size);
 471             }
 472             m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release);
 473          }
 474       }
 475
 476       for(std::size_t i = (std::max)(static_cast<std::size_t>(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
 477       {
 478          *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)];
 479          ++out;
 480       }
 481
 482    #endif
 483       return out;
 484    }
 485
 486    template <class OutputIterator>
 487    OutputIterator copy_tangent_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol)
 488    {
 489       //
 490       // There are basically 3 thread safety options:
 491       //
 492       // 1) There are no threads (BOOST_HAS_THREADS is not defined).
 493       // 2) There are threads, but we do not have a true atomic integer type,
 494       //    in this case we just use a mutex to guard against race conditions.
 495       // 3) There are threads, and we have an atomic integer: in this case we can
 496       //    use the double-checked locking pattern to avoid thread synchronisation
 497       //    when accessing values already in the cache.
 498       //
 499       //
 500       // First off handle the common case for overflow and/or asymptotic expansion:
 501       //
 502       if(start + n > bn.capacity())
 503       {
 504          if(start < bn.capacity())
 505          {
 506             out = copy_tangent_numbers(out, start, bn.capacity() - start, pol);
 507             n -= bn.capacity() - start;
 508             start = static_cast<std::size_t>(bn.capacity());
 509          }
 510          if(start < b2n_overflow_limit<T, Policy>() + 2u)
 511          {
 512             for(; n; ++start, --n)
 513             {
 514                *out = t2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start));
 515                ++out;
 516             }
 517          }
 518          for(; n; ++start, --n)
 519          {
 520             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(start), pol);
 521             ++out;
 522          }
 523          return out;
 524       }
 525    #if !defined(BOOST_HAS_THREADS)
 526       //
 527       // Single threaded code, very simple:
 528       //
 529       if(m_current_precision < boost::math::tools::digits<T>())
 530       {
 531          bn.clear();
 532          tn.clear();
 533          m_intermediates.clear();
 534          m_current_precision = boost::math::tools::digits<T>();
 535       }
 536       if(start + n >= bn.size())
 537       {
 538          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()));
 539          tangent_numbers_series(new_size);
 540       }
 541
 542       for(std::size_t i = start; i < start + n; ++i)
 543       {
 544          if(i >= m_overflow_limit)
 545             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
 546          else
 547          {
 548             if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
 549                *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
 550             else
 551                *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
 552          }
 553          ++out;
 554       }
 555    #elif defined(BOOST_MATH_NO_ATOMIC_INT)
 556       //
 557       // We need to grab a mutex every time we get here, for both readers and writers:
 558       //
 559       boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
 560       if(m_current_precision < boost::math::tools::digits<T>())
 561       {
 562          bn.clear();
 563          tn.clear();
 564          m_intermediates.clear();
 565          m_current_precision = boost::math::tools::digits<T>();
 566       }
 567       if(start + n >= bn.size())
 568       {
 569          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()));
 570          tangent_numbers_series(new_size);
 571       }
 572
 573       for(std::size_t i = start; i < start + n; ++i)
 574       {
 575          if(i >= m_overflow_limit)
 576             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
 577          else
 578          {
 579             if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
 580                *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
 581             else
 582                *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
 583          }
 584          ++out;
 585       }
 586
 587    #else
 588       //
 589       // Double-checked locking pattern, lets us access cached already cached values
 590       // without locking:
 591       //
 592       // Get the counter and see if we need to calculate more constants:
 593       //
 594       if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
 595          || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()))
 596       {
 597          boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
 598
 599          if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
 600             || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()))
 601          {
 602             if(static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())
 603             {
 604                bn.clear();
 605                tn.clear();
 606                m_intermediates.clear();
 607                m_counter.store(0, BOOST_MATH_ATOMIC_NS::memory_order_release);
 608                m_current_precision = boost::math::tools::digits<T>();
 609             }
 610             if(start + n >= bn.size())
 611             {
 612                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()));
 613                tangent_numbers_series(new_size);
 614             }
 615             m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release);
 616          }
 617       }
 618
 619       for(std::size_t i = start; i < start + n; ++i)
 620       {
 621          if(i >= m_overflow_limit)
 622             *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
 623          else
 624          {
 625             if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
 626                *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
 627             else
 628                *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
 629          }
 630          ++out;
 631       }
 632
 633    #endif
 634       return out;
 635    }
 636
 637 private:
 638    //
 639    // The caches for Bernoulli and tangent numbers, once allocated,
 640    // these must NEVER EVER reallocate as it breaks our thread
 641    // safety guarantees:
 642    //
 643    fixed_vector<T> bn, tn;
 644    std::vector<T> m_intermediates;
 645    // The value at which we know overflow has already occurred for the Bn:
 646    std::size_t m_overflow_limit;
 647 #if !defined(BOOST_HAS_THREADS)
 648    int m_current_precision;
 649 #elif defined(BOOST_MATH_NO_ATOMIC_INT)
 650    boost::detail::lightweight_mutex m_mutex;
 651    int m_current_precision;
 652 #else
 653    boost::detail::lightweight_mutex m_mutex;
 654    atomic_counter_type m_counter, m_current_precision;
 655 #endif
 656 };
 657
 658 template <class T, class Policy>
 659 inline bernoulli_numbers_cache<T, Policy>& get_bernoulli_numbers_cache()
 660 {
 661    //
 662    // Force this function to be called at program startup so all the static variables
 663    // get initailzed then (thread safety).
 664    //
 665    bernoulli_initializer<T, Policy>::force_instantiate();
 666    static bernoulli_numbers_cache<T, Policy> data;
 667    return data;
 668 }
 669
 670 }}}
 671
 672 #endif // BOOST_MATH_BERNOULLI_DETAIL_HPP