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

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