]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright John Maddock 2010, 2012. |
2 | // Copyright Paul A. Bristow 2011, 2012. | |
3 | ||
4 | // Use, modification and distribution are subject to the | |
5 | // Boost Software License, Version 1.0. (See accompanying file | |
6 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
7 | ||
8 | #ifndef BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED | |
9 | #define BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED | |
1e59de90 | 10 | #include <type_traits> |
7c673cae FG |
11 | |
12 | namespace boost{ namespace math{ namespace constants{ namespace detail{ | |
13 | ||
14 | template <class T> | |
15 | template<int N> | |
1e59de90 | 16 | inline T constant_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
17 | { |
18 | BOOST_MATH_STD_USING | |
19 | ||
20 | return ldexp(acos(T(0)), 1); | |
21 | ||
22 | /* | |
23 | // Although this code works well, it's usually more accurate to just call acos | |
24 | // and access the number types own representation of PI which is usually calculated | |
25 | // at slightly higher precision... | |
26 | ||
27 | T result; | |
28 | T a = 1; | |
29 | T b; | |
30 | T A(a); | |
31 | T B = 0.5f; | |
32 | T D = 0.25f; | |
33 | ||
34 | T lim; | |
35 | lim = boost::math::tools::epsilon<T>(); | |
36 | ||
37 | unsigned k = 1; | |
38 | ||
39 | do | |
40 | { | |
41 | result = A + B; | |
42 | result = ldexp(result, -2); | |
43 | b = sqrt(B); | |
44 | a += b; | |
45 | a = ldexp(a, -1); | |
46 | A = a * a; | |
47 | B = A - result; | |
48 | B = ldexp(B, 1); | |
49 | result = A - B; | |
50 | bool neg = boost::math::sign(result) < 0; | |
51 | if(neg) | |
52 | result = -result; | |
53 | if(result <= lim) | |
54 | break; | |
55 | if(neg) | |
56 | result = -result; | |
57 | result = ldexp(result, k - 1); | |
58 | D -= result; | |
59 | ++k; | |
60 | lim = ldexp(lim, 1); | |
61 | } | |
62 | while(true); | |
63 | ||
64 | result = B / D; | |
65 | return result; | |
66 | */ | |
67 | } | |
68 | ||
69 | template <class T> | |
70 | template<int N> | |
1e59de90 | 71 | inline T constant_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
72 | { |
73 | return 2 * pi<T, policies::policy<policies::digits2<N> > >(); | |
74 | } | |
75 | ||
76 | template <class T> // 2 / pi | |
77 | template<int N> | |
1e59de90 | 78 | inline T constant_two_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
79 | { |
80 | return 2 / pi<T, policies::policy<policies::digits2<N> > >(); | |
81 | } | |
82 | ||
83 | template <class T> // sqrt(2/pi) | |
84 | template <int N> | |
1e59de90 | 85 | inline T constant_root_two_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
86 | { |
87 | BOOST_MATH_STD_USING | |
88 | return sqrt((2 / pi<T, policies::policy<policies::digits2<N> > >())); | |
89 | } | |
90 | ||
91 | template <class T> | |
92 | template<int N> | |
1e59de90 | 93 | inline T constant_one_div_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
94 | { |
95 | return 1 / two_pi<T, policies::policy<policies::digits2<N> > >(); | |
96 | } | |
97 | ||
98 | template <class T> | |
99 | template<int N> | |
1e59de90 | 100 | inline T constant_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
101 | { |
102 | BOOST_MATH_STD_USING | |
103 | return sqrt(pi<T, policies::policy<policies::digits2<N> > >()); | |
104 | } | |
105 | ||
106 | template <class T> | |
107 | template<int N> | |
1e59de90 | 108 | inline T constant_root_half_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
109 | { |
110 | BOOST_MATH_STD_USING | |
111 | return sqrt(pi<T, policies::policy<policies::digits2<N> > >() / 2); | |
112 | } | |
113 | ||
114 | template <class T> | |
115 | template<int N> | |
1e59de90 | 116 | inline T constant_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
117 | { |
118 | BOOST_MATH_STD_USING | |
119 | return sqrt(two_pi<T, policies::policy<policies::digits2<N> > >()); | |
120 | } | |
121 | ||
122 | template <class T> | |
123 | template<int N> | |
1e59de90 | 124 | inline T constant_log_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
125 | { |
126 | BOOST_MATH_STD_USING | |
127 | return log(root_two_pi<T, policies::policy<policies::digits2<N> > >()); | |
128 | } | |
129 | ||
130 | template <class T> | |
131 | template<int N> | |
1e59de90 | 132 | inline T constant_root_ln_four<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
133 | { |
134 | BOOST_MATH_STD_USING | |
135 | return sqrt(log(static_cast<T>(4))); | |
136 | } | |
137 | ||
138 | template <class T> | |
139 | template<int N> | |
1e59de90 | 140 | inline T constant_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
141 | { |
142 | // | |
143 | // Although we can clearly calculate this from first principles, this hooks into | |
144 | // T's own notion of e, which hopefully will more accurate than one calculated to | |
145 | // a few epsilon: | |
146 | // | |
147 | BOOST_MATH_STD_USING | |
148 | return exp(static_cast<T>(1)); | |
149 | } | |
150 | ||
151 | template <class T> | |
152 | template<int N> | |
1e59de90 | 153 | inline T constant_half<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
154 | { |
155 | return static_cast<T>(1) / static_cast<T>(2); | |
156 | } | |
157 | ||
158 | template <class T> | |
159 | template<int M> | |
1e59de90 | 160 | inline T constant_euler<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, M>))) |
7c673cae FG |
161 | { |
162 | BOOST_MATH_STD_USING | |
163 | // | |
164 | // This is the method described in: | |
165 | // "Some New Algorithms for High-Precision Computation of Euler's Constant" | |
166 | // Richard P Brent and Edwin M McMillan. | |
167 | // Mathematics of Computation, Volume 34, Number 149, Jan 1980, pages 305-312. | |
168 | // See equation 17 with p = 2. | |
169 | // | |
170 | T n = 3 + (M ? (std::min)(M, tools::digits<T>()) : tools::digits<T>()) / 4; | |
171 | T lim = M ? ldexp(T(1), 1 - (std::min)(M, tools::digits<T>())) : tools::epsilon<T>(); | |
172 | T lnn = log(n); | |
173 | T term = 1; | |
174 | T N = -lnn; | |
175 | T D = 1; | |
176 | T Hk = 0; | |
177 | T one = 1; | |
178 | ||
179 | for(unsigned k = 1;; ++k) | |
180 | { | |
181 | term *= n * n; | |
182 | term /= k * k; | |
183 | Hk += one / k; | |
184 | N += term * (Hk - lnn); | |
185 | D += term; | |
186 | ||
187 | if(term < D * lim) | |
188 | break; | |
189 | } | |
190 | return N / D; | |
191 | } | |
192 | ||
193 | template <class T> | |
194 | template<int N> | |
1e59de90 | 195 | inline T constant_euler_sqr<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
196 | { |
197 | BOOST_MATH_STD_USING | |
198 | return euler<T, policies::policy<policies::digits2<N> > >() | |
199 | * euler<T, policies::policy<policies::digits2<N> > >(); | |
200 | } | |
201 | ||
202 | template <class T> | |
203 | template<int N> | |
1e59de90 | 204 | inline T constant_one_div_euler<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
205 | { |
206 | BOOST_MATH_STD_USING | |
207 | return static_cast<T>(1) | |
208 | / euler<T, policies::policy<policies::digits2<N> > >(); | |
209 | } | |
210 | ||
211 | ||
212 | template <class T> | |
213 | template<int N> | |
1e59de90 | 214 | inline T constant_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
215 | { |
216 | BOOST_MATH_STD_USING | |
217 | return sqrt(static_cast<T>(2)); | |
218 | } | |
219 | ||
220 | ||
221 | template <class T> | |
222 | template<int N> | |
1e59de90 | 223 | inline T constant_root_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
224 | { |
225 | BOOST_MATH_STD_USING | |
226 | return sqrt(static_cast<T>(3)); | |
227 | } | |
228 | ||
229 | template <class T> | |
230 | template<int N> | |
1e59de90 | 231 | inline T constant_half_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
232 | { |
233 | BOOST_MATH_STD_USING | |
234 | return sqrt(static_cast<T>(2)) / 2; | |
235 | } | |
236 | ||
237 | template <class T> | |
238 | template<int N> | |
1e59de90 | 239 | inline T constant_ln_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
240 | { |
241 | // | |
242 | // Although there are good ways to calculate this from scratch, this hooks into | |
243 | // T's own notion of log(2) which will hopefully be accurate to the full precision | |
244 | // of T: | |
245 | // | |
246 | BOOST_MATH_STD_USING | |
247 | return log(static_cast<T>(2)); | |
248 | } | |
249 | ||
250 | template <class T> | |
251 | template<int N> | |
1e59de90 | 252 | inline T constant_ln_ten<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
253 | { |
254 | BOOST_MATH_STD_USING | |
255 | return log(static_cast<T>(10)); | |
256 | } | |
257 | ||
258 | template <class T> | |
259 | template<int N> | |
1e59de90 | 260 | inline T constant_ln_ln_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
261 | { |
262 | BOOST_MATH_STD_USING | |
263 | return log(log(static_cast<T>(2))); | |
264 | } | |
265 | ||
266 | template <class T> | |
267 | template<int N> | |
1e59de90 | 268 | inline T constant_third<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
269 | { |
270 | BOOST_MATH_STD_USING | |
271 | return static_cast<T>(1) / static_cast<T>(3); | |
272 | } | |
273 | ||
274 | template <class T> | |
275 | template<int N> | |
1e59de90 | 276 | inline T constant_twothirds<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
277 | { |
278 | BOOST_MATH_STD_USING | |
279 | return static_cast<T>(2) / static_cast<T>(3); | |
280 | } | |
281 | ||
282 | template <class T> | |
283 | template<int N> | |
1e59de90 | 284 | inline T constant_two_thirds<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
285 | { |
286 | BOOST_MATH_STD_USING | |
287 | return static_cast<T>(2) / static_cast<T>(3); | |
288 | } | |
289 | ||
290 | template <class T> | |
291 | template<int N> | |
1e59de90 | 292 | inline T constant_three_quarters<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
293 | { |
294 | BOOST_MATH_STD_USING | |
295 | return static_cast<T>(3) / static_cast<T>(4); | |
296 | } | |
297 | ||
92f5a8d4 TL |
298 | template <class T> |
299 | template<int N> | |
1e59de90 | 300 | inline T constant_sixth<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
92f5a8d4 TL |
301 | { |
302 | BOOST_MATH_STD_USING | |
303 | return static_cast<T>(1) / static_cast<T>(6); | |
304 | } | |
305 | ||
306 | // Pi and related constants. | |
7c673cae FG |
307 | template <class T> |
308 | template<int N> | |
1e59de90 | 309 | inline T constant_pi_minus_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
310 | { |
311 | return pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(3); | |
312 | } | |
313 | ||
314 | template <class T> | |
315 | template<int N> | |
1e59de90 | 316 | inline T constant_four_minus_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
317 | { |
318 | return static_cast<T>(4) - pi<T, policies::policy<policies::digits2<N> > >(); | |
319 | } | |
320 | ||
7c673cae FG |
321 | template <class T> |
322 | template<int N> | |
1e59de90 | 323 | inline T constant_exp_minus_half<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
324 | { |
325 | BOOST_MATH_STD_USING | |
326 | return exp(static_cast<T>(-0.5)); | |
327 | } | |
328 | ||
92f5a8d4 TL |
329 | template <class T> |
330 | template<int N> | |
1e59de90 | 331 | inline T constant_exp_minus_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
92f5a8d4 TL |
332 | { |
333 | BOOST_MATH_STD_USING | |
334 | return exp(static_cast<T>(-1.)); | |
335 | } | |
336 | ||
7c673cae FG |
337 | template <class T> |
338 | template<int N> | |
1e59de90 | 339 | inline T constant_one_div_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
340 | { |
341 | return static_cast<T>(1) / root_two<T, policies::policy<policies::digits2<N> > >(); | |
342 | } | |
343 | ||
344 | template <class T> | |
345 | template<int N> | |
1e59de90 | 346 | inline T constant_one_div_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
347 | { |
348 | return static_cast<T>(1) / root_pi<T, policies::policy<policies::digits2<N> > >(); | |
349 | } | |
350 | ||
351 | template <class T> | |
352 | template<int N> | |
1e59de90 | 353 | inline T constant_one_div_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
354 | { |
355 | return static_cast<T>(1) / root_two_pi<T, policies::policy<policies::digits2<N> > >(); | |
356 | } | |
357 | ||
358 | template <class T> | |
359 | template<int N> | |
1e59de90 | 360 | inline T constant_root_one_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
361 | { |
362 | BOOST_MATH_STD_USING | |
363 | return sqrt(static_cast<T>(1) / pi<T, policies::policy<policies::digits2<N> > >()); | |
364 | } | |
365 | ||
7c673cae FG |
366 | template <class T> |
367 | template<int N> | |
1e59de90 | 368 | inline T constant_four_thirds_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
369 | { |
370 | BOOST_MATH_STD_USING | |
371 | return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(4) / static_cast<T>(3); | |
372 | } | |
373 | ||
374 | template <class T> | |
375 | template<int N> | |
1e59de90 | 376 | inline T constant_half_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
377 | { |
378 | BOOST_MATH_STD_USING | |
379 | return pi<T, policies::policy<policies::digits2<N> > >() / static_cast<T>(2); | |
380 | } | |
381 | ||
7c673cae FG |
382 | template <class T> |
383 | template<int N> | |
1e59de90 | 384 | inline T constant_third_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
385 | { |
386 | BOOST_MATH_STD_USING | |
387 | return pi<T, policies::policy<policies::digits2<N> > >() / static_cast<T>(3); | |
388 | } | |
389 | ||
390 | template <class T> | |
391 | template<int N> | |
1e59de90 | 392 | inline T constant_sixth_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
393 | { |
394 | BOOST_MATH_STD_USING | |
395 | return pi<T, policies::policy<policies::digits2<N> > >() / static_cast<T>(6); | |
396 | } | |
397 | ||
398 | template <class T> | |
399 | template<int N> | |
1e59de90 | 400 | inline T constant_two_thirds_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
401 | { |
402 | BOOST_MATH_STD_USING | |
403 | return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(2) / static_cast<T>(3); | |
404 | } | |
405 | ||
406 | template <class T> | |
407 | template<int N> | |
1e59de90 | 408 | inline T constant_three_quarters_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
409 | { |
410 | BOOST_MATH_STD_USING | |
411 | return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(3) / static_cast<T>(4); | |
412 | } | |
413 | ||
414 | template <class T> | |
415 | template<int N> | |
1e59de90 | 416 | inline T constant_pi_pow_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
417 | { |
418 | BOOST_MATH_STD_USING | |
419 | return pow(pi<T, policies::policy<policies::digits2<N> > >(), e<T, policies::policy<policies::digits2<N> > >()); // | |
420 | } | |
421 | ||
422 | template <class T> | |
423 | template<int N> | |
1e59de90 | 424 | inline T constant_pi_sqr<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
425 | { |
426 | BOOST_MATH_STD_USING | |
427 | return pi<T, policies::policy<policies::digits2<N> > >() | |
428 | * pi<T, policies::policy<policies::digits2<N> > >() ; // | |
429 | } | |
430 | ||
431 | template <class T> | |
432 | template<int N> | |
1e59de90 | 433 | inline T constant_pi_sqr_div_six<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
434 | { |
435 | BOOST_MATH_STD_USING | |
436 | return pi<T, policies::policy<policies::digits2<N> > >() | |
437 | * pi<T, policies::policy<policies::digits2<N> > >() | |
438 | / static_cast<T>(6); // | |
439 | } | |
440 | ||
7c673cae FG |
441 | template <class T> |
442 | template<int N> | |
1e59de90 | 443 | inline T constant_pi_cubed<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
444 | { |
445 | BOOST_MATH_STD_USING | |
446 | return pi<T, policies::policy<policies::digits2<N> > >() | |
447 | * pi<T, policies::policy<policies::digits2<N> > >() | |
448 | * pi<T, policies::policy<policies::digits2<N> > >() | |
449 | ; // | |
450 | } | |
451 | ||
452 | template <class T> | |
453 | template<int N> | |
1e59de90 | 454 | inline T constant_cbrt_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
455 | { |
456 | BOOST_MATH_STD_USING | |
457 | return pow(pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(1)/ static_cast<T>(3)); | |
458 | } | |
459 | ||
460 | template <class T> | |
461 | template<int N> | |
1e59de90 | 462 | inline T constant_one_div_cbrt_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
463 | { |
464 | BOOST_MATH_STD_USING | |
465 | return static_cast<T>(1) | |
466 | / pow(pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(1)/ static_cast<T>(3)); | |
467 | } | |
468 | ||
469 | // Euler's e | |
470 | ||
471 | template <class T> | |
472 | template<int N> | |
1e59de90 | 473 | inline T constant_e_pow_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
474 | { |
475 | BOOST_MATH_STD_USING | |
476 | return pow(e<T, policies::policy<policies::digits2<N> > >(), pi<T, policies::policy<policies::digits2<N> > >()); // | |
477 | } | |
478 | ||
479 | template <class T> | |
480 | template<int N> | |
1e59de90 | 481 | inline T constant_root_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
482 | { |
483 | BOOST_MATH_STD_USING | |
484 | return sqrt(e<T, policies::policy<policies::digits2<N> > >()); | |
485 | } | |
486 | ||
487 | template <class T> | |
488 | template<int N> | |
1e59de90 | 489 | inline T constant_log10_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
490 | { |
491 | BOOST_MATH_STD_USING | |
492 | return log10(e<T, policies::policy<policies::digits2<N> > >()); | |
493 | } | |
494 | ||
495 | template <class T> | |
496 | template<int N> | |
1e59de90 | 497 | inline T constant_one_div_log10_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
498 | { |
499 | BOOST_MATH_STD_USING | |
500 | return static_cast<T>(1) / | |
501 | log10(e<T, policies::policy<policies::digits2<N> > >()); | |
502 | } | |
503 | ||
504 | // Trigonometric | |
505 | ||
506 | template <class T> | |
507 | template<int N> | |
1e59de90 | 508 | inline T constant_degree<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
509 | { |
510 | BOOST_MATH_STD_USING | |
511 | return pi<T, policies::policy<policies::digits2<N> > >() | |
512 | / static_cast<T>(180) | |
513 | ; // | |
514 | } | |
515 | ||
516 | template <class T> | |
517 | template<int N> | |
1e59de90 | 518 | inline T constant_radian<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
519 | { |
520 | BOOST_MATH_STD_USING | |
521 | return static_cast<T>(180) | |
522 | / pi<T, policies::policy<policies::digits2<N> > >() | |
523 | ; // | |
524 | } | |
525 | ||
526 | template <class T> | |
527 | template<int N> | |
1e59de90 | 528 | inline T constant_sin_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
529 | { |
530 | BOOST_MATH_STD_USING | |
531 | return sin(static_cast<T>(1)) ; // | |
532 | } | |
533 | ||
534 | template <class T> | |
535 | template<int N> | |
1e59de90 | 536 | inline T constant_cos_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
537 | { |
538 | BOOST_MATH_STD_USING | |
539 | return cos(static_cast<T>(1)) ; // | |
540 | } | |
541 | ||
542 | template <class T> | |
543 | template<int N> | |
1e59de90 | 544 | inline T constant_sinh_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
545 | { |
546 | BOOST_MATH_STD_USING | |
547 | return sinh(static_cast<T>(1)) ; // | |
548 | } | |
549 | ||
550 | template <class T> | |
551 | template<int N> | |
1e59de90 | 552 | inline T constant_cosh_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
553 | { |
554 | BOOST_MATH_STD_USING | |
555 | return cosh(static_cast<T>(1)) ; // | |
556 | } | |
557 | ||
558 | template <class T> | |
559 | template<int N> | |
1e59de90 | 560 | inline T constant_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
561 | { |
562 | BOOST_MATH_STD_USING | |
563 | return (static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) ; // | |
564 | } | |
565 | ||
566 | template <class T> | |
567 | template<int N> | |
1e59de90 | 568 | inline T constant_ln_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
569 | { |
570 | BOOST_MATH_STD_USING | |
7c673cae FG |
571 | return log((static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) ); |
572 | } | |
20effc67 | 573 | |
7c673cae FG |
574 | template <class T> |
575 | template<int N> | |
1e59de90 | 576 | inline T constant_one_div_ln_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
577 | { |
578 | BOOST_MATH_STD_USING | |
579 | return static_cast<T>(1) / | |
580 | log((static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) ); | |
581 | } | |
582 | ||
583 | // Zeta | |
584 | ||
585 | template <class T> | |
586 | template<int N> | |
1e59de90 | 587 | inline T constant_zeta_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
588 | { |
589 | BOOST_MATH_STD_USING | |
590 | ||
591 | return pi<T, policies::policy<policies::digits2<N> > >() | |
592 | * pi<T, policies::policy<policies::digits2<N> > >() | |
593 | /static_cast<T>(6); | |
594 | } | |
595 | ||
596 | template <class T> | |
597 | template<int N> | |
1e59de90 | 598 | inline T constant_zeta_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
599 | { |
600 | // http://mathworld.wolfram.com/AperysConstant.html | |
601 | // http://en.wikipedia.org/wiki/Mathematical_constant | |
602 | ||
603 | // http://oeis.org/A002117/constant | |
604 | //T zeta3("1.20205690315959428539973816151144999076" | |
605 | // "4986292340498881792271555341838205786313" | |
606 | // "09018645587360933525814619915"); | |
607 | ||
608 | //"1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915" A002117 | |
609 | // 1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915780, +00); | |
610 | //"1.2020569031595942 double | |
611 | // http://www.spaennare.se/SSPROG/ssnum.pdf // section 11, Algorithm for Apery's constant zeta(3). | |
612 | // Programs to Calculate some Mathematical Constants to Large Precision, Document Version 1.50 | |
613 | ||
614 | // by Stefan Spannare September 19, 2007 | |
615 | // zeta(3) = 1/64 * sum | |
616 | BOOST_MATH_STD_USING | |
617 | T n_fact=static_cast<T>(1); // build n! for n = 0. | |
618 | T sum = static_cast<double>(77); // Start with n = 0 case. | |
619 | // for n = 0, (77/1) /64 = 1.203125 | |
620 | //double lim = std::numeric_limits<double>::epsilon(); | |
621 | T lim = N ? ldexp(T(1), 1 - (std::min)(N, tools::digits<T>())) : tools::epsilon<T>(); | |
622 | for(unsigned int n = 1; n < 40; ++n) | |
623 | { // three to five decimal digits per term, so 40 should be plenty for 100 decimal digits. | |
624 | //cout << "n = " << n << endl; | |
625 | n_fact *= n; // n! | |
626 | T n_fact_p10 = n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact; // (n!)^10 | |
627 | T num = ((205 * n * n) + (250 * n) + 77) * n_fact_p10; // 205n^2 + 250n + 77 | |
628 | // int nn = (2 * n + 1); | |
629 | // T d = factorial(nn); // inline factorial. | |
630 | T d = 1; | |
631 | for(unsigned int i = 1; i <= (n+n + 1); ++i) // (2n + 1) | |
632 | { | |
633 | d *= i; | |
634 | } | |
635 | T den = d * d * d * d * d; // [(2n+1)!]^5 | |
636 | //cout << "den = " << den << endl; | |
637 | T term = num/den; | |
638 | if (n % 2 != 0) | |
639 | { //term *= -1; | |
640 | sum -= term; | |
641 | } | |
642 | else | |
643 | { | |
644 | sum += term; | |
645 | } | |
646 | //cout << "term = " << term << endl; | |
647 | //cout << "sum/64 = " << sum/64 << endl; | |
648 | if(abs(term) < lim) | |
649 | { | |
650 | break; | |
651 | } | |
652 | } | |
653 | return sum / 64; | |
654 | } | |
655 | ||
656 | template <class T> | |
657 | template<int N> | |
1e59de90 | 658 | inline T constant_catalan<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
659 | { // http://oeis.org/A006752/constant |
660 | //T c("0.915965594177219015054603514932384110774" | |
661 | //"149374281672134266498119621763019776254769479356512926115106248574"); | |
662 | ||
663 | // 9.159655941772190150546035149323841107, 74149374281672134266498119621763019776254769479356512926115106248574422619, -01); | |
664 | ||
665 | // This is equation (entry) 31 from | |
666 | // http://www-2.cs.cmu.edu/~adamchik/articles/catalan/catalan.htm | |
667 | // See also http://www.mpfr.org/algorithms.pdf | |
668 | BOOST_MATH_STD_USING | |
669 | T k_fact = 1; | |
670 | T tk_fact = 1; | |
671 | T sum = 1; | |
672 | T term; | |
673 | T lim = N ? ldexp(T(1), 1 - (std::min)(N, tools::digits<T>())) : tools::epsilon<T>(); | |
674 | ||
675 | for(unsigned k = 1;; ++k) | |
676 | { | |
677 | k_fact *= k; | |
678 | tk_fact *= (2 * k) * (2 * k - 1); | |
679 | term = k_fact * k_fact / (tk_fact * (2 * k + 1) * (2 * k + 1)); | |
680 | sum += term; | |
681 | if(term < lim) | |
682 | { | |
683 | break; | |
684 | } | |
685 | } | |
686 | return boost::math::constants::pi<T, boost::math::policies::policy<> >() | |
687 | * log(2 + boost::math::constants::root_three<T, boost::math::policies::policy<> >()) | |
688 | / 8 | |
689 | + 3 * sum / 8; | |
690 | } | |
691 | ||
692 | namespace khinchin_detail{ | |
693 | ||
694 | template <class T> | |
695 | T zeta_polynomial_series(T s, T sc, int digits) | |
696 | { | |
697 | BOOST_MATH_STD_USING | |
698 | // | |
699 | // This is algorithm 3 from: | |
700 | // | |
701 | // "An Efficient Algorithm for the Riemann Zeta Function", P. Borwein, | |
702 | // Canadian Mathematical Society, Conference Proceedings, 2000. | |
703 | // See: http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf | |
704 | // | |
705 | BOOST_MATH_STD_USING | |
706 | int n = (digits * 19) / 53; | |
707 | T sum = 0; | |
708 | T two_n = ldexp(T(1), n); | |
709 | int ej_sign = 1; | |
710 | for(int j = 0; j < n; ++j) | |
711 | { | |
712 | sum += ej_sign * -two_n / pow(T(j + 1), s); | |
713 | ej_sign = -ej_sign; | |
714 | } | |
715 | T ej_sum = 1; | |
716 | T ej_term = 1; | |
717 | for(int j = n; j <= 2 * n - 1; ++j) | |
718 | { | |
719 | sum += ej_sign * (ej_sum - two_n) / pow(T(j + 1), s); | |
720 | ej_sign = -ej_sign; | |
721 | ej_term *= 2 * n - j; | |
722 | ej_term /= j - n + 1; | |
723 | ej_sum += ej_term; | |
724 | } | |
725 | return -sum / (two_n * (1 - pow(T(2), sc))); | |
726 | } | |
727 | ||
728 | template <class T> | |
729 | T khinchin(int digits) | |
730 | { | |
731 | BOOST_MATH_STD_USING | |
732 | T sum = 0; | |
733 | T term; | |
734 | T lim = ldexp(T(1), 1-digits); | |
735 | T factor = 0; | |
736 | unsigned last_k = 1; | |
737 | T num = 1; | |
738 | for(unsigned n = 1;; ++n) | |
739 | { | |
740 | for(unsigned k = last_k; k <= 2 * n - 1; ++k) | |
741 | { | |
742 | factor += num / k; | |
743 | num = -num; | |
744 | } | |
745 | last_k = 2 * n; | |
746 | term = (zeta_polynomial_series(T(2 * n), T(1 - T(2 * n)), digits) - 1) * factor / n; | |
747 | sum += term; | |
748 | if(term < lim) | |
749 | break; | |
750 | } | |
751 | return exp(sum / boost::math::constants::ln_two<T, boost::math::policies::policy<> >()); | |
752 | } | |
753 | ||
754 | } | |
755 | ||
756 | template <class T> | |
757 | template<int N> | |
1e59de90 | 758 | inline T constant_khinchin<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
759 | { |
760 | int n = N ? (std::min)(N, tools::digits<T>()) : tools::digits<T>(); | |
761 | return khinchin_detail::khinchin<T>(n); | |
762 | } | |
763 | ||
764 | template <class T> | |
765 | template<int N> | |
1e59de90 TL |
766 | inline T constant_extreme_value_skewness<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
767 | { // N[12 Sqrt[6] Zeta[3]/Pi^3, 1101] | |
7c673cae FG |
768 | BOOST_MATH_STD_USING |
769 | T ev(12 * sqrt(static_cast<T>(6)) * zeta_three<T, policies::policy<policies::digits2<N> > >() | |
770 | / pi_cubed<T, policies::policy<policies::digits2<N> > >() ); | |
771 | ||
772 | //T ev( | |
773 | //"1.1395470994046486574927930193898461120875997958365518247216557100852480077060706857071875468869385150" | |
774 | //"1894272048688553376986765366075828644841024041679714157616857834895702411080704529137366329462558680" | |
775 | //"2015498788776135705587959418756809080074611906006528647805347822929577145038743873949415294942796280" | |
776 | //"0895597703063466053535550338267721294164578901640163603544404938283861127819804918174973533694090594" | |
777 | //"3094963822672055237678432023017824416203652657301470473548274848068762500300316769691474974950757965" | |
778 | //"8640779777748741897542093874605477776538884083378029488863880220988107155275203245233994097178778984" | |
779 | //"3488995668362387892097897322246698071290011857605809901090220903955815127463328974447572119951192970" | |
780 | //"3684453635456559086126406960279692862247058250100678008419431185138019869693206366891639436908462809" | |
781 | //"9756051372711251054914491837034685476095423926553367264355374652153595857163724698198860485357368964" | |
782 | //"3807049634423621246870868566707915720704996296083373077647528285782964567312903914752617978405994377" | |
783 | //"9064157147206717895272199736902453130842229559980076472936976287378945035706933650987259357729800315"); | |
784 | ||
785 | return ev; | |
786 | } | |
787 | ||
788 | namespace detail{ | |
789 | // | |
790 | // Calculation of the Glaisher constant depends upon calculating the | |
791 | // derivative of the zeta function at 2, we can then use the relation: | |
792 | // zeta'(2) = 1/6 pi^2 [euler + ln(2pi)-12ln(A)] | |
793 | // To get the constant A. | |
794 | // See equation 45 at http://mathworld.wolfram.com/RiemannZetaFunction.html. | |
795 | // | |
796 | // The derivative of the zeta function is computed by direct differentiation | |
797 | // of the relation: | |
798 | // (1-2^(1-s))zeta(s) = SUM(n=0, INF){ (-n)^n / (n+1)^s } | |
799 | // Which gives us 2 slowly converging but alternating sums to compute, | |
800 | // for this we use Algorithm 1 from "Convergent Acceleration of Alternating Series", | |
801 | // Henri Cohen, Fernando Rodriguez Villegas and Don Zagier, Experimental Mathematics 9:1 (1999). | |
802 | // See http://www.math.utexas.edu/users/villegas/publications/conv-accel.pdf | |
803 | // | |
804 | template <class T> | |
805 | T zeta_series_derivative_2(unsigned digits) | |
806 | { | |
807 | // Derivative of the series part, evaluated at 2: | |
808 | BOOST_MATH_STD_USING | |
809 | int n = digits * 301 * 13 / 10000; | |
7c673cae FG |
810 | T d = pow(3 + sqrt(T(8)), n); |
811 | d = (d + 1 / d) / 2; | |
812 | T b = -1; | |
813 | T c = -d; | |
814 | T s = 0; | |
815 | for(int k = 0; k < n; ++k) | |
816 | { | |
817 | T a = -log(T(k+1)) / ((k+1) * (k+1)); | |
818 | c = b - c; | |
819 | s = s + c * a; | |
820 | b = (k + n) * (k - n) * b / ((k + T(0.5f)) * (k + 1)); | |
821 | } | |
822 | return s / d; | |
823 | } | |
824 | ||
825 | template <class T> | |
826 | T zeta_series_2(unsigned digits) | |
827 | { | |
828 | // Series part of zeta at 2: | |
829 | BOOST_MATH_STD_USING | |
830 | int n = digits * 301 * 13 / 10000; | |
831 | T d = pow(3 + sqrt(T(8)), n); | |
832 | d = (d + 1 / d) / 2; | |
833 | T b = -1; | |
834 | T c = -d; | |
835 | T s = 0; | |
836 | for(int k = 0; k < n; ++k) | |
837 | { | |
838 | T a = T(1) / ((k + 1) * (k + 1)); | |
839 | c = b - c; | |
840 | s = s + c * a; | |
841 | b = (k + n) * (k - n) * b / ((k + T(0.5f)) * (k + 1)); | |
842 | } | |
843 | return s / d; | |
844 | } | |
845 | ||
846 | template <class T> | |
847 | inline T zeta_series_lead_2() | |
848 | { | |
849 | // lead part at 2: | |
850 | return 2; | |
851 | } | |
852 | ||
853 | template <class T> | |
854 | inline T zeta_series_derivative_lead_2() | |
855 | { | |
856 | // derivative of lead part at 2: | |
857 | return -2 * boost::math::constants::ln_two<T>(); | |
858 | } | |
859 | ||
860 | template <class T> | |
861 | inline T zeta_derivative_2(unsigned n) | |
862 | { | |
863 | // zeta derivative at 2: | |
864 | return zeta_series_derivative_2<T>(n) * zeta_series_lead_2<T>() | |
865 | + zeta_series_derivative_lead_2<T>() * zeta_series_2<T>(n); | |
866 | } | |
867 | ||
868 | } // namespace detail | |
869 | ||
870 | template <class T> | |
871 | template<int N> | |
1e59de90 | 872 | inline T constant_glaisher<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
873 | { |
874 | ||
875 | BOOST_MATH_STD_USING | |
876 | typedef policies::policy<policies::digits2<N> > forwarding_policy; | |
877 | int n = N ? (std::min)(N, tools::digits<T>()) : tools::digits<T>(); | |
878 | T v = detail::zeta_derivative_2<T>(n); | |
879 | v *= 6; | |
880 | v /= boost::math::constants::pi<T, forwarding_policy>() * boost::math::constants::pi<T, forwarding_policy>(); | |
881 | v -= boost::math::constants::euler<T, forwarding_policy>(); | |
882 | v -= log(2 * boost::math::constants::pi<T, forwarding_policy>()); | |
883 | v /= -12; | |
884 | return exp(v); | |
885 | ||
1e59de90 | 886 | /* |
7c673cae FG |
887 | // from http://mpmath.googlecode.com/svn/data/glaisher.txt |
888 | // 20,000 digits of the Glaisher-Kinkelin constant A = exp(1/2 - zeta'(-1)) | |
889 | // Computed using A = exp((6 (-zeta'(2))/pi^2 + log 2 pi + gamma)/12) | |
1e59de90 TL |
890 | // with Euler-Maclaurin summation for zeta'(2). |
891 | T g( | |
892 | "1.282427129100622636875342568869791727767688927325001192063740021740406308858826" | |
893 | "46112973649195820237439420646120399000748933157791362775280404159072573861727522" | |
894 | "14334327143439787335067915257366856907876561146686449997784962754518174312394652" | |
895 | "76128213808180219264516851546143919901083573730703504903888123418813674978133050" | |
896 | "93770833682222494115874837348064399978830070125567001286994157705432053927585405" | |
897 | "81731588155481762970384743250467775147374600031616023046613296342991558095879293" | |
898 | "36343887288701988953460725233184702489001091776941712153569193674967261270398013" | |
899 | "52652668868978218897401729375840750167472114895288815996668743164513890306962645" | |
900 | "59870469543740253099606800842447417554061490189444139386196089129682173528798629" | |
901 | "88434220366989900606980888785849587494085307347117090132667567503310523405221054" | |
902 | "14176776156308191919997185237047761312315374135304725819814797451761027540834943" | |
903 | "14384965234139453373065832325673954957601692256427736926358821692159870775858274" | |
904 | "69575162841550648585890834128227556209547002918593263079373376942077522290940187"); | |
905 | ||
906 | return g; | |
907 | */ | |
7c673cae FG |
908 | } |
909 | ||
910 | template <class T> | |
911 | template<int N> | |
1e59de90 TL |
912 | inline T constant_rayleigh_skewness<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
913 | { // 1100 digits of the Rayleigh distribution skewness | |
914 | // N[2 Sqrt[Pi] (Pi - 3)/((4 - Pi)^(3/2)), 1100] | |
7c673cae FG |
915 | |
916 | BOOST_MATH_STD_USING | |
917 | T rs(2 * root_pi<T, policies::policy<policies::digits2<N> > >() | |
918 | * pi_minus_three<T, policies::policy<policies::digits2<N> > >() | |
919 | / pow(four_minus_pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(3./2)) | |
920 | ); | |
1e59de90 | 921 | // 6.31110657818937138191899351544227779844042203134719497658094585692926819617473725459905027032537306794400047264, |
7c673cae | 922 | |
1e59de90 TL |
923 | //"0.6311106578189371381918993515442277798440422031347194976580945856929268196174737254599050270325373067" |
924 | //"9440004726436754739597525250317640394102954301685809920213808351450851396781817932734836994829371322" | |
925 | //"5797376021347531983451654130317032832308462278373358624120822253764532674177325950686466133508511968" | |
926 | //"2389168716630349407238090652663422922072397393006683401992961569208109477307776249225072042971818671" | |
927 | //"4058887072693437217879039875871765635655476241624825389439481561152126886932506682176611183750503553" | |
928 | //"1218982627032068396407180216351425758181396562859085306247387212297187006230007438534686340210168288" | |
929 | //"8956816965453815849613622117088096547521391672977226658826566757207615552041767516828171274858145957" | |
930 | //"6137539156656005855905288420585194082284972984285863898582313048515484073396332610565441264220790791" | |
931 | //"0194897267890422924599776483890102027823328602965235306539844007677157873140562950510028206251529523" | |
932 | //"7428049693650605954398446899724157486062545281504433364675815915402937209673727753199567661561209251" | |
933 | //"4695589950526053470201635372590001578503476490223746511106018091907936826431407434894024396366284848"); ; | |
934 | return rs; | |
7c673cae FG |
935 | } |
936 | ||
937 | template <class T> | |
938 | template<int N> | |
1e59de90 | 939 | inline T constant_rayleigh_kurtosis_excess<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
940 | { // - (6 Pi^2 - 24 Pi + 16)/((Pi - 4)^2) |
941 | // Might provide and calculate this using pi_minus_four. | |
942 | BOOST_MATH_STD_USING | |
943 | return - (((static_cast<T>(6) * pi<T, policies::policy<policies::digits2<N> > >() | |
944 | * pi<T, policies::policy<policies::digits2<N> > >()) | |
945 | - (static_cast<T>(24) * pi<T, policies::policy<policies::digits2<N> > >()) + static_cast<T>(16) ) | |
946 | / | |
947 | ((pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4)) | |
948 | * (pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4))) | |
949 | ); | |
950 | } | |
951 | ||
952 | template <class T> | |
953 | template<int N> | |
1e59de90 | 954 | inline T constant_rayleigh_kurtosis<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
7c673cae FG |
955 | { // 3 - (6 Pi^2 - 24 Pi + 16)/((Pi - 4)^2) |
956 | // Might provide and calculate this using pi_minus_four. | |
957 | BOOST_MATH_STD_USING | |
958 | return static_cast<T>(3) - (((static_cast<T>(6) * pi<T, policies::policy<policies::digits2<N> > >() | |
959 | * pi<T, policies::policy<policies::digits2<N> > >()) | |
960 | - (static_cast<T>(24) * pi<T, policies::policy<policies::digits2<N> > >()) + static_cast<T>(16) ) | |
961 | / | |
962 | ((pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4)) | |
963 | * (pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4))) | |
964 | ); | |
965 | } | |
966 | ||
92f5a8d4 TL |
967 | template <class T> |
968 | template<int N> | |
1e59de90 | 969 | inline T constant_log2_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
92f5a8d4 TL |
970 | { |
971 | return 1 / boost::math::constants::ln_two<T>(); | |
972 | } | |
973 | ||
974 | template <class T> | |
975 | template<int N> | |
1e59de90 | 976 | inline T constant_quarter_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
92f5a8d4 TL |
977 | { |
978 | return boost::math::constants::pi<T>() / 4; | |
979 | } | |
980 | ||
981 | template <class T> | |
982 | template<int N> | |
1e59de90 | 983 | inline T constant_one_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
92f5a8d4 TL |
984 | { |
985 | return 1 / boost::math::constants::pi<T>(); | |
986 | } | |
987 | ||
988 | template <class T> | |
989 | template<int N> | |
1e59de90 | 990 | inline T constant_two_div_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
92f5a8d4 TL |
991 | { |
992 | return 2 * boost::math::constants::one_div_root_pi<T>(); | |
993 | } | |
994 | ||
20effc67 TL |
995 | #if __cplusplus >= 201103L || (defined(_MSC_VER) && _MSC_VER >= 1900) |
996 | template <class T> | |
997 | template<int N> | |
1e59de90 | 998 | inline T constant_first_feigenbaum<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
20effc67 TL |
999 | { |
1000 | // We know the constant to 1018 decimal digits. | |
1001 | // See: http://www.plouffe.fr/simon/constants/feigenbaum.txt | |
1002 | // Also: https://oeis.org/A006890 | |
1003 | // N is in binary digits; so we multiply by log_2(10) | |
1004 | ||
1e59de90 | 1005 | static_assert(N < 3.321*1018, "\nThe first Feigenbaum constant cannot be computed at runtime; it is too expensive. It is known to 1018 decimal digits; you must request less than that."); |
20effc67 TL |
1006 | T alpha{"4.6692016091029906718532038204662016172581855774757686327456513430041343302113147371386897440239480138171659848551898151344086271420279325223124429888908908599449354632367134115324817142199474556443658237932020095610583305754586176522220703854106467494942849814533917262005687556659523398756038256372256480040951071283890611844702775854285419801113440175002428585382498335715522052236087250291678860362674527213399057131606875345083433934446103706309452019115876972432273589838903794946257251289097948986768334611626889116563123474460575179539122045562472807095202198199094558581946136877445617396074115614074243754435499204869180982648652368438702799649017397793425134723808737136211601860128186102056381818354097598477964173900328936171432159878240789776614391395764037760537119096932066998361984288981837003229412030210655743295550388845849737034727532121925706958414074661841981961006129640161487712944415901405467941800198133253378592493365883070459999938375411726563553016862529032210862320550634510679399023341675"}; |
1007 | return alpha; | |
1008 | } | |
1009 | ||
1010 | template <class T> | |
1011 | template<int N> | |
1e59de90 | 1012 | inline T constant_plastic<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
20effc67 TL |
1013 | { |
1014 | using std::cbrt; | |
1015 | using std::sqrt; | |
1016 | return (cbrt(9-sqrt(T(69))) + cbrt(9+sqrt(T(69))))/cbrt(T(18)); | |
1017 | } | |
1018 | ||
1019 | ||
1020 | template <class T> | |
1021 | template<int N> | |
1e59de90 | 1022 | inline T constant_gauss<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
20effc67 TL |
1023 | { |
1024 | using std::sqrt; | |
1025 | T a = sqrt(T(2)); | |
1026 | T g = 1; | |
1027 | const T scale = sqrt(std::numeric_limits<T>::epsilon())/512; | |
1028 | while (a-g > scale*g) | |
1029 | { | |
1030 | T anp1 = (a + g)/2; | |
1031 | g = sqrt(a*g); | |
1032 | a = anp1; | |
1033 | } | |
1034 | ||
1035 | return 2/(a + g); | |
1036 | } | |
1037 | ||
1038 | template <class T> | |
1039 | template<int N> | |
1e59de90 | 1040 | inline T constant_dottie<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
20effc67 TL |
1041 | { |
1042 | // Error analysis: cos(x(1+d)) - x(1+d) = -(sin(x)+1)xd; plug in x = 0.739 gives -1.236d; take d as half an ulp gives the termination criteria we want. | |
1043 | using std::cos; | |
1044 | using std::abs; | |
1045 | using std::sin; | |
1046 | T x{".739085133215160641655312087673873404013411758900757464965680635773284654883547594599376106931766531849801246"}; | |
1047 | T residual = cos(x) - x; | |
1048 | do { | |
1049 | x += residual/(sin(x)+1); | |
1050 | residual = cos(x) - x; | |
1051 | } while(abs(residual) > std::numeric_limits<T>::epsilon()); | |
1052 | return x; | |
1053 | } | |
1054 | ||
1055 | ||
1056 | template <class T> | |
1057 | template<int N> | |
1e59de90 | 1058 | inline T constant_reciprocal_fibonacci<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
20effc67 TL |
1059 | { |
1060 | // Wikipedia says Gosper has deviced a faster algorithm for this, but I read the linked paper and couldn't see it! | |
1061 | // In any case, k bits per iteration is fine, though it would be better to sum from smallest to largest. | |
1062 | // That said, the condition number is unity, so it should be fine. | |
1063 | T x0 = 1; | |
1064 | T x1 = 1; | |
1065 | T sum = 2; | |
1066 | T diff = 1; | |
1067 | while (diff > std::numeric_limits<T>::epsilon()) { | |
1068 | T tmp = x1 + x0; | |
1069 | diff = 1/tmp; | |
1070 | sum += diff; | |
1071 | x0 = x1; | |
1072 | x1 = tmp; | |
1073 | } | |
1074 | return sum; | |
1075 | } | |
1076 | ||
1077 | template <class T> | |
1078 | template<int N> | |
1e59de90 | 1079 | inline T constant_laplace_limit<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC((std::integral_constant<int, N>))) |
20effc67 TL |
1080 | { |
1081 | // If x is the exact root, then the approximate root is given by x(1+delta). | |
1082 | // Plugging this into the equation for the Laplace limit gives the residual of approximately | |
1083 | // 2.6389delta. Take delta as half an epsilon and give some leeway so we don't get caught in an infinite loop, | |
1084 | // gives a termination condition as 2eps. | |
1085 | using std::abs; | |
1086 | using std::exp; | |
1087 | using std::sqrt; | |
1088 | T x{"0.66274341934918158097474209710925290705623354911502241752039253499097185308651127724965480259895818168"}; | |
1089 | T tmp = sqrt(1+x*x); | |
1090 | T etmp = exp(tmp); | |
1091 | T residual = x*exp(tmp) - 1 - tmp; | |
1092 | T df = etmp -x/tmp + etmp*x*x/tmp; | |
1093 | do { | |
1094 | x -= residual/df; | |
1095 | tmp = sqrt(1+x*x); | |
1096 | etmp = exp(tmp); | |
1097 | residual = x*exp(tmp) - 1 - tmp; | |
1098 | df = etmp -x/tmp + etmp*x*x/tmp; | |
1099 | } while(abs(residual) > 2*std::numeric_limits<T>::epsilon()); | |
1100 | return x; | |
1101 | } | |
1102 | ||
1103 | #endif | |
1104 | ||
92f5a8d4 TL |
1105 | } |
1106 | } | |
1107 | } | |
1108 | } // namespaces | |
7c673cae FG |
1109 | |
1110 | #endif // BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED |