boost/multiprecision
2
0
Fork 0
mirror of https://github.com/boostorg/multiprecision.git synced 2026-02-01 08:42:09 +00:00
Code Issues Projects Releases Wiki Activity
Files
73a6ec1da430d4e80a16ce804ebc6a49dc54596f
multiprecision/example/hypergeometric_luke_algorithms.cpp

914 lines
30 KiB
C++
Raw Blame History

///////////////////////////////////////////////////////////////////////////////
// Copyright Christopher Kormanyos 2013.
// Copyright John Maddock 2013.
// Distributed under the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
//
// This work is based on an earlier work:
// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
//
#include <algorithm>
#include <array>
#include <cstdint>
#include <deque>
#include <functional>
#include <iostream>
#include <limits>
#include <numeric>
#include <vector>
#include <boost/math/constants/constants.hpp>
#include <boost/noncopyable.hpp>
//#define USE_CPP_BIN_FLOAT
#define USE_CPP_DEC_FLOAT
//#define USE_MPFR
#if !defined(DIGIT_COUNT)
#define DIGIT_COUNT 100
#endif
#define HAS_STD_CHRONO
#if defined(HAS_STD_CHRONO)
#include <chrono>
#define STD_CHRONO std::chrono
#else
#include <boost/chrono.hpp>
#define STD_CHRONO boost::chrono
#endif
#if defined(USE_CPP_BIN_FLOAT)
#include <boost/multiprecision/cpp_bin_float.hpp>
typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<DIGIT_COUNT + 10> > mp_type;
#elif defined(USE_CPP_DEC_FLOAT)
#include <boost/multiprecision/cpp_dec_float.hpp>
typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<DIGIT_COUNT + 10> > mp_type;
#elif defined(USE_MPFR)
#include <boost/multiprecision/mpfr.hpp>
typedef boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<DIGIT_COUNT + 10> > mp_type;
#else
#error no multiprecision floating type is defined
#endif
template <class clock_type>
struct stopwatch
{
public:
typedef typename clock_type::duration duration_type;
stopwatch() : m_start(clock_type::now()) { }
duration_type elapsed() const
{
return clock_type::now() - m_start;
}
void reset()
{
m_start = clock_type::now();
}
private:
typename clock_type::time_point m_start;
};
namespace my_math
{
mp_type chebyshev_t(const std::int32_t n, const mp_type& x);
mp_type chebyshev_u(const std::int32_t n, const mp_type& x);
mp_type hermite (const std::int32_t n, const mp_type& x);
mp_type laguerre (const std::int32_t n, const mp_type& x);
mp_type legendre_p (const std::int32_t n, const mp_type& x);
mp_type legendre_q (const std::int32_t n, const mp_type& x);
mp_type chebyshev_t(const std::uint32_t n, const mp_type& x, std::vector<mp_type>* vp);
mp_type chebyshev_u(const std::uint32_t n, const mp_type& x, std::vector<mp_type>* vp);
mp_type hermite (const std::uint32_t n, const mp_type& x, std::vector<mp_type>* vp);
mp_type laguerre (const std::uint32_t n, const mp_type& x, std::vector<mp_type>* vp);
bool isneg(const mp_type& x) { return (x < mp_type(0)); }
const mp_type& zero() { static const mp_type value_zero(0); return value_zero; }
const mp_type& one () { static const mp_type value_one (1); return value_one; }
const mp_type& two () { static const mp_type value_two (2); return value_two; }
}
namespace orthogonal_polynomial_series
{
typedef enum enum_polynomial_type
{
chebyshev_t_type = 1,
chebyshev_u_type = 2,
laguerre_l_type = 3,
hermite_h_type = 4
}
polynomial_type;
template<typename T> static inline T orthogonal_polynomial_template(const T& x, const std::uint32_t n, const polynomial_type type, std::vector<T>* const vp = static_cast<std::vector<T>*>(0u))
{
// Compute the value of an orthogonal polynomial of one of the following types:
// Chebyshev 1st, Chebyshev 2nd, Laguerre, or Hermite
if(vp != nullptr)
{
vp->clear();
vp->reserve(static_cast<std::size_t>(n + 1u));
}
T y0 = my_math::one();
if(vp != nullptr)
{
vp->push_back(y0);
}
if(n == static_cast<std::uint32_t>(0u))
{
return y0;
}
T y1;
if(type == chebyshev_t_type)
{
y1 = x;
}
else if(type == laguerre_l_type)
{
y1 = my_math::one() - x;
}
else
{
y1 = x * static_cast<std::uint32_t>(2u);
}
if(vp != nullptr)
{
vp->push_back(y1);
}
if(n == static_cast<std::uint32_t>(1u))
{
return y1;
}
T a = my_math::two();
T b = my_math::zero();
T c = my_math::one();
T yk;
// Calculate higher orders using the recurrence relation.
// The direction of stability is upward recursion.
for(std::int32_t k = static_cast<std::int32_t>(2); k <= static_cast<std::int32_t>(n); ++k)
{
if(type == laguerre_l_type)
{
a = -my_math::one() / k;
b = my_math::two() + a;
c = my_math::one() + a;
}
else if(type == hermite_h_type)
{
c = my_math::two() * (k - my_math::one());
}
yk = (((a * x) + b) * y1) - (c * y0);
y0 = y1;
y1 = yk;
if(vp != nullptr)
{
vp->push_back(yk);
}
}
return yk;
}
}
mp_type my_math::chebyshev_t(const std::int32_t n, const mp_type& x)
{
if(my_math::isneg(x))
{
const bool b_negate = ((n % static_cast<std::int32_t>(2)) != static_cast<std::int32_t>(0));
const mp_type y = chebyshev_t(n, -x);
return (!b_negate ? y : -y);
}
if(n < static_cast<std::int32_t>(0))
{
const std::int32_t nn = static_cast<std::int32_t>(-n);
return chebyshev_t(nn, x);
}
else
{
return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::uint32_t>(n), orthogonal_polynomial_series::chebyshev_t_type);
}
}
mp_type my_math::chebyshev_u(const std::int32_t n, const mp_type& x)
{
if(my_math::isneg(x))
{
const bool b_negate = ((n % static_cast<std::int32_t>(2)) != static_cast<std::int32_t>(0));
const mp_type y = chebyshev_u(n, -x);
return ((!b_negate) ? y : -y);
}
if(n < static_cast<std::int32_t>(0))
{
if(n == static_cast<std::int32_t>(-2))
{
return my_math::one();
}
else if(n == static_cast<std::int32_t>(-1))
{
return my_math::zero();
}
else
{
const std::int32_t n_minus_two = static_cast<std::int32_t>(static_cast<std::int32_t>(-n) - static_cast<std::int32_t>(2));
return -chebyshev_u(n_minus_two, x);
}
}
else
{
return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::uint32_t>(n), orthogonal_polynomial_series::chebyshev_u_type);
}
}
mp_type my_math::hermite(const std::int32_t n, const mp_type& x)
{
if(n < static_cast<std::int32_t>(0))
{
// Negative order is not supported.
return my_math::zero();
}
else
{
return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::uint32_t>(n), orthogonal_polynomial_series::hermite_h_type);
}
}
mp_type my_math::laguerre(const std::int32_t n, const mp_type& x)
{
if(n < static_cast<std::int32_t>(0))
{
// Negative order is not supported.
return my_math::zero();
}
else if(n == static_cast<std::int32_t>(0))
{
return my_math::one();
}
else if(n == static_cast<std::int32_t>(1))
{
return my_math::one() - x;
}
if(my_math::isneg(x))
{
// Negative argument is not supported.
return my_math::zero();
}
else
{
return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::uint32_t>(n), orthogonal_polynomial_series::laguerre_l_type);
}
}
mp_type my_math::chebyshev_t(const std::uint32_t n, const mp_type& x, std::vector<mp_type>* const vp) { return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::int32_t>(n), orthogonal_polynomial_series::chebyshev_t_type, vp); }
mp_type my_math::chebyshev_u(const std::uint32_t n, const mp_type& x, std::vector<mp_type>* const vp) { return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::int32_t>(n), orthogonal_polynomial_series::chebyshev_u_type, vp); }
mp_type my_math::hermite (const std::uint32_t n, const mp_type& x, std::vector<mp_type>* const vp) { return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::int32_t>(n), orthogonal_polynomial_series::hermite_h_type, vp); }
mp_type my_math::laguerre (const std::uint32_t n, const mp_type& x, std::vector<mp_type>* const vp) { return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::int32_t>(n), orthogonal_polynomial_series::laguerre_l_type, vp); }
namespace util
{
double digit_scale()
{
return static_cast<double>((std::max)(std::numeric_limits<mp_type>::digits10, 15)) / 300.0;
}
}
namespace examples
{
namespace nr_006
{
template<typename T> class hypergeometric_pfq_base : private boost::noncopyable
{
public:
virtual ~hypergeometric_pfq_base() { }
virtual void ccoef() const = 0;
virtual T series() const
{
using my_math::chebyshev_t;
// Compute the Chebyshev coefficients.
// Get the values of the shifted Chebyshev polynomials.
std::vector<T> chebyshev_t_shifted_values;
const T z_shifted = ((Z / W) * static_cast<std::int32_t>(2)) - static_cast<std::int32_t>(1);
chebyshev_t(static_cast<std::uint32_t>(C.size()),
z_shifted,
&chebyshev_t_shifted_values);
// Luke: C ---------- COMPUTE SCALE FACTOR ----------
// Luke: C
// Luke: C ---------- SCALE THE COEFFICIENTS ----------
// Luke: C
// The coefficient scaling is preformed after the Chebyshev summation,
// and it is carried out with a single division operation.
bool b_neg = false;
const T scale = std::accumulate(C.begin(),
C.end(),
T(0),
[&b_neg](T& scale_sum, const T& ck) -> T
{
((!b_neg) ? (scale_sum += ck) : (scale_sum -= ck));
b_neg = (!b_neg);
return scale_sum;
});
// Compute the result of the series expansion using unscaled coefficients.
const T sum = std::inner_product(C.begin(),
C.end(),
chebyshev_t_shifted_values.begin(),
T(0));
// Return the properly scaled result.
return sum / scale;
}
protected:
const T Z;
const mp_type W;
mutable std::deque<T> C;
hypergeometric_pfq_base(const T& z,
const mp_type& w) : Z(z),
W(w),
C(0u) { }
virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale() * 500.0); }
};
template<typename T> class ccoef4_hypergeometric_0f1 : public hypergeometric_pfq_base<T>
{
public:
ccoef4_hypergeometric_0f1(const T& c,
const T& z,
const mp_type& w) : hypergeometric_pfq_base<T>(z, w),
CP(c) { }
virtual ~ccoef4_hypergeometric_0f1() { }
virtual void ccoef() const
{
// See Luke 1977 page 80.
const std::int32_t N1 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
const std::int32_t N2 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2));
// Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
// Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
// Luke: C
T A3(0);
T A2(0);
T A1(1);
hypergeometric_pfq_base<T>::C.resize(1u, A1);
std::int32_t X1 = N2;
T C1 = T(1) - CP;
const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
{
const T DIVFAC = T(1) / X1;
--X1;
// The terms have been slightly re-arranged resulting in lower complexity.
// Parentheses have been added to avoid reliance on operator precedence.
const T term = (A2 - ((A3 * DIVFAC) * X1))
+ ((A2 * X1) * ((1 + (C1 + X1)) * Z1))
+ ((A1 * X1) * ((DIVFAC - (C1 * Z1)) + (X1 * Z1)));
hypergeometric_pfq_base<T>::C.push_front(term);
A3 = A2;
A2 = A1;
A1 = hypergeometric_pfq_base<T>::C.front();
}
hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
}
private:
const T CP;
};
template<typename T> class ccoef1_hypergeometric_1f0 : public hypergeometric_pfq_base<T>
{
public:
ccoef1_hypergeometric_1f0(const T& a,
const T& z,
const mp_type& w) : hypergeometric_pfq_base<T>(z, w),
AP(a) { }
virtual ~ccoef1_hypergeometric_1f0() { }
virtual void ccoef() const
{
// See Luke 1977 page 67.
const std::int32_t N1 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
const std::int32_t N2 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2));
// Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
// Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
// Luke: C
T A2(0);
T A1(1);
hypergeometric_pfq_base<T>::C.resize(1u, A1);
std::int32_t X1 = N2;
T V1 = T(1) - AP;
// Here, we have corrected what appears to be an error in
// Luke's code. Luke has "AFAC = 2 + FOUR/W". But it appears
// as though "AFAC = 2 - FOUR/W" is correct.
const T AFAC = 2 - (T(4) / hypergeometric_pfq_base<T>::W);
for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
{
--X1;
// The terms have been slightly re-arranged resulting in lower complexity.
// Parentheses have been added to avoid reliance on operator precedence.
const T term = -(X1 * AFAC * A1 + (X1 + V1) * A2) / (X1 - V1);
hypergeometric_pfq_base<T>::C.push_front(term);
A2 = A1;
A1 = hypergeometric_pfq_base<T>::C.front();
}
hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
}
private:
const T AP;
virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale() * 1600.0); }
};
template<typename T> class ccoef3_hypergeometric_1f1 : public hypergeometric_pfq_base<T>
{
public:
ccoef3_hypergeometric_1f1(const T& a,
const T& c,
const T& z,
const mp_type& w) : hypergeometric_pfq_base<T>(z, w),
AP(a),
CP(c) { }
virtual ~ccoef3_hypergeometric_1f1() { }
virtual void ccoef() const
{
// See Luke 1977 page 74.
const std::int32_t N1 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
const std::int32_t N2 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2));
// Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
// Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
// Luke: C
T A3(0);
T A2(0);
T A1(1);
hypergeometric_pfq_base<T>::C.resize(1u, A1);
std::int32_t X = N1;
std::int32_t X1 = N2;
T XA = X + AP;
T X3A = (X + 3) - AP;
const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
{
--X;
--X1;
--XA;
--X3A;
const T X3A_over_X2 = X3A / static_cast<std::int32_t>(X + 2);
// The terms have been slightly re-arranged resulting in lower complexity.
// Parentheses have been added to avoid reliance on operator precedence.
const T PART1 = A1 * (((X + CP) * Z1) - X3A_over_X2);
const T PART2 = A2 * (Z1 * ((X + 3) - CP) + (XA / X1));
const T PART3 = A3 * X3A_over_X2;
const T term = (((PART1 + PART2) + PART3) * X1) / XA;
hypergeometric_pfq_base<T>::C.push_front(term);
A3 = A2;
A2 = A1;
A1 = hypergeometric_pfq_base<T>::C.front();
}
hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
}
private:
const T AP;
const T CP;
};
template<typename T> class ccoef6_hypergeometric_1f2 : public hypergeometric_pfq_base<T>
{
public:
ccoef6_hypergeometric_1f2(const T& a,
const T& b,
const T& c,
const T& z,
const mp_type& w) : hypergeometric_pfq_base<T>(z, w),
AP(a),
BP(b),
CP(c) { }
virtual ~ccoef6_hypergeometric_1f2() { }
virtual void ccoef() const
{
// See Luke 1977 page 85.
const std::int32_t N1 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
// Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
// Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
// Luke: C
T A4(0);
T A3(0);
T A2(0);
T A1(1);
hypergeometric_pfq_base<T>::C.resize(1u, A1);
std::int32_t X = N1;
T PP = X + AP;
const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
{
--X;
--PP;
const std::int32_t TWO_X = static_cast<std::int32_t>(X * 2);
const std::int32_t X_PLUS_1 = static_cast<std::int32_t>(X + 1);
const std::int32_t X_PLUS_3 = static_cast<std::int32_t>(X + 3);
const std::int32_t X_PLUS_4 = static_cast<std::int32_t>(X + 4);
const T QQ = T(TWO_X + 3) / static_cast<std::int32_t>(TWO_X + static_cast<std::int32_t>(5));
const T SS = (X + BP) * (X + CP);
// The terms have been slightly re-arranged resulting in lower complexity.
// Parentheses have been added to avoid reliance on operator precedence.
const T PART1 = A1 * (((PP - (QQ * (PP + 1))) * 2) + (SS * Z1));
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))));
const T PART3 = A3 * ((((X_PLUS_3 - AP) - (QQ * (X_PLUS_4 - AP))) * 2) + (((QQ * Z1) * (X_PLUS_4 - BP)) * (X_PLUS_4 - CP)));
const T PART4 = ((A4 * QQ) * (X_PLUS_4 - AP)) / X_PLUS_3;
const T term = (((PART1 - PART2) + (PART3 - PART4)) * X_PLUS_1) / PP;
hypergeometric_pfq_base<T>::C.push_front(term);
A4 = A3;
A3 = A2;
A2 = A1;
A1 = hypergeometric_pfq_base<T>::C.front();
}
hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
}
private:
const T AP;
const T BP;
const T CP;
};
template<typename T> class ccoef2_hypergeometric_2f1 : public hypergeometric_pfq_base<T>
{
public:
ccoef2_hypergeometric_2f1(const T& a,
const T& b,
const T& c,
const T& z,
const mp_type& w) : hypergeometric_pfq_base<T>(z, w),
AP(a),
BP(b),
CP(c) { }
virtual ~ccoef2_hypergeometric_2f1() { }
virtual void ccoef() const
{
// See Luke 1977 page 59.
const std::int32_t N1 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(1));
const std::int32_t N2 = static_cast<std::int32_t>(N() + static_cast<std::int32_t>(2));
// Luke: C ---------- START COMPUTING COEFFICIENTS USING ----------
// Luke: C ---------- BACKWARD RECURRENCE SCHEME ----------
// Luke: C
T A3(0);
T A2(0);
T A1(1);
hypergeometric_pfq_base<T>::C.resize(1u, A1);
std::int32_t X = N1;
std::int32_t X1 = N2;
std::int32_t X3 = static_cast<std::int32_t>((X * 2) + 3);
T X3A = (X + 3) - AP;
T X3B = (X + 3) - BP;
const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
{
--X;
--X1;
--X3A;
--X3B;
X3 -= 2;
const std::int32_t X_PLUS_2 = static_cast<std::int32_t>(X + 2);
const T XAB = T(1) / ((X + AP) * (X + BP));
// The terms have been slightly re-arranged resulting in lower complexity.
// Parentheses have been added to avoid reliance on operator precedence.
const T PART1 = (A1 * X1) * (2 - (((AP + X1) * (BP + X1)) * ((T(X3) / X_PLUS_2) * XAB)) + ((CP + X) * (XAB * Z1)));
const T PART2 = (A2 * XAB) * ((X3A * X3B) - (X3 * ((X3A + X3B) - 1)) + (((3 - CP) + X) * (X1 * Z1)));
const T PART3 = (A3 * X1) * (X3A / X_PLUS_2) * (X3B * XAB);
const T term = (PART1 + PART2) - PART3;
hypergeometric_pfq_base<T>::C.push_front(term);
A3 = A2;
A2 = A1;
A1 = hypergeometric_pfq_base<T>::C.front();
}
hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
}
private:
const T AP;
const T BP;
const T CP;
virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale() * 1600.0); }
};
mp_type luke_ccoef4_hypergeometric_0f1(const mp_type& a, const mp_type& x);
mp_type luke_ccoef1_hypergeometric_1f0(const mp_type& a, const mp_type& x);
mp_type luke_ccoef3_hypergeometric_1f1(const mp_type& a, const mp_type& b, const mp_type& x);
mp_type luke_ccoef6_hypergeometric_1f2(const mp_type& a, const mp_type& b, const mp_type& c, const mp_type& x);
mp_type luke_ccoef2_hypergeometric_2f1(const mp_type& a, const mp_type& b, const mp_type& c, const mp_type& x);
}
}
mp_type examples::nr_006::luke_ccoef4_hypergeometric_0f1(const mp_type& a, const mp_type& x)
{
const ccoef4_hypergeometric_0f1<mp_type> hypergeometric_0f1_object(a, x, mp_type(-20));
hypergeometric_0f1_object.ccoef();
return hypergeometric_0f1_object.series();
}
mp_type examples::nr_006::luke_ccoef1_hypergeometric_1f0(const mp_type& a, const mp_type& x)
{
const ccoef1_hypergeometric_1f0<mp_type> hypergeometric_1f0_object(a, x, mp_type(-20));
hypergeometric_1f0_object.ccoef();
return hypergeometric_1f0_object.series();
}
mp_type examples::nr_006::luke_ccoef3_hypergeometric_1f1(const mp_type& a, const mp_type& b, const mp_type& x)
{
const ccoef3_hypergeometric_1f1<mp_type> hypergeometric_1f1_object(a, b, x, mp_type(-20));
hypergeometric_1f1_object.ccoef();
return hypergeometric_1f1_object.series();
}
mp_type examples::nr_006::luke_ccoef6_hypergeometric_1f2(const mp_type& a, const mp_type& b, const mp_type& c, const mp_type& x)
{
const ccoef6_hypergeometric_1f2<mp_type> hypergeometric_1f2_object(a, b, c, x, mp_type(-20));
hypergeometric_1f2_object.ccoef();
return hypergeometric_1f2_object.series();
}
mp_type examples::nr_006::luke_ccoef2_hypergeometric_2f1(const mp_type& a, const mp_type& b, const mp_type& c, const mp_type& x)
{
const ccoef2_hypergeometric_2f1<mp_type> hypergeometric_2f1_object(a, b, c, x, mp_type(-20));
hypergeometric_2f1_object.ccoef();
return hypergeometric_2f1_object.series();
}
int main()
{
stopwatch<STD_CHRONO::high_resolution_clock> my_stopwatch;
double total_time = 0.0;
std::vector<mp_type> hypergeometric_0f1_results(20U);
std::vector<mp_type> hypergeometric_1f0_results(20U);
std::vector<mp_type> hypergeometric_1f1_results(20U);
std::vector<mp_type> hypergeometric_2f1_results(20U);
std::vector<mp_type> hypergeometric_1f2_results(20U);
const mp_type a(mp_type(3) / 7);
const mp_type b(mp_type(2) / 3);
const mp_type c(mp_type(1) / 4);
std::int_least16_t i;
std::cout << "test hypergeometric_0f1." << std::endl;
i = 1U;
my_stopwatch.reset();
// Generate a table of values of Hypergeometric0F1.
// Compare with the Mathematica command:
// Table[N[HypergeometricPFQ[{}, {3/7}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
std::for_each(hypergeometric_0f1_results.begin(),
hypergeometric_0f1_results.end(),
[&i, &a](mp_type& new_value)
{
const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
new_value = examples::nr_006::luke_ccoef4_hypergeometric_0f1(a, x);
++i;
});
total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<double> >(my_stopwatch.elapsed()).count();
// Print the values of Hypergeometric0F1.
std::for_each(hypergeometric_0f1_results.begin(),
hypergeometric_0f1_results.end(),
[](const mp_type& h)
{
std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
});
std::cout << "test hypergeometric_1f0." << std::endl;
i = 1U;
my_stopwatch.reset();
// Generate a table of values of Hypergeometric1F0.
// Compare with the Mathematica command:
// Table[N[HypergeometricPFQ[{3/7}, {}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
std::for_each(hypergeometric_1f0_results.begin(),
hypergeometric_1f0_results.end(),
[&i, &a](mp_type& new_value)
{
const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
new_value = examples::nr_006::luke_ccoef1_hypergeometric_1f0(a, x);
++i;
});
total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<double> >(my_stopwatch.elapsed()).count();
// Print the values of Hypergeometric1F0.
std::for_each(hypergeometric_1f0_results.begin(),
hypergeometric_1f0_results.end(),
[](const mp_type& h)
{
std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
});
std::cout << "test hypergeometric_1f1." << std::endl;
i = 1U;
my_stopwatch.reset();
// Generate a table of values of Hypergeometric1F1.
// Compare with the Mathematica command:
// Table[N[HypergeometricPFQ[{3/7}, {2/3}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
std::for_each(hypergeometric_1f1_results.begin(),
hypergeometric_1f1_results.end(),
[&i, &a, &b](mp_type& new_value)
{
const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
new_value = examples::nr_006::luke_ccoef3_hypergeometric_1f1(a, b, x);
++i;
});
total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<double> >(my_stopwatch.elapsed()).count();
// Print the values of Hypergeometric1F1.
std::for_each(hypergeometric_1f1_results.begin(),
hypergeometric_1f1_results.end(),
[](const mp_type& h)
{
std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
});
std::cout << "test hypergeometric_1f2." << std::endl;
i = 1U;
my_stopwatch.reset();
// Generate a table of values of Hypergeometric1F2.
// Compare with the Mathematica command:
// Table[N[HypergeometricPFQ[{3/7}, {2/3, 1/4}, -(i*EulerGamma)], 100], {i, 1, 20, 1}]
std::for_each(hypergeometric_1f2_results.begin(),
hypergeometric_1f2_results.end(),
[&i, &a, &b, &c](mp_type& new_value)
{
const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
new_value = examples::nr_006::luke_ccoef6_hypergeometric_1f2(a, b, c, x);
++i;
});
total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<double> >(my_stopwatch.elapsed()).count();
// Print the values of Hypergeometric1F2.
std::for_each(hypergeometric_1f2_results.begin(),
hypergeometric_1f2_results.end(),
[](const mp_type& h)
{
std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
});
std::cout << "test hypergeometric_2f1." << std::endl;
i = 1U;
my_stopwatch.reset();
// Generate a table of values of Hypergeometric2F1.
// Compare with the Mathematica command:
// Table[N[HypergeometricPFQ[{3/7, 2/3}, {1/4}, -(i * EulerGamma)], 100], {i, 1, 20, 1}]
std::for_each(hypergeometric_2f1_results.begin(),
hypergeometric_2f1_results.end(),
[&i, &a, &b, &c](mp_type& new_value)
{
const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
new_value = examples::nr_006::luke_ccoef2_hypergeometric_2f1(a, b, c, x);
++i;
});
total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<double> >(my_stopwatch.elapsed()).count();
// Print the values of Hypergeometric2F1.
std::for_each(hypergeometric_2f1_results.begin(),
hypergeometric_2f1_results.end(),
[](const mp_type& h)
{
std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
});
std::cout << "Total execution time = " << std::setprecision(3) << total_time << "s" << std::endl;
return 0;
}
Reference in New Issue View Git Blame Copy Permalink