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