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

source code of include/boost/math/special_functions/detail/bernoulli_details.hpp