2 ///////////////////////////////////////////////////////////////////////////////
3 // Copyright Christopher Kormanyos 2013 - 2014.
4 // Copyright John Maddock 2013.
5 // Distributed under the Boost Software License,
6 // Version 1.0. (See accompanying file LICENSE_1_0.txt
7 // or copy at http://www.boost.org/LICENSE_1_0.txt)
9 // This work is based on an earlier work:
10 // "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
11 // in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
22 #include <boost/math/constants/constants.hpp>
23 #include <boost/noncopyable.hpp>
25 //#define USE_CPP_BIN_FLOAT
26 #define USE_CPP_DEC_FLOAT
29 #if !defined(DIGIT_COUNT)
30 #define DIGIT_COUNT 100
33 #if !defined(BOOST_NO_CXX11_HDR_CHRONO)
35 #define STD_CHRONO std::chrono
37 #include <boost/chrono.hpp>
38 #define STD_CHRONO boost::chrono
41 #if defined(USE_CPP_BIN_FLOAT)
42 #include <boost/multiprecision/cpp_bin_float.hpp>
43 typedef boost::multiprecision::number
<boost::multiprecision::cpp_bin_float
<DIGIT_COUNT
+ 10> > mp_type
;
44 #elif defined(USE_CPP_DEC_FLOAT)
45 #include <boost/multiprecision/cpp_dec_float.hpp>
46 typedef boost::multiprecision::number
<boost::multiprecision::cpp_dec_float
<DIGIT_COUNT
+ 10> > mp_type
;
47 #elif defined(USE_MPFR)
48 #include <boost/multiprecision/mpfr.hpp>
49 typedef boost::multiprecision::number
<boost::multiprecision::mpfr_float_backend
<DIGIT_COUNT
+ 10> > mp_type
;
51 #error no multiprecision floating type is defined
54 template <class clock_type
>
58 typedef typename
clock_type::duration duration_type
;
60 stopwatch() : m_start(clock_type::now()) { }
62 stopwatch(const stopwatch
& other
) : m_start(other
.m_start
) { }
64 stopwatch
& operator=(const stopwatch
& other
)
66 m_start
= other
.m_start
;
72 duration_type
elapsed() const
74 return (clock_type::now() - m_start
);
79 m_start
= clock_type::now();
83 typename
clock_type::time_point m_start
;
88 template<class T
> T
chebyshev_t(const std::int32_t n
, const T
& x
);
90 template<class T
> T
chebyshev_t(const std::uint32_t n
, const T
& x
, std::vector
<T
>* vp
);
92 template<class T
> bool isneg(const T
& x
) { return (x
< T(0)); }
94 template<class T
> const T
& zero() { static const T
value_zero(0); return value_zero
; }
95 template<class T
> const T
& one () { static const T
value_one (1); return value_one
; }
96 template<class T
> const T
& two () { static const T
value_two (2); return value_two
; }
99 namespace orthogonal_polynomial_series
101 template<typename T
> static inline T
orthogonal_polynomial_template(const T
& x
, const std::uint32_t n
, std::vector
<T
>* const vp
= static_cast<std::vector
<T
>*>(0u))
103 // Compute the value of an orthogonal chebyshev polinomial.
104 // Use stable upward recursion.
109 vp
->reserve(static_cast<std::size_t>(n
+ 1u));
112 T y0
= my_math::one
<T
>();
114 if(vp
!= nullptr) { vp
->push_back(y0
); }
116 if(n
== static_cast<std::uint32_t>(0u))
123 if(vp
!= nullptr) { vp
->push_back(y1
); }
125 if(n
== static_cast<std::uint32_t>(1u))
130 T a
= my_math::two
<T
>();
131 T b
= my_math::zero
<T
>();
132 T c
= my_math::one
<T
>();
136 // Calculate higher orders using the recurrence relation.
137 // The direction of stability is upward recursion.
138 for(std::int32_t k
= static_cast<std::int32_t>(2); k
<= static_cast<std::int32_t>(n
); ++k
)
140 yk
= (((a
* x
) + b
) * y1
) - (c
* y0
);
145 if(vp
!= nullptr) { vp
->push_back(yk
); }
152 template<class T
> T
my_math::chebyshev_t(const std::int32_t n
, const T
& x
)
154 if(my_math::isneg(x
))
156 const bool b_negate
= ((n
% static_cast<std::int32_t>(2)) != static_cast<std::int32_t>(0));
158 const T y
= chebyshev_t(n
, -x
);
160 return (!b_negate
? y
: -y
);
163 if(n
< static_cast<std::int32_t>(0))
165 const std::int32_t nn
= static_cast<std::int32_t>(-n
);
167 return chebyshev_t(nn
, x
);
171 return orthogonal_polynomial_series::orthogonal_polynomial_template(x
, static_cast<std::uint32_t>(n
));
175 template<class T
> T
my_math::chebyshev_t(const std::uint32_t n
, const T
& x
, std::vector
<T
>* const vp
) { return orthogonal_polynomial_series::orthogonal_polynomial_template(x
, static_cast<std::int32_t>(n
), vp
); }
179 template <class T
> float digit_scale()
181 const int d
= ((std::max
)(std::numeric_limits
<T
>::digits10
, 15));
182 return static_cast<float>(d
) / 300.0F
;
190 template<typename T
> class hypergeometric_pfq_base
: private boost::noncopyable
193 virtual ~hypergeometric_pfq_base() { }
195 virtual void ccoef() const = 0;
197 virtual T
series() const
199 using my_math::chebyshev_t
;
201 // Compute the Chebyshev coefficients.
202 // Get the values of the shifted Chebyshev polynomials.
203 std::vector
<T
> chebyshev_t_shifted_values
;
204 const T z_shifted
= ((Z
/ W
) * static_cast<std::int32_t>(2)) - static_cast<std::int32_t>(1);
206 chebyshev_t(static_cast<std::uint32_t>(C
.size()),
208 &chebyshev_t_shifted_values
);
210 // Luke: C ---------- COMPUTE SCALE FACTOR ----------
212 // Luke: C ---------- SCALE THE COEFFICIENTS ----------
215 // The coefficient scaling is preformed after the Chebyshev summation,
216 // and it is carried out with a single division operation.
219 const T scale
= std::accumulate(C
.begin(),
222 [&b_neg
](T
& scale_sum
, const T
& ck
) -> T
224 ((!b_neg
) ? (scale_sum
+= ck
) : (scale_sum
-= ck
));
229 // Compute the result of the series expansion using unscaled coefficients.
230 const T sum
= std::inner_product(C
.begin(),
232 chebyshev_t_shifted_values
.begin(),
235 // Return the properly scaled result.
242 mutable std::deque
<T
> C
;
244 hypergeometric_pfq_base(const T
& z
,
249 virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale
<T
>() * 500.0F
); }
252 template<typename T
> class ccoef4_hypergeometric_0f1
: public hypergeometric_pfq_base
<T
>
255 ccoef4_hypergeometric_0f1(const T
& c
,
257 const T
& w
) : hypergeometric_pfq_base
<T
>(z
, w
),
260 virtual ~ccoef4_hypergeometric_0f1() { }
262 virtual void ccoef() const
264 // See Luke 1977 page 80.
265 const std::int32_t N1
= static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1));
266 const std::int32_t N2
= static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(2));
268 // Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
269 // Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
273 T
A1(boost::math::tools::root_epsilon
<T
>());
275 hypergeometric_pfq_base
<T
>::C
.resize(1u, A1
);
277 std::int32_t X1
= N2
;
281 const T Z1
= T(4) / hypergeometric_pfq_base
<T
>::W
;
283 for(std::int32_t k
= static_cast<std::int32_t>(0); k
< N1
; ++k
)
285 const T DIVFAC
= T(1) / X1
;
289 // The terms have been slightly re-arranged resulting in lower complexity.
290 // Parentheses have been added to avoid reliance on operator precedence.
291 const T term
= (A2
- ((A3
* DIVFAC
) * X1
))
292 + ((A2
* X1
) * ((1 + (C1
+ X1
)) * Z1
))
293 + ((A1
* X1
) * ((DIVFAC
- (C1
* Z1
)) + (X1
* Z1
)));
295 hypergeometric_pfq_base
<T
>::C
.push_front(term
);
299 A1
= hypergeometric_pfq_base
<T
>::C
.front();
302 hypergeometric_pfq_base
<T
>::C
.front() /= static_cast<std::int32_t>(2);
309 template<typename T
> class ccoef1_hypergeometric_1f0
: public hypergeometric_pfq_base
<T
>
312 ccoef1_hypergeometric_1f0(const T
& a
,
314 const T
& w
) : hypergeometric_pfq_base
<T
>(z
, w
),
317 virtual ~ccoef1_hypergeometric_1f0() { }
319 virtual void ccoef() const
321 // See Luke 1977 page 67.
322 const std::int32_t N1
= static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
323 const std::int32_t N2
= static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2));
325 // Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
326 // Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
329 T
A1(boost::math::tools::root_epsilon
<T
>());
331 hypergeometric_pfq_base
<T
>::C
.resize(1u, A1
);
333 std::int32_t X1
= N2
;
337 // Here, we have corrected what appears to be an error in Luke's code.
339 // Luke's original code listing has:
341 // But it appears as though the correct form is:
342 // AFAC = 2 - FOUR/W.
344 const T AFAC
= 2 - (T(4) / hypergeometric_pfq_base
<T
>::W
);
346 for(std::int32_t k
= static_cast<std::int32_t>(0); k
< N1
; ++k
)
350 // The terms have been slightly re-arranged resulting in lower complexity.
351 // Parentheses have been added to avoid reliance on operator precedence.
352 const T term
= -(((X1
* AFAC
) * A1
) + ((X1
+ V1
) * A2
)) / (X1
- V1
);
354 hypergeometric_pfq_base
<T
>::C
.push_front(term
);
357 A1
= hypergeometric_pfq_base
<T
>::C
.front();
360 hypergeometric_pfq_base
<T
>::C
.front() /= static_cast<std::int32_t>(2);
366 virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale
<T
>() * 1600.0F
); }
369 template<typename T
> class ccoef3_hypergeometric_1f1
: public hypergeometric_pfq_base
<T
>
372 ccoef3_hypergeometric_1f1(const T
& a
,
375 const T
& w
) : hypergeometric_pfq_base
<T
>(z
, w
),
379 virtual ~ccoef3_hypergeometric_1f1() { }
381 virtual void ccoef() const
383 // See Luke 1977 page 74.
384 const std::int32_t N1
= static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1));
385 const std::int32_t N2
= static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(2));
387 // Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
388 // Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
392 T
A1(boost::math::tools::root_epsilon
<T
>());
394 hypergeometric_pfq_base
<T
>::C
.resize(1u, A1
);
397 std::int32_t X1
= N2
;
400 T X3A
= (X
+ 3) - AP
;
402 const T Z1
= T(4) / hypergeometric_pfq_base
<T
>::W
;
404 for(std::int32_t k
= static_cast<std::int32_t>(0); k
< N1
; ++k
)
411 const T X3A_over_X2
= X3A
/ static_cast<std::int32_t>(X
+ 2);
413 // The terms have been slightly re-arranged resulting in lower complexity.
414 // Parentheses have been added to avoid reliance on operator precedence.
415 const T PART1
= A1
* (((X
+ CP
) * Z1
) - X3A_over_X2
);
416 const T PART2
= A2
* (Z1
* ((X
+ 3) - CP
) + (XA
/ X1
));
417 const T PART3
= A3
* X3A_over_X2
;
419 const T term
= (((PART1
+ PART2
) + PART3
) * X1
) / XA
;
421 hypergeometric_pfq_base
<T
>::C
.push_front(term
);
425 A1
= hypergeometric_pfq_base
<T
>::C
.front();
428 hypergeometric_pfq_base
<T
>::C
.front() /= static_cast<std::int32_t>(2);
436 template<typename T
> class ccoef6_hypergeometric_1f2
: public hypergeometric_pfq_base
<T
>
439 ccoef6_hypergeometric_1f2(const T
& a
,
443 const T
& w
) : hypergeometric_pfq_base
<T
>(z
, w
),
448 virtual ~ccoef6_hypergeometric_1f2() { }
450 virtual void ccoef() const
452 // See Luke 1977 page 85.
453 const std::int32_t N1
= static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1));
455 // Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
456 // Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
461 T
A1(boost::math::tools::root_epsilon
<T
>());
463 hypergeometric_pfq_base
<T
>::C
.resize(1u, A1
);
468 const T Z1
= T(4) / hypergeometric_pfq_base
<T
>::W
;
470 for(std::int32_t k
= static_cast<std::int32_t>(0); k
< N1
; ++k
)
475 const std::int32_t TWO_X
= static_cast<std::int32_t>(X
* 2);
476 const std::int32_t X_PLUS_1
= static_cast<std::int32_t>(X
+ 1);
477 const std::int32_t X_PLUS_3
= static_cast<std::int32_t>(X
+ 3);
478 const std::int32_t X_PLUS_4
= static_cast<std::int32_t>(X
+ 4);
480 const T QQ
= T(TWO_X
+ 3) / static_cast<std::int32_t>(TWO_X
+ static_cast<std::int32_t>(5));
481 const T SS
= (X
+ BP
) * (X
+ CP
);
483 // The terms have been slightly re-arranged resulting in lower complexity.
484 // Parentheses have been added to avoid reliance on operator precedence.
485 const T PART1
= A1
* (((PP
- (QQ
* (PP
+ 1))) * 2) + (SS
* Z1
));
486 const T PART2
= (A2
* (X
+ 2)) * ((((TWO_X
+ 1) * PP
) / X_PLUS_1
) - ((QQ
* 4) * (PP
+ 1)) + (((TWO_X
+ 3) * (PP
+ 2)) / X_PLUS_3
) + ((Z1
* 2) * (SS
- (QQ
* (X_PLUS_1
+ BP
)) * (X_PLUS_1
+ CP
))));
487 const T PART3
= A3
* ((((X_PLUS_3
- AP
) - (QQ
* (X_PLUS_4
- AP
))) * 2) + (((QQ
* Z1
) * (X_PLUS_4
- BP
)) * (X_PLUS_4
- CP
)));
488 const T PART4
= ((A4
* QQ
) * (X_PLUS_4
- AP
)) / X_PLUS_3
;
490 const T term
= (((PART1
- PART2
) + (PART3
- PART4
)) * X_PLUS_1
) / PP
;
492 hypergeometric_pfq_base
<T
>::C
.push_front(term
);
497 A1
= hypergeometric_pfq_base
<T
>::C
.front();
500 hypergeometric_pfq_base
<T
>::C
.front() /= static_cast<std::int32_t>(2);
509 template<typename T
> class ccoef2_hypergeometric_2f1
: public hypergeometric_pfq_base
<T
>
512 ccoef2_hypergeometric_2f1(const T
& a
,
516 const T
& w
) : hypergeometric_pfq_base
<T
>(z
, w
),
521 virtual ~ccoef2_hypergeometric_2f1() { }
523 virtual void ccoef() const
525 // See Luke 1977 page 59.
526 const std::int32_t N1
= static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
527 const std::int32_t N2
= static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2));
529 // Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
530 // Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
534 T
A1(boost::math::tools::root_epsilon
<T
>());
536 hypergeometric_pfq_base
<T
>::C
.resize(1u, A1
);
539 std::int32_t X1
= N2
;
540 std::int32_t X3
= static_cast<std::int32_t>((X
* 2) + 3);
542 T X3A
= (X
+ 3) - AP
;
543 T X3B
= (X
+ 3) - BP
;
545 const T Z1
= T(4) / hypergeometric_pfq_base
<T
>::W
;
547 for(std::int32_t k
= static_cast<std::int32_t>(0); k
< N1
; ++k
)
555 const std::int32_t X_PLUS_2
= static_cast<std::int32_t>(X
+ 2);
557 const T XAB
= T(1) / ((X
+ AP
) * (X
+ BP
));
559 // The terms have been slightly re-arranged resulting in lower complexity.
560 // Parentheses have been added to avoid reliance on operator precedence.
561 const T PART1
= (A1
* X1
) * (2 - (((AP
+ X1
) * (BP
+ X1
)) * ((T(X3
) / X_PLUS_2
) * XAB
)) + ((CP
+ X
) * (XAB
* Z1
)));
562 const T PART2
= (A2
* XAB
) * ((X3A
* X3B
) - (X3
* ((X3A
+ X3B
) - 1)) + (((3 - CP
) + X
) * (X1
* Z1
)));
563 const T PART3
= (A3
* X1
) * (X3A
/ X_PLUS_2
) * (X3B
* XAB
);
565 const T term
= (PART1
+ PART2
) - PART3
;
567 hypergeometric_pfq_base
<T
>::C
.push_front(term
);
571 A1
= hypergeometric_pfq_base
<T
>::C
.front();
574 hypergeometric_pfq_base
<T
>::C
.front() /= static_cast<std::int32_t>(2);
582 virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale
<T
>() * 1600.0F
); }
585 template<class T
> T
luke_ccoef4_hypergeometric_0f1(const T
& a
, const T
& x
);
586 template<class T
> T
luke_ccoef1_hypergeometric_1f0(const T
& a
, const T
& x
);
587 template<class T
> T
luke_ccoef3_hypergeometric_1f1(const T
& a
, const T
& b
, const T
& x
);
588 template<class T
> T
luke_ccoef6_hypergeometric_1f2(const T
& a
, const T
& b
, const T
& c
, const T
& x
);
589 template<class T
> T
luke_ccoef2_hypergeometric_2f1(const T
& a
, const T
& b
, const T
& c
, const T
& x
);
594 T
examples::nr_006::luke_ccoef4_hypergeometric_0f1(const T
& a
, const T
& x
)
596 const ccoef4_hypergeometric_0f1
<T
> hypergeometric_0f1_object(a
, x
, T(-20));
598 hypergeometric_0f1_object
.ccoef();
600 return hypergeometric_0f1_object
.series();
604 T
examples::nr_006::luke_ccoef1_hypergeometric_1f0(const T
& a
, const T
& x
)
606 const ccoef1_hypergeometric_1f0
<T
> hypergeometric_1f0_object(a
, x
, T(-20));
608 hypergeometric_1f0_object
.ccoef();
610 return hypergeometric_1f0_object
.series();
614 T
examples::nr_006::luke_ccoef3_hypergeometric_1f1(const T
& a
, const T
& b
, const T
& x
)
616 const ccoef3_hypergeometric_1f1
<T
> hypergeometric_1f1_object(a
, b
, x
, T(-20));
618 hypergeometric_1f1_object
.ccoef();
620 return hypergeometric_1f1_object
.series();
624 T
examples::nr_006::luke_ccoef6_hypergeometric_1f2(const T
& a
, const T
& b
, const T
& c
, const T
& x
)
626 const ccoef6_hypergeometric_1f2
<T
> hypergeometric_1f2_object(a
, b
, c
, x
, T(-20));
628 hypergeometric_1f2_object
.ccoef();
630 return hypergeometric_1f2_object
.series();
634 T
examples::nr_006::luke_ccoef2_hypergeometric_2f1(const T
& a
, const T
& b
, const T
& c
, const T
& x
)
636 const ccoef2_hypergeometric_2f1
<T
> hypergeometric_2f1_object(a
, b
, c
, x
, T(-20));
638 hypergeometric_2f1_object
.ccoef();
640 return hypergeometric_2f1_object
.series();
645 stopwatch
<STD_CHRONO::high_resolution_clock
> my_stopwatch
;
646 float total_time
= 0.0F
;
648 std::vector
<mp_type
> hypergeometric_0f1_results(20U);
649 std::vector
<mp_type
> hypergeometric_1f0_results(20U);
650 std::vector
<mp_type
> hypergeometric_1f1_results(20U);
651 std::vector
<mp_type
> hypergeometric_2f1_results(20U);
652 std::vector
<mp_type
> hypergeometric_1f2_results(20U);
654 const mp_type
a(mp_type(3) / 7);
655 const mp_type
b(mp_type(2) / 3);
656 const mp_type
c(mp_type(1) / 4);
658 std::int_least16_t i
;
660 std::cout
<< "test hypergeometric_0f1." << std::endl
;
662 my_stopwatch
.reset();
664 // Generate a table of values of Hypergeometric0F1.
665 // Compare with the Mathematica command:
666 // Table[N[HypergeometricPFQ[{}, {3/7}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
667 std::for_each(hypergeometric_0f1_results
.begin(),
668 hypergeometric_0f1_results
.end(),
669 [&i
, &a
](mp_type
& new_value
)
671 const mp_type
x(-(boost::math::constants::euler
<mp_type
>() * i
));
673 new_value
= examples::nr_006::luke_ccoef4_hypergeometric_0f1(a
, x
);
678 total_time
+= STD_CHRONO::duration_cast
<STD_CHRONO::duration
<float> >(my_stopwatch
.elapsed()).count();
680 // Print the values of Hypergeometric0F1.
681 std::for_each(hypergeometric_0f1_results
.begin(),
682 hypergeometric_0f1_results
.end(),
685 std::cout
<< std::setprecision(DIGIT_COUNT
) << h
<< std::endl
;
688 std::cout
<< "test hypergeometric_1f0." << std::endl
;
690 my_stopwatch
.reset();
692 // Generate a table of values of Hypergeometric1F0.
693 // Compare with the Mathematica command:
694 // Table[N[HypergeometricPFQ[{3/7}, {}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
695 std::for_each(hypergeometric_1f0_results
.begin(),
696 hypergeometric_1f0_results
.end(),
697 [&i
, &a
](mp_type
& new_value
)
699 const mp_type
x(-(boost::math::constants::euler
<mp_type
>() * i
));
701 new_value
= examples::nr_006::luke_ccoef1_hypergeometric_1f0(a
, x
);
706 total_time
+= STD_CHRONO::duration_cast
<STD_CHRONO::duration
<float> >(my_stopwatch
.elapsed()).count();
708 // Print the values of Hypergeometric1F0.
709 std::for_each(hypergeometric_1f0_results
.begin(),
710 hypergeometric_1f0_results
.end(),
713 std::cout
<< std::setprecision(DIGIT_COUNT
) << h
<< std::endl
;
716 std::cout
<< "test hypergeometric_1f1." << std::endl
;
718 my_stopwatch
.reset();
720 // Generate a table of values of Hypergeometric1F1.
721 // Compare with the Mathematica command:
722 // Table[N[HypergeometricPFQ[{3/7}, {2/3}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
723 std::for_each(hypergeometric_1f1_results
.begin(),
724 hypergeometric_1f1_results
.end(),
725 [&i
, &a
, &b
](mp_type
& new_value
)
727 const mp_type
x(-(boost::math::constants::euler
<mp_type
>() * i
));
729 new_value
= examples::nr_006::luke_ccoef3_hypergeometric_1f1(a
, b
, x
);
734 total_time
+= STD_CHRONO::duration_cast
<STD_CHRONO::duration
<float> >(my_stopwatch
.elapsed()).count();
736 // Print the values of Hypergeometric1F1.
737 std::for_each(hypergeometric_1f1_results
.begin(),
738 hypergeometric_1f1_results
.end(),
741 std::cout
<< std::setprecision(DIGIT_COUNT
) << h
<< std::endl
;
744 std::cout
<< "test hypergeometric_1f2." << std::endl
;
746 my_stopwatch
.reset();
748 // Generate a table of values of Hypergeometric1F2.
749 // Compare with the Mathematica command:
750 // Table[N[HypergeometricPFQ[{3/7}, {2/3, 1/4}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
751 std::for_each(hypergeometric_1f2_results
.begin(),
752 hypergeometric_1f2_results
.end(),
753 [&i
, &a
, &b
, &c
](mp_type
& new_value
)
755 const mp_type
x(-(boost::math::constants::euler
<mp_type
>() * i
));
757 new_value
= examples::nr_006::luke_ccoef6_hypergeometric_1f2(a
, b
, c
, x
);
762 total_time
+= STD_CHRONO::duration_cast
<STD_CHRONO::duration
<float> >(my_stopwatch
.elapsed()).count();
764 // Print the values of Hypergeometric1F2.
765 std::for_each(hypergeometric_1f2_results
.begin(),
766 hypergeometric_1f2_results
.end(),
769 std::cout
<< std::setprecision(DIGIT_COUNT
) << h
<< std::endl
;
772 std::cout
<< "test hypergeometric_2f1." << std::endl
;
774 my_stopwatch
.reset();
776 // Generate a table of values of Hypergeometric2F1.
777 // Compare with the Mathematica command:
778 // Table[N[HypergeometricPFQ[{3/7, 2/3}, {1/4}, -(i * EulerGamma)], 100], {i, 1, 20, 1}]
779 std::for_each(hypergeometric_2f1_results
.begin(),
780 hypergeometric_2f1_results
.end(),
781 [&i
, &a
, &b
, &c
](mp_type
& new_value
)
783 const mp_type
x(-(boost::math::constants::euler
<mp_type
>() * i
));
785 new_value
= examples::nr_006::luke_ccoef2_hypergeometric_2f1(a
, b
, c
, x
);
790 total_time
+= STD_CHRONO::duration_cast
<STD_CHRONO::duration
<float> >(my_stopwatch
.elapsed()).count();
792 // Print the values of Hypergeometric2F1.
793 std::for_each(hypergeometric_2f1_results
.begin(),
794 hypergeometric_2f1_results
.end(),
797 std::cout
<< std::setprecision(DIGIT_COUNT
) << h
<< std::endl
;
800 std::cout
<< "Total execution time = " << std::setprecision(3) << total_time
<< "s" << std::endl
;