]>
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 | { | |
195 | for(unsigned i = 0; i < m_used; ++i) | |
196 | this->destroy(&m_data[i]); | |
197 | this->deallocate(m_data, m_capacity); | |
198 | } | |
199 | T& operator[](unsigned n) { BOOST_ASSERT(n < m_used); return m_data[n]; } | |
200 | const T& operator[](unsigned n)const { BOOST_ASSERT(n < m_used); return m_data[n]; } | |
201 | unsigned size()const { return m_used; } | |
202 | unsigned size() { return m_used; } | |
203 | void resize(unsigned n, const T& val) | |
204 | { | |
205 | if(n > m_capacity) | |
206 | { | |
207 | BOOST_THROW_EXCEPTION(std::runtime_error("Exhausted storage for Bernoulli numbers.")); | |
208 | } | |
209 | for(unsigned i = m_used; i < n; ++i) | |
210 | new (m_data + i) T(val); | |
211 | m_used = n; | |
212 | } | |
213 | void resize(unsigned n) { resize(n, T()); } | |
214 | T* begin() { return m_data; } | |
215 | T* end() { return m_data + m_used; } | |
216 | T* begin()const { return m_data; } | |
217 | T* end()const { return m_data + m_used; } | |
218 | unsigned capacity()const { return m_capacity; } | |
219 | void clear() { m_used = 0; } | |
220 | private: | |
221 | T* m_data; | |
222 | unsigned m_used, m_capacity; | |
223 | }; | |
224 | ||
225 | template <class T, class Policy> | |
226 | class bernoulli_numbers_cache | |
227 | { | |
228 | public: | |
229 | bernoulli_numbers_cache() : m_overflow_limit((std::numeric_limits<std::size_t>::max)()) | |
230 | #if defined(BOOST_HAS_THREADS) && !defined(BOOST_MATH_NO_ATOMIC_INT) | |
231 | , m_counter(0) | |
232 | #endif | |
233 | , m_current_precision(boost::math::tools::digits<T>()) | |
234 | {} | |
235 | ||
236 | typedef fixed_vector<T> container_type; | |
237 | ||
238 | void tangent(std::size_t m) | |
239 | { | |
240 | static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1; | |
241 | tn.resize(static_cast<typename container_type::size_type>(m), T(0U)); | |
242 | ||
243 | BOOST_MATH_INSTRUMENT_VARIABLE(min_overflow_index); | |
244 | ||
245 | std::size_t prev_size = m_intermediates.size(); | |
246 | m_intermediates.resize(m, T(0U)); | |
247 | ||
248 | if(prev_size == 0) | |
249 | { | |
250 | m_intermediates[1] = tangent_scale_factor<T>() /*T(1U)*/; | |
251 | tn[0U] = T(0U); | |
252 | tn[1U] = tangent_scale_factor<T>()/* T(1U)*/; | |
253 | BOOST_MATH_INSTRUMENT_VARIABLE(tn[0]); | |
254 | BOOST_MATH_INSTRUMENT_VARIABLE(tn[1]); | |
255 | } | |
256 | ||
257 | for(std::size_t i = std::max<size_t>(2, prev_size); i < m; i++) | |
258 | { | |
259 | bool overflow_check = false; | |
260 | if(i >= min_overflow_index && (boost::math::tools::max_value<T>() / (i-1) < m_intermediates[1]) ) | |
261 | { | |
262 | std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>()); | |
263 | break; | |
264 | } | |
265 | m_intermediates[1] = m_intermediates[1] * (i-1); | |
266 | for(std::size_t j = 2; j <= i; j++) | |
267 | { | |
268 | overflow_check = | |
269 | (i >= min_overflow_index) && ( | |
270 | (boost::math::tools::max_value<T>() / (i - j) < m_intermediates[j]) | |
271 | || (boost::math::tools::max_value<T>() / (i - j + 2) < m_intermediates[j-1]) | |
272 | || (boost::math::tools::max_value<T>() - m_intermediates[j] * (i - j) < m_intermediates[j-1] * (i - j + 2)) | |
273 | || ((boost::math::isinf)(m_intermediates[j])) | |
274 | ); | |
275 | ||
276 | if(overflow_check) | |
277 | { | |
278 | std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>()); | |
279 | break; | |
280 | } | |
281 | m_intermediates[j] = m_intermediates[j] * (i - j) + m_intermediates[j-1] * (i - j + 2); | |
282 | } | |
283 | if(overflow_check) | |
284 | break; // already filled the tn... | |
285 | tn[static_cast<typename container_type::size_type>(i)] = m_intermediates[i]; | |
286 | BOOST_MATH_INSTRUMENT_VARIABLE(i); | |
287 | BOOST_MATH_INSTRUMENT_VARIABLE(tn[static_cast<typename container_type::size_type>(i)]); | |
288 | } | |
289 | } | |
290 | ||
291 | void tangent_numbers_series(const std::size_t m) | |
292 | { | |
293 | BOOST_MATH_STD_USING | |
294 | static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1; | |
295 | ||
296 | typename container_type::size_type old_size = bn.size(); | |
297 | ||
298 | tangent(m); | |
299 | bn.resize(static_cast<typename container_type::size_type>(m)); | |
300 | ||
301 | if(!old_size) | |
302 | { | |
303 | bn[0] = 1; | |
304 | old_size = 1; | |
305 | } | |
306 | ||
307 | T power_two(ldexp(T(1), static_cast<int>(2 * old_size))); | |
308 | ||
309 | for(std::size_t i = old_size; i < m; i++) | |
310 | { | |
311 | T b(static_cast<T>(i * 2)); | |
312 | // | |
313 | // Not only do we need to take care to avoid spurious over/under flow in | |
314 | // the calculation, but we also need to avoid overflow altogether in case | |
315 | // we're calculating with a type where "bad things" happen in that case: | |
316 | // | |
317 | b = b / (power_two * tangent_scale_factor<T>()); | |
318 | b /= (power_two - 1); | |
319 | bool overflow_check = (i >= min_overflow_index) && (tools::max_value<T>() / tn[static_cast<typename container_type::size_type>(i)] < b); | |
320 | if(overflow_check) | |
321 | { | |
322 | m_overflow_limit = i; | |
323 | while(i < m) | |
324 | { | |
325 | b = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : tools::max_value<T>(); | |
326 | bn[static_cast<typename container_type::size_type>(i)] = ((i % 2U) ? b : T(-b)); | |
327 | ++i; | |
328 | } | |
329 | break; | |
330 | } | |
331 | else | |
332 | { | |
333 | b *= tn[static_cast<typename container_type::size_type>(i)]; | |
334 | } | |
335 | ||
336 | power_two = ldexp(power_two, 2); | |
337 | ||
338 | const bool b_neg = i % 2 == 0; | |
339 | ||
340 | bn[static_cast<typename container_type::size_type>(i)] = ((!b_neg) ? b : T(-b)); | |
341 | } | |
342 | } | |
343 | ||
344 | template <class OutputIterator> | |
345 | OutputIterator copy_bernoulli_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol) | |
346 | { | |
347 | // | |
348 | // There are basically 3 thread safety options: | |
349 | // | |
350 | // 1) There are no threads (BOOST_HAS_THREADS is not defined). | |
351 | // 2) There are threads, but we do not have a true atomic integer type, | |
352 | // in this case we just use a mutex to guard against race conditions. | |
353 | // 3) There are threads, and we have an atomic integer: in this case we can | |
354 | // use the double-checked locking pattern to avoid thread synchronisation | |
355 | // when accessing values already in the cache. | |
356 | // | |
357 | // First off handle the common case for overflow and/or asymptotic expansion: | |
358 | // | |
359 | if(start + n > bn.capacity()) | |
360 | { | |
361 | if(start < bn.capacity()) | |
362 | { | |
363 | out = copy_bernoulli_numbers(out, start, bn.capacity() - start, pol); | |
364 | n -= bn.capacity() - start; | |
365 | start = static_cast<std::size_t>(bn.capacity()); | |
366 | } | |
367 | if(start < b2n_overflow_limit<T, Policy>() + 2u) | |
368 | { | |
369 | for(; n; ++start, --n) | |
370 | { | |
371 | *out = b2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start * 2U)); | |
372 | ++out; | |
373 | } | |
374 | } | |
375 | for(; n; ++start, --n) | |
376 | { | |
377 | *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(start), pol); | |
378 | ++out; | |
379 | } | |
380 | return out; | |
381 | } | |
382 | #if !defined(BOOST_HAS_THREADS) | |
383 | // | |
384 | // Single threaded code, very simple: | |
385 | // | |
386 | if(m_current_precision < boost::math::tools::digits<T>()) | |
387 | { | |
388 | bn.clear(); | |
389 | tn.clear(); | |
390 | m_intermediates.clear(); | |
391 | m_current_precision = boost::math::tools::digits<T>(); | |
392 | } | |
393 | if(start + n >= bn.size()) | |
394 | { | |
395 | 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())); | |
396 | tangent_numbers_series(new_size); | |
397 | } | |
398 | ||
399 | for(std::size_t i = (std::max)(std::size_t(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i) | |
400 | { | |
401 | *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i]; | |
402 | ++out; | |
403 | } | |
404 | #elif defined(BOOST_MATH_NO_ATOMIC_INT) | |
405 | // | |
406 | // We need to grab a mutex every time we get here, for both readers and writers: | |
407 | // | |
408 | boost::detail::lightweight_mutex::scoped_lock l(m_mutex); | |
409 | if(m_current_precision < boost::math::tools::digits<T>()) | |
410 | { | |
411 | bn.clear(); | |
412 | tn.clear(); | |
413 | m_intermediates.clear(); | |
414 | m_current_precision = boost::math::tools::digits<T>(); | |
415 | } | |
416 | if(start + n >= bn.size()) | |
417 | { | |
418 | 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())); | |
419 | tangent_numbers_series(new_size); | |
420 | } | |
421 | ||
422 | for(std::size_t i = (std::max)(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[i]; | |
425 | ++out; | |
426 | } | |
427 | ||
428 | #else | |
429 | // | |
430 | // Double-checked locking pattern, lets us access cached already cached values | |
431 | // without locking: | |
432 | // | |
433 | // Get the counter and see if we need to calculate more constants: | |
434 | // | |
435 | if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n) | |
436 | || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())) | |
437 | { | |
438 | boost::detail::lightweight_mutex::scoped_lock l(m_mutex); | |
439 | ||
440 | if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n) | |
441 | || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())) | |
442 | { | |
443 | if(static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()) | |
444 | { | |
445 | bn.clear(); | |
446 | tn.clear(); | |
447 | m_intermediates.clear(); | |
448 | m_counter.store(0, BOOST_MATH_ATOMIC_NS::memory_order_release); | |
449 | m_current_precision = boost::math::tools::digits<T>(); | |
450 | } | |
451 | if(start + n >= bn.size()) | |
452 | { | |
453 | 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())); | |
454 | tangent_numbers_series(new_size); | |
455 | } | |
456 | m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release); | |
457 | } | |
458 | } | |
459 | ||
460 | for(std::size_t i = (std::max)(static_cast<std::size_t>(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i) | |
461 | { | |
462 | *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)]; | |
463 | ++out; | |
464 | } | |
465 | ||
466 | #endif | |
467 | return out; | |
468 | } | |
469 | ||
470 | template <class OutputIterator> | |
471 | OutputIterator copy_tangent_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol) | |
472 | { | |
473 | // | |
474 | // There are basically 3 thread safety options: | |
475 | // | |
476 | // 1) There are no threads (BOOST_HAS_THREADS is not defined). | |
477 | // 2) There are threads, but we do not have a true atomic integer type, | |
478 | // in this case we just use a mutex to guard against race conditions. | |
479 | // 3) There are threads, and we have an atomic integer: in this case we can | |
480 | // use the double-checked locking pattern to avoid thread synchronisation | |
481 | // when accessing values already in the cache. | |
482 | // | |
483 | // | |
484 | // First off handle the common case for overflow and/or asymptotic expansion: | |
485 | // | |
486 | if(start + n > bn.capacity()) | |
487 | { | |
488 | if(start < bn.capacity()) | |
489 | { | |
490 | out = copy_tangent_numbers(out, start, bn.capacity() - start, pol); | |
491 | n -= bn.capacity() - start; | |
492 | start = static_cast<std::size_t>(bn.capacity()); | |
493 | } | |
494 | if(start < b2n_overflow_limit<T, Policy>() + 2u) | |
495 | { | |
496 | for(; n; ++start, --n) | |
497 | { | |
498 | *out = t2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start)); | |
499 | ++out; | |
500 | } | |
501 | } | |
502 | for(; n; ++start, --n) | |
503 | { | |
504 | *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(start), pol); | |
505 | ++out; | |
506 | } | |
507 | return out; | |
508 | } | |
509 | #if !defined(BOOST_HAS_THREADS) | |
510 | // | |
511 | // Single threaded code, very simple: | |
512 | // | |
513 | if(m_current_precision < boost::math::tools::digits<T>()) | |
514 | { | |
515 | bn.clear(); | |
516 | tn.clear(); | |
517 | m_intermediates.clear(); | |
518 | m_current_precision = boost::math::tools::digits<T>(); | |
519 | } | |
520 | if(start + n >= bn.size()) | |
521 | { | |
522 | 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())); | |
523 | tangent_numbers_series(new_size); | |
524 | } | |
525 | ||
526 | for(std::size_t i = start; i < start + n; ++i) | |
527 | { | |
528 | if(i >= m_overflow_limit) | |
529 | *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol); | |
530 | else | |
531 | { | |
532 | if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)]) | |
533 | *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol); | |
534 | else | |
535 | *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>(); | |
536 | } | |
537 | ++out; | |
538 | } | |
539 | #elif defined(BOOST_MATH_NO_ATOMIC_INT) | |
540 | // | |
541 | // We need to grab a mutex every time we get here, for both readers and writers: | |
542 | // | |
543 | boost::detail::lightweight_mutex::scoped_lock l(m_mutex); | |
544 | if(m_current_precision < boost::math::tools::digits<T>()) | |
545 | { | |
546 | bn.clear(); | |
547 | tn.clear(); | |
548 | m_intermediates.clear(); | |
549 | m_current_precision = boost::math::tools::digits<T>(); | |
550 | } | |
551 | if(start + n >= bn.size()) | |
552 | { | |
553 | 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())); | |
554 | tangent_numbers_series(new_size); | |
555 | } | |
556 | ||
557 | for(std::size_t i = start; i < start + n; ++i) | |
558 | { | |
559 | if(i >= m_overflow_limit) | |
560 | *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol); | |
561 | else | |
562 | { | |
563 | if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)]) | |
564 | *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol); | |
565 | else | |
566 | *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>(); | |
567 | } | |
568 | ++out; | |
569 | } | |
570 | ||
571 | #else | |
572 | // | |
573 | // Double-checked locking pattern, lets us access cached already cached values | |
574 | // without locking: | |
575 | // | |
576 | // Get the counter and see if we need to calculate more constants: | |
577 | // | |
578 | if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n) | |
579 | || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())) | |
580 | { | |
581 | boost::detail::lightweight_mutex::scoped_lock l(m_mutex); | |
582 | ||
583 | if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n) | |
584 | || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())) | |
585 | { | |
586 | if(static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()) | |
587 | { | |
588 | bn.clear(); | |
589 | tn.clear(); | |
590 | m_intermediates.clear(); | |
591 | m_counter.store(0, BOOST_MATH_ATOMIC_NS::memory_order_release); | |
592 | m_current_precision = boost::math::tools::digits<T>(); | |
593 | } | |
594 | if(start + n >= bn.size()) | |
595 | { | |
596 | 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())); | |
597 | tangent_numbers_series(new_size); | |
598 | } | |
599 | m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release); | |
600 | } | |
601 | } | |
602 | ||
603 | for(std::size_t i = start; i < start + n; ++i) | |
604 | { | |
605 | if(i >= m_overflow_limit) | |
606 | *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol); | |
607 | else | |
608 | { | |
609 | if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)]) | |
610 | *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol); | |
611 | else | |
612 | *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>(); | |
613 | } | |
614 | ++out; | |
615 | } | |
616 | ||
617 | #endif | |
618 | return out; | |
619 | } | |
620 | ||
621 | private: | |
622 | // | |
623 | // The caches for Bernoulli and tangent numbers, once allocated, | |
624 | // these must NEVER EVER reallocate as it breaks our thread | |
625 | // safety guarantees: | |
626 | // | |
627 | fixed_vector<T> bn, tn; | |
628 | std::vector<T> m_intermediates; | |
629 | // The value at which we know overflow has already occurred for the Bn: | |
630 | std::size_t m_overflow_limit; | |
631 | #if !defined(BOOST_HAS_THREADS) | |
632 | int m_current_precision; | |
633 | #elif defined(BOOST_MATH_NO_ATOMIC_INT) | |
634 | boost::detail::lightweight_mutex m_mutex; | |
635 | int m_current_precision; | |
636 | #else | |
637 | boost::detail::lightweight_mutex m_mutex; | |
638 | atomic_counter_type m_counter, m_current_precision; | |
639 | #endif | |
640 | }; | |
641 | ||
642 | template <class T, class Policy> | |
643 | inline bernoulli_numbers_cache<T, Policy>& get_bernoulli_numbers_cache() | |
644 | { | |
645 | // | |
646 | // Force this function to be called at program startup so all the static variables | |
647 | // get initailzed then (thread safety). | |
648 | // | |
649 | bernoulli_initializer<T, Policy>::force_instantiate(); | |
650 | static bernoulli_numbers_cache<T, Policy> data; | |
651 | return data; | |
652 | } | |
653 | ||
654 | }}} | |
655 | ||
656 | #endif // BOOST_MATH_BERNOULLI_DETAIL_HPP |