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