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
35 #if defined(USE_CPP_BIN_FLOAT)
36 #include <boost/multiprecision/cpp_bin_float.hpp>
37 typedef boost::multiprecision::number
<boost::multiprecision::cpp_bin_float
<DIGIT_COUNT
+ 10> > mp_type
;
38 #elif defined(USE_CPP_DEC_FLOAT)
39 #include <boost/multiprecision/cpp_dec_float.hpp>
40 typedef boost::multiprecision::number
<boost::multiprecision::cpp_dec_float
<DIGIT_COUNT
+ 10> > mp_type
;
41 #elif defined(USE_MPFR)
42 #include <boost/multiprecision/mpfr.hpp>
43 typedef boost::multiprecision::number
<boost::multiprecision::mpfr_float_backend
<DIGIT_COUNT
+ 10> > mp_type
;
45 #error no multiprecision floating type is defined
48 template <class clock_type
>
52 typedef typename
clock_type::duration duration_type
;
54 stopwatch() : m_start(clock_type::now()) { }
56 stopwatch(const stopwatch
& other
) : m_start(other
.m_start
) { }
58 stopwatch
& operator=(const stopwatch
& other
)
60 m_start
= other
.m_start
;
66 duration_type
elapsed() const
68 return (clock_type::now() - m_start
);
73 m_start
= clock_type::now();
77 typename
clock_type::time_point m_start
;
82 template<class T
> T
chebyshev_t(const std::int32_t n
, const T
& x
);
84 template<class T
> T
chebyshev_t(const std::uint32_t n
, const T
& x
, std::vector
<T
>* vp
);
86 template<class T
> bool isneg(const T
& x
) { return (x
< T(0)); }
88 template<class T
> const T
& zero() { static const T
value_zero(0); return value_zero
; }
89 template<class T
> const T
& one () { static const T
value_one (1); return value_one
; }
90 template<class T
> const T
& two () { static const T
value_two (2); return value_two
; }
93 namespace orthogonal_polynomial_series
95 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))
97 // Compute the value of an orthogonal chebyshev polinomial.
98 // Use stable upward recursion.
103 vp
->reserve(static_cast<std::size_t>(n
+ 1u));
106 T y0
= my_math::one
<T
>();
108 if(vp
!= nullptr) { vp
->push_back(y0
); }
110 if(n
== static_cast<std::uint32_t>(0u))
117 if(vp
!= nullptr) { vp
->push_back(y1
); }
119 if(n
== static_cast<std::uint32_t>(1u))
124 T a
= my_math::two
<T
>();
125 T b
= my_math::zero
<T
>();
126 T c
= my_math::one
<T
>();
130 // Calculate higher orders using the recurrence relation.
131 // The direction of stability is upward recursion.
132 for(std::int32_t k
= static_cast<std::int32_t>(2); k
<= static_cast<std::int32_t>(n
); ++k
)
134 yk
= (((a
* x
) + b
) * y1
) - (c
* y0
);
139 if(vp
!= nullptr) { vp
->push_back(yk
); }
146 template<class T
> T
my_math::chebyshev_t(const std::int32_t n
, const T
& x
)
148 if(my_math::isneg(x
))
150 const bool b_negate
= ((n
% static_cast<std::int32_t>(2)) != static_cast<std::int32_t>(0));
152 const T y
= chebyshev_t(n
, -x
);
154 return (!b_negate
? y
: -y
);
157 if(n
< static_cast<std::int32_t>(0))
159 const std::int32_t nn
= static_cast<std::int32_t>(-n
);
161 return chebyshev_t(nn
, x
);
165 return orthogonal_polynomial_series::orthogonal_polynomial_template(x
, static_cast<std::uint32_t>(n
));
169 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
); }
173 template <class T
> float digit_scale()
175 const int d
= ((std::max
)(std::numeric_limits
<T
>::digits10
, 15));
176 return static_cast<float>(d
) / 300.0F
;
184 template<typename T
> class hypergeometric_pfq_base
: private boost::noncopyable
187 virtual ~hypergeometric_pfq_base() { }
189 virtual void ccoef() const = 0;
191 virtual T
series() const
193 using my_math::chebyshev_t
;
195 // Compute the Chebyshev coefficients.
196 // Get the values of the shifted Chebyshev polynomials.
197 std::vector
<T
> chebyshev_t_shifted_values
;
198 const T z_shifted
= ((Z
/ W
) * static_cast<std::int32_t>(2)) - static_cast<std::int32_t>(1);
200 chebyshev_t(static_cast<std::uint32_t>(C
.size()),
202 &chebyshev_t_shifted_values
);
204 // Luke: C ---------- COMPUTE SCALE FACTOR ----------
206 // Luke: C ---------- SCALE THE COEFFICIENTS ----------
209 // The coefficient scaling is preformed after the Chebyshev summation,
210 // and it is carried out with a single division operation.
213 const T scale
= std::accumulate(C
.begin(),
216 [&b_neg
](T scale_sum
, const T
& ck
) -> T
218 ((!b_neg
) ? (scale_sum
+= ck
) : (scale_sum
-= ck
));
223 // Compute the result of the series expansion using unscaled coefficients.
224 const T sum
= std::inner_product(C
.begin(),
226 chebyshev_t_shifted_values
.begin(),
229 // Return the properly scaled result.
236 mutable std::deque
<T
> C
;
238 hypergeometric_pfq_base(const T
& z
,
243 virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale
<T
>() * 500.0F
); }
246 template<typename T
> class ccoef4_hypergeometric_0f1
: public hypergeometric_pfq_base
<T
>
249 ccoef4_hypergeometric_0f1(const T
& c
,
251 const T
& w
) : hypergeometric_pfq_base
<T
>(z
, w
),
254 virtual ~ccoef4_hypergeometric_0f1() { }
256 virtual void ccoef() const
258 // See Luke 1977 page 80.
259 const std::int32_t N1
= static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1));
260 const std::int32_t N2
= static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(2));
262 // Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
263 // Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
267 T
A1(boost::math::tools::root_epsilon
<T
>());
269 hypergeometric_pfq_base
<T
>::C
.resize(1u, A1
);
271 std::int32_t X1
= N2
;
275 const T Z1
= T(4) / hypergeometric_pfq_base
<T
>::W
;
277 for(std::int32_t k
= static_cast<std::int32_t>(0); k
< N1
; ++k
)
279 const T DIVFAC
= T(1) / X1
;
283 // The terms have been slightly re-arranged resulting in lower complexity.
284 // Parentheses have been added to avoid reliance on operator precedence.
285 const T term
= (A2
- ((A3
* DIVFAC
) * X1
))
286 + ((A2
* X1
) * ((1 + (C1
+ X1
)) * Z1
))
287 + ((A1
* X1
) * ((DIVFAC
- (C1
* Z1
)) + (X1
* Z1
)));
289 hypergeometric_pfq_base
<T
>::C
.push_front(term
);
293 A1
= hypergeometric_pfq_base
<T
>::C
.front();
296 hypergeometric_pfq_base
<T
>::C
.front() /= static_cast<std::int32_t>(2);
303 template<typename T
> class ccoef1_hypergeometric_1f0
: public hypergeometric_pfq_base
<T
>
306 ccoef1_hypergeometric_1f0(const T
& a
,
308 const T
& w
) : hypergeometric_pfq_base
<T
>(z
, w
),
311 virtual ~ccoef1_hypergeometric_1f0() { }
313 virtual void ccoef() const
315 // See Luke 1977 page 67.
316 const std::int32_t N1
= static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
317 const std::int32_t N2
= static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2));
319 // Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
320 // Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
323 T
A1(boost::math::tools::root_epsilon
<T
>());
325 hypergeometric_pfq_base
<T
>::C
.resize(1u, A1
);
327 std::int32_t X1
= N2
;
331 // Here, we have corrected what appears to be an error in Luke's code.
333 // Luke's original code listing has:
335 // But it appears as though the correct form is:
336 // AFAC = 2 - FOUR/W.
338 const T AFAC
= 2 - (T(4) / hypergeometric_pfq_base
<T
>::W
);
340 for(std::int32_t k
= static_cast<std::int32_t>(0); k
< N1
; ++k
)
344 // The terms have been slightly re-arranged resulting in lower complexity.
345 // Parentheses have been added to avoid reliance on operator precedence.
346 const T term
= -(((X1
* AFAC
) * A1
) + ((X1
+ V1
) * A2
)) / (X1
- V1
);
348 hypergeometric_pfq_base
<T
>::C
.push_front(term
);
351 A1
= hypergeometric_pfq_base
<T
>::C
.front();
354 hypergeometric_pfq_base
<T
>::C
.front() /= static_cast<std::int32_t>(2);
360 virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale
<T
>() * 1600.0F
); }
363 template<typename T
> class ccoef3_hypergeometric_1f1
: public hypergeometric_pfq_base
<T
>
366 ccoef3_hypergeometric_1f1(const T
& a
,
369 const T
& w
) : hypergeometric_pfq_base
<T
>(z
, w
),
373 virtual ~ccoef3_hypergeometric_1f1() { }
375 virtual void ccoef() const
377 // See Luke 1977 page 74.
378 const std::int32_t N1
= static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1));
379 const std::int32_t N2
= static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(2));
381 // Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
382 // Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
386 T
A1(boost::math::tools::root_epsilon
<T
>());
388 hypergeometric_pfq_base
<T
>::C
.resize(1u, A1
);
391 std::int32_t X1
= N2
;
394 T X3A
= (X
+ 3) - AP
;
396 const T Z1
= T(4) / hypergeometric_pfq_base
<T
>::W
;
398 for(std::int32_t k
= static_cast<std::int32_t>(0); k
< N1
; ++k
)
405 const T X3A_over_X2
= X3A
/ static_cast<std::int32_t>(X
+ 2);
407 // The terms have been slightly re-arranged resulting in lower complexity.
408 // Parentheses have been added to avoid reliance on operator precedence.
409 const T PART1
= A1
* (((X
+ CP
) * Z1
) - X3A_over_X2
);
410 const T PART2
= A2
* (Z1
* ((X
+ 3) - CP
) + (XA
/ X1
));
411 const T PART3
= A3
* X3A_over_X2
;
413 const T term
= (((PART1
+ PART2
) + PART3
) * X1
) / XA
;
415 hypergeometric_pfq_base
<T
>::C
.push_front(term
);
419 A1
= hypergeometric_pfq_base
<T
>::C
.front();
422 hypergeometric_pfq_base
<T
>::C
.front() /= static_cast<std::int32_t>(2);
430 template<typename T
> class ccoef6_hypergeometric_1f2
: public hypergeometric_pfq_base
<T
>
433 ccoef6_hypergeometric_1f2(const T
& a
,
437 const T
& w
) : hypergeometric_pfq_base
<T
>(z
, w
),
442 virtual ~ccoef6_hypergeometric_1f2() { }
444 virtual void ccoef() const
446 // See Luke 1977 page 85.
447 const std::int32_t N1
= static_cast<std::int32_t>(this->N() + static_cast<std::int32_t>(1));
449 // Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
450 // Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
455 T
A1(boost::math::tools::root_epsilon
<T
>());
457 hypergeometric_pfq_base
<T
>::C
.resize(1u, A1
);
462 const T Z1
= T(4) / hypergeometric_pfq_base
<T
>::W
;
464 for(std::int32_t k
= static_cast<std::int32_t>(0); k
< N1
; ++k
)
469 const std::int32_t TWO_X
= static_cast<std::int32_t>(X
* 2);
470 const std::int32_t X_PLUS_1
= static_cast<std::int32_t>(X
+ 1);
471 const std::int32_t X_PLUS_3
= static_cast<std::int32_t>(X
+ 3);
472 const std::int32_t X_PLUS_4
= static_cast<std::int32_t>(X
+ 4);
474 const T QQ
= T(TWO_X
+ 3) / static_cast<std::int32_t>(TWO_X
+ static_cast<std::int32_t>(5));
475 const T SS
= (X
+ BP
) * (X
+ CP
);
477 // The terms have been slightly re-arranged resulting in lower complexity.
478 // Parentheses have been added to avoid reliance on operator precedence.
479 const T PART1
= A1
* (((PP
- (QQ
* (PP
+ 1))) * 2) + (SS
* Z1
));
480 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
))));
481 const T PART3
= A3
* ((((X_PLUS_3
- AP
) - (QQ
* (X_PLUS_4
- AP
))) * 2) + (((QQ
* Z1
) * (X_PLUS_4
- BP
)) * (X_PLUS_4
- CP
)));
482 const T PART4
= ((A4
* QQ
) * (X_PLUS_4
- AP
)) / X_PLUS_3
;
484 const T term
= (((PART1
- PART2
) + (PART3
- PART4
)) * X_PLUS_1
) / PP
;
486 hypergeometric_pfq_base
<T
>::C
.push_front(term
);
491 A1
= hypergeometric_pfq_base
<T
>::C
.front();
494 hypergeometric_pfq_base
<T
>::C
.front() /= static_cast<std::int32_t>(2);
503 template<typename T
> class ccoef2_hypergeometric_2f1
: public hypergeometric_pfq_base
<T
>
506 ccoef2_hypergeometric_2f1(const T
& a
,
510 const T
& w
) : hypergeometric_pfq_base
<T
>(z
, w
),
515 virtual ~ccoef2_hypergeometric_2f1() { }
517 virtual void ccoef() const
519 // See Luke 1977 page 59.
520 const std::int32_t N1
= static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
521 const std::int32_t N2
= static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2));
523 // Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
524 // Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
528 T
A1(boost::math::tools::root_epsilon
<T
>());
530 hypergeometric_pfq_base
<T
>::C
.resize(1u, A1
);
533 std::int32_t X1
= N2
;
534 std::int32_t X3
= static_cast<std::int32_t>((X
* 2) + 3);
536 T X3A
= (X
+ 3) - AP
;
537 T X3B
= (X
+ 3) - BP
;
539 const T Z1
= T(4) / hypergeometric_pfq_base
<T
>::W
;
541 for(std::int32_t k
= static_cast<std::int32_t>(0); k
< N1
; ++k
)
549 const std::int32_t X_PLUS_2
= static_cast<std::int32_t>(X
+ 2);
551 const T XAB
= T(1) / ((X
+ AP
) * (X
+ BP
));
553 // The terms have been slightly re-arranged resulting in lower complexity.
554 // Parentheses have been added to avoid reliance on operator precedence.
555 const T PART1
= (A1
* X1
) * (2 - (((AP
+ X1
) * (BP
+ X1
)) * ((T(X3
) / X_PLUS_2
) * XAB
)) + ((CP
+ X
) * (XAB
* Z1
)));
556 const T PART2
= (A2
* XAB
) * ((X3A
* X3B
) - (X3
* ((X3A
+ X3B
) - 1)) + (((3 - CP
) + X
) * (X1
* Z1
)));
557 const T PART3
= (A3
* X1
) * (X3A
/ X_PLUS_2
) * (X3B
* XAB
);
559 const T term
= (PART1
+ PART2
) - PART3
;
561 hypergeometric_pfq_base
<T
>::C
.push_front(term
);
565 A1
= hypergeometric_pfq_base
<T
>::C
.front();
568 hypergeometric_pfq_base
<T
>::C
.front() /= static_cast<std::int32_t>(2);
576 virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale
<T
>() * 1600.0F
); }
579 template<class T
> T
luke_ccoef4_hypergeometric_0f1(const T
& a
, const T
& x
);
580 template<class T
> T
luke_ccoef1_hypergeometric_1f0(const T
& a
, const T
& x
);
581 template<class T
> T
luke_ccoef3_hypergeometric_1f1(const T
& a
, const T
& b
, const T
& x
);
582 template<class T
> T
luke_ccoef6_hypergeometric_1f2(const T
& a
, const T
& b
, const T
& c
, const T
& x
);
583 template<class T
> T
luke_ccoef2_hypergeometric_2f1(const T
& a
, const T
& b
, const T
& c
, const T
& x
);
588 T
examples::nr_006::luke_ccoef4_hypergeometric_0f1(const T
& a
, const T
& x
)
590 const ccoef4_hypergeometric_0f1
<T
> hypergeometric_0f1_object(a
, x
, T(-20));
592 hypergeometric_0f1_object
.ccoef();
594 return hypergeometric_0f1_object
.series();
598 T
examples::nr_006::luke_ccoef1_hypergeometric_1f0(const T
& a
, const T
& x
)
600 const ccoef1_hypergeometric_1f0
<T
> hypergeometric_1f0_object(a
, x
, T(-20));
602 hypergeometric_1f0_object
.ccoef();
604 return hypergeometric_1f0_object
.series();
608 T
examples::nr_006::luke_ccoef3_hypergeometric_1f1(const T
& a
, const T
& b
, const T
& x
)
610 const ccoef3_hypergeometric_1f1
<T
> hypergeometric_1f1_object(a
, b
, x
, T(-20));
612 hypergeometric_1f1_object
.ccoef();
614 return hypergeometric_1f1_object
.series();
618 T
examples::nr_006::luke_ccoef6_hypergeometric_1f2(const T
& a
, const T
& b
, const T
& c
, const T
& x
)
620 const ccoef6_hypergeometric_1f2
<T
> hypergeometric_1f2_object(a
, b
, c
, x
, T(-20));
622 hypergeometric_1f2_object
.ccoef();
624 return hypergeometric_1f2_object
.series();
628 T
examples::nr_006::luke_ccoef2_hypergeometric_2f1(const T
& a
, const T
& b
, const T
& c
, const T
& x
)
630 const ccoef2_hypergeometric_2f1
<T
> hypergeometric_2f1_object(a
, b
, c
, x
, T(-20));
632 hypergeometric_2f1_object
.ccoef();
634 return hypergeometric_2f1_object
.series();
639 stopwatch
<std::chrono::high_resolution_clock
> my_stopwatch
;
640 float total_time
= 0.0F
;
642 std::vector
<mp_type
> hypergeometric_0f1_results(20U);
643 std::vector
<mp_type
> hypergeometric_1f0_results(20U);
644 std::vector
<mp_type
> hypergeometric_1f1_results(20U);
645 std::vector
<mp_type
> hypergeometric_2f1_results(20U);
646 std::vector
<mp_type
> hypergeometric_1f2_results(20U);
648 const mp_type
a(mp_type(3) / 7);
649 const mp_type
b(mp_type(2) / 3);
650 const mp_type
c(mp_type(1) / 4);
652 std::int_least16_t i
;
654 std::cout
<< "test hypergeometric_0f1." << std::endl
;
656 my_stopwatch
.reset();
658 // Generate a table of values of Hypergeometric0F1.
659 // Compare with the Mathematica command:
660 // Table[N[HypergeometricPFQ[{}, {3/7}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
661 std::for_each(hypergeometric_0f1_results
.begin(),
662 hypergeometric_0f1_results
.end(),
663 [&i
, &a
](mp_type
& new_value
)
665 const mp_type
x(-(boost::math::constants::euler
<mp_type
>() * i
));
667 new_value
= examples::nr_006::luke_ccoef4_hypergeometric_0f1(a
, x
);
672 total_time
+= std::chrono::duration_cast
<std::chrono::duration
<float> >(my_stopwatch
.elapsed()).count();
674 // Print the values of Hypergeometric0F1.
675 std::for_each(hypergeometric_0f1_results
.begin(),
676 hypergeometric_0f1_results
.end(),
679 std::cout
<< std::setprecision(DIGIT_COUNT
) << h
<< std::endl
;
682 std::cout
<< "test hypergeometric_1f0." << std::endl
;
684 my_stopwatch
.reset();
686 // Generate a table of values of Hypergeometric1F0.
687 // Compare with the Mathematica command:
688 // Table[N[HypergeometricPFQ[{3/7}, {}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
689 std::for_each(hypergeometric_1f0_results
.begin(),
690 hypergeometric_1f0_results
.end(),
691 [&i
, &a
](mp_type
& new_value
)
693 const mp_type
x(-(boost::math::constants::euler
<mp_type
>() * i
));
695 new_value
= examples::nr_006::luke_ccoef1_hypergeometric_1f0(a
, x
);
700 total_time
+= std::chrono::duration_cast
<std::chrono::duration
<float> >(my_stopwatch
.elapsed()).count();
702 // Print the values of Hypergeometric1F0.
703 std::for_each(hypergeometric_1f0_results
.begin(),
704 hypergeometric_1f0_results
.end(),
707 std::cout
<< std::setprecision(DIGIT_COUNT
) << h
<< std::endl
;
710 std::cout
<< "test hypergeometric_1f1." << std::endl
;
712 my_stopwatch
.reset();
714 // Generate a table of values of Hypergeometric1F1.
715 // Compare with the Mathematica command:
716 // Table[N[HypergeometricPFQ[{3/7}, {2/3}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
717 std::for_each(hypergeometric_1f1_results
.begin(),
718 hypergeometric_1f1_results
.end(),
719 [&i
, &a
, &b
](mp_type
& new_value
)
721 const mp_type
x(-(boost::math::constants::euler
<mp_type
>() * i
));
723 new_value
= examples::nr_006::luke_ccoef3_hypergeometric_1f1(a
, b
, x
);
728 total_time
+= std::chrono::duration_cast
<std::chrono::duration
<float> >(my_stopwatch
.elapsed()).count();
730 // Print the values of Hypergeometric1F1.
731 std::for_each(hypergeometric_1f1_results
.begin(),
732 hypergeometric_1f1_results
.end(),
735 std::cout
<< std::setprecision(DIGIT_COUNT
) << h
<< std::endl
;
738 std::cout
<< "test hypergeometric_1f2." << std::endl
;
740 my_stopwatch
.reset();
742 // Generate a table of values of Hypergeometric1F2.
743 // Compare with the Mathematica command:
744 // Table[N[HypergeometricPFQ[{3/7}, {2/3, 1/4}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
745 std::for_each(hypergeometric_1f2_results
.begin(),
746 hypergeometric_1f2_results
.end(),
747 [&i
, &a
, &b
, &c
](mp_type
& new_value
)
749 const mp_type
x(-(boost::math::constants::euler
<mp_type
>() * i
));
751 new_value
= examples::nr_006::luke_ccoef6_hypergeometric_1f2(a
, b
, c
, x
);
756 total_time
+= std::chrono::duration_cast
<std::chrono::duration
<float> >(my_stopwatch
.elapsed()).count();
758 // Print the values of Hypergeometric1F2.
759 std::for_each(hypergeometric_1f2_results
.begin(),
760 hypergeometric_1f2_results
.end(),
763 std::cout
<< std::setprecision(DIGIT_COUNT
) << h
<< std::endl
;
766 std::cout
<< "test hypergeometric_2f1." << std::endl
;
768 my_stopwatch
.reset();
770 // Generate a table of values of Hypergeometric2F1.
771 // Compare with the Mathematica command:
772 // Table[N[HypergeometricPFQ[{3/7, 2/3}, {1/4}, -(i * EulerGamma)], 100], {i, 1, 20, 1}]
773 std::for_each(hypergeometric_2f1_results
.begin(),
774 hypergeometric_2f1_results
.end(),
775 [&i
, &a
, &b
, &c
](mp_type
& new_value
)
777 const mp_type
x(-(boost::math::constants::euler
<mp_type
>() * i
));
779 new_value
= examples::nr_006::luke_ccoef2_hypergeometric_2f1(a
, b
, c
, x
);
784 total_time
+= std::chrono::duration_cast
<std::chrono::duration
<float> >(my_stopwatch
.elapsed()).count();
786 // Print the values of Hypergeometric2F1.
787 std::for_each(hypergeometric_2f1_results
.begin(),
788 hypergeometric_2f1_results
.end(),
791 std::cout
<< std::setprecision(DIGIT_COUNT
) << h
<< std::endl
;
794 std::cout
<< "Total execution time = " << std::setprecision(3) << total_time
<< "s" << std::endl
;