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 | |
18 | namespace 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 | // |
23 | template <class T, class Policy> |
24 | T 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 | |
40 | template <class T, class Policy> |
41 | T 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 | // |
71 | struct 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 | } |
90 | private: |
91 | double target; |
92 | }; |
93 | |
94 | template <class T, class Policy> |
95 | inline 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 | |
111 | template <class T, class Policy> |
112 | inline 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 | |
117 | template <class T, class Policy> |
118 | std::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 | // |
132 | template <class T> |
133 | inline 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 | } |
138 | template <class T> |
139 | inline 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 | // |
146 | template <class T, class Policy> |
147 | struct 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 | |
176 | template <class T, class Policy> |
177 | const 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 | // |
189 | template <class T> |
190 | struct 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; } |
237 | private: |
238 | T* m_data; |
239 | unsigned m_used, m_capacity; |
240 | }; |
241 | |
242 | template <class T, class Policy> |
243 | class bernoulli_numbers_cache |
244 | { |
245 | public: |
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 | |
638 | private: |
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 | |
659 | template <class T, class Policy> |
660 | inline 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 | |