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