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

source code of include/boost/math/constants/calculate_constants.hpp