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