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