libs/multiprecision/example/hypergeometric_luke_algorithms.cpp

   1
   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)
   8 //
   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
  12 //
  13
  14 #include <algorithm>
  15 #include <cstdint>
  16 #include <deque>
  17 #include <functional>
  18 #include <iostream>
  19 #include <limits>
  20 #include <numeric>
  21 #include <vector>
  22 #include <boost/math/constants/constants.hpp>
  23 #include <boost/noncopyable.hpp>
  24
  25 //#define USE_CPP_BIN_FLOAT
  26 #define USE_CPP_DEC_FLOAT
  27 //#define USE_MPFR
  28
  29 #if !defined(DIGIT_COUNT)
  30 #define DIGIT_COUNT 100
  31 #endif
  32
  33 #if !defined(BOOST_NO_CXX11_HDR_CHRONO)
  34   #include <chrono>
  35   #define STD_CHRONO std::chrono
  36 #else
  37   #include <boost/chrono.hpp>
  38   #define STD_CHRONO boost::chrono
  39 #endif
  40
  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;
  50 #else
  51   #error no multiprecision floating type is defined
  52 #endif
  53
  54 template <class clock_type>
  55 struct stopwatch
  56 {
  57 public:
  58   typedef typename clock_type::duration duration_type;
  59
  60   stopwatch() : m_start(clock_type::now()) { }
  61
  62   stopwatch(const stopwatch& other) : m_start(other.m_start) { }
  63
  64   stopwatch& operator=(const stopwatch& other)
  65   {
  66     m_start = other.m_start;
  67     return *this;
  68   }
  69
  70   ~stopwatch() { }
  71
  72   duration_type elapsed() const
  73   {
  74     return (clock_type::now() - m_start);
  75   }
  76
  77   void reset()
  78   {
  79     m_start = clock_type::now();
  80   }
  81
  82 private:
  83   typename clock_type::time_point m_start;
  84 };
  85
  86 namespace my_math
  87 {
  88   template<class T> T chebyshev_t(const std::int32_t n, const T& x);
  89
  90   template<class T> T chebyshev_t(const std::uint32_t n, const T& x, std::vector<T>* vp);
  91
  92   template<class T> bool isneg(const T& x) { return (x < T(0)); }
  93
  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; }
  97 }
  98
  99 namespace orthogonal_polynomial_series
 100 {
 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))
 102   {
 103     // Compute the value of an orthogonal chebyshev polinomial.
 104     // Use stable upward recursion.
 105
 106     if(vp != nullptr)
 107     {
 108       vp->clear();
 109       vp->reserve(static_cast<std::size_t>(n + 1u));
 110     }
 111
 112     T y0 = my_math::one<T>();
 113
 114     if(vp != nullptr) { vp->push_back(y0); }
 115
 116     if(n == static_cast<std::uint32_t>(0u))
 117     {
 118       return y0;
 119     }
 120
 121     T y1 = x;
 122
 123     if(vp != nullptr) { vp->push_back(y1); }
 124
 125     if(n == static_cast<std::uint32_t>(1u))
 126     {
 127       return y1;
 128     }
 129
 130     T a = my_math::two <T>();
 131     T b = my_math::zero<T>();
 132     T c = my_math::one <T>();
 133
 134     T yk;
 135
 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)
 139     {
 140       yk = (((a * x) + b) * y1) - (c * y0);
 141
 142       y0 = y1;
 143       y1 = yk;
 144
 145       if(vp != nullptr) { vp->push_back(yk); }
 146     }
 147
 148     return yk;
 149   }
 150 }
 151
 152 template<class T> T my_math::chebyshev_t(const std::int32_t n, const T& x)
 153 {
 154   if(my_math::isneg(x))
 155   {
 156     const bool b_negate = ((n % static_cast<std::int32_t>(2)) != static_cast<std::int32_t>(0));
 157
 158     const T y = chebyshev_t(n, -x);
 159
 160     return (!b_negate ? y : -y);
 161   }
 162
 163   if(n < static_cast<std::int32_t>(0))
 164   {
 165     const std::int32_t nn = static_cast<std::int32_t>(-n);
 166
 167     return chebyshev_t(nn, x);
 168   }
 169   else
 170   {
 171     return orthogonal_polynomial_series::orthogonal_polynomial_template(x, static_cast<std::uint32_t>(n));
 172   }
 173 }
 174
 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); }
 176
 177 namespace util
 178 {
 179   template <class T> float digit_scale()
 180   {
 181     const int d = ((std::max)(std::numeric_limits<T>::digits10, 15));
 182     return static_cast<float>(d) / 300.0F;
 183   }
 184 }
 185
 186 namespace examples
 187 {
 188   namespace nr_006
 189   {
 190     template<typename T> class hypergeometric_pfq_base : private boost::noncopyable
 191     {
 192     public:
 193       virtual ~hypergeometric_pfq_base() { }
 194
 195       virtual void ccoef() const = 0;
 196
 197       virtual T series() const
 198       {
 199         using my_math::chebyshev_t;
 200
 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);
 205
 206         chebyshev_t(static_cast<std::uint32_t>(C.size()),
 207                     z_shifted,
 208                     &chebyshev_t_shifted_values);
 209
 210         // Luke: C     ---------- COMPUTE SCALE FACTOR                       ----------
 211         // Luke: C
 212         // Luke: C     ---------- SCALE THE COEFFICIENTS                     ----------
 213         // Luke: C
 214
 215         // The coefficient scaling is preformed after the Chebyshev summation,
 216         // and it is carried out with a single division operation.
 217         bool b_neg = false;
 218
 219         const T scale = std::accumulate(C.begin(),
 220                                         C.end(),
 221                                         T(0),
 222                                         [&b_neg](T scale_sum, const T& ck) -> T
 223                                         {
 224                                           ((!b_neg) ? (scale_sum += ck) : (scale_sum -= ck));
 225                                           b_neg = (!b_neg);
 226                                           return scale_sum;
 227                                         });
 228
 229         // Compute the result of the series expansion using unscaled coefficients.
 230         const T sum = std::inner_product(C.begin(),
 231                                          C.end(),
 232                                          chebyshev_t_shifted_values.begin(),
 233                                          T(0));
 234
 235         // Return the properly scaled result.
 236         return sum / scale;
 237       }
 238
 239     protected:
 240       const   T             Z;
 241       const   T             W;
 242       mutable std::deque<T> C;
 243
 244       hypergeometric_pfq_base(const T& z,
 245                               const T& w) : Z(z),
 246                                             W(w),
 247                                             C(0u) { }
 248
 249       virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale<T>() * 500.0F); }
 250     };
 251
 252     template<typename T> class ccoef4_hypergeometric_0f1 : public hypergeometric_pfq_base<T>
 253     {
 254     public:
 255       ccoef4_hypergeometric_0f1(const T& c,
 256                                 const T& z,
 257                                 const T& w) : hypergeometric_pfq_base<T>(z, w),
 258                                               CP(c) { }
 259
 260       virtual ~ccoef4_hypergeometric_0f1() { }
 261
 262       virtual void ccoef() const
 263       {
 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));
 267
 268         // Luke: C     ---------- START COMPUTING COEFFICIENTS USING         ----------
 269         // Luke: C     ---------- BACKWARD RECURRENCE SCHEME                 ----------
 270         // Luke: C
 271         T A3(0);
 272         T A2(0);
 273         T A1(boost::math::tools::root_epsilon<T>());
 274
 275         hypergeometric_pfq_base<T>::C.resize(1u, A1);
 276
 277         std::int32_t X1 = N2;
 278
 279         T C1 = T(1) - CP;
 280
 281         const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
 282
 283         for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
 284         {
 285           const T DIVFAC = T(1) / X1;
 286
 287           --X1;
 288
 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)));
 294
 295           hypergeometric_pfq_base<T>::C.push_front(term);
 296
 297           A3 = A2;
 298           A2 = A1;
 299           A1 = hypergeometric_pfq_base<T>::C.front();
 300         }
 301
 302         hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
 303       }
 304
 305     private:
 306       const T CP;
 307     };
 308
 309     template<typename T> class ccoef1_hypergeometric_1f0 : public hypergeometric_pfq_base<T>
 310     {
 311     public:
 312       ccoef1_hypergeometric_1f0(const T& a,
 313                                 const T& z,
 314                                 const T& w) : hypergeometric_pfq_base<T>(z, w),
 315                                               AP(a) { }
 316
 317       virtual ~ccoef1_hypergeometric_1f0() { }
 318
 319       virtual void ccoef() const
 320       {
 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));
 324
 325         // Luke: C     ---------- START COMPUTING COEFFICIENTS USING         ----------
 326         // Luke: C     ---------- BACKWARD RECURRENCE SCHEME                 ----------
 327         // Luke: C
 328         T A2(0);
 329         T A1(boost::math::tools::root_epsilon<T>());
 330
 331         hypergeometric_pfq_base<T>::C.resize(1u, A1);
 332
 333         std::int32_t X1 = N2;
 334
 335         T V1 = T(1) - AP;
 336
 337         // Here, we have corrected what appears to be an error in Luke's code.
 338
 339         // Luke's original code listing has:
 340         //  AFAC = 2 + FOUR/W
 341         // But it appears as though the correct form is:
 342         //  AFAC = 2 - FOUR/W.
 343
 344         const T AFAC = 2 - (T(4) / hypergeometric_pfq_base<T>::W);
 345
 346         for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
 347         {
 348           --X1;
 349
 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);
 353
 354           hypergeometric_pfq_base<T>::C.push_front(term);
 355
 356           A2 = A1;
 357           A1 = hypergeometric_pfq_base<T>::C.front();
 358         }
 359
 360         hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
 361       }
 362
 363     private:
 364       const T AP;
 365
 366       virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale<T>() * 1600.0F); }
 367     };
 368
 369     template<typename T> class ccoef3_hypergeometric_1f1 : public hypergeometric_pfq_base<T>
 370     {
 371     public:
 372       ccoef3_hypergeometric_1f1(const T& a,
 373                                 const T& c,
 374                                 const T& z,
 375                                 const T& w) : hypergeometric_pfq_base<T>(z, w),
 376                                               AP(a),
 377                                               CP(c) { }
 378
 379       virtual ~ccoef3_hypergeometric_1f1() { }
 380
 381       virtual void ccoef() const
 382       {
 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));
 386
 387         // Luke: C     ---------- START COMPUTING COEFFICIENTS USING         ----------
 388         // Luke: C     ---------- BACKWARD RECURRENCE SCHEME                 ----------
 389         // Luke: C
 390         T A3(0);
 391         T A2(0);
 392         T A1(boost::math::tools::root_epsilon<T>());
 393
 394         hypergeometric_pfq_base<T>::C.resize(1u, A1);
 395
 396         std::int32_t X  = N1;
 397         std::int32_t X1 = N2;
 398
 399         T XA  =  X + AP;
 400         T X3A = (X + 3) - AP;
 401
 402         const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
 403
 404         for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
 405         {
 406           --X;
 407           --X1;
 408           --XA;
 409           --X3A;
 410
 411           const T X3A_over_X2 = X3A / static_cast<std::int32_t>(X + 2);
 412
 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;
 418
 419           const T term = (((PART1 + PART2) + PART3) * X1) / XA;
 420
 421           hypergeometric_pfq_base<T>::C.push_front(term);
 422
 423           A3 = A2;
 424           A2 = A1;
 425           A1 = hypergeometric_pfq_base<T>::C.front();
 426         }
 427
 428         hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
 429       }
 430
 431     private:
 432       const T AP;
 433       const T CP;
 434     };
 435
 436     template<typename T> class ccoef6_hypergeometric_1f2 : public hypergeometric_pfq_base<T>
 437     {
 438     public:
 439       ccoef6_hypergeometric_1f2(const T& a,
 440                                 const T& b,
 441                                 const T& c,
 442                                 const T& z,
 443                                 const T& w) : hypergeometric_pfq_base<T>(z, w),
 444                                               AP(a),
 445                                               BP(b),
 446                                               CP(c) { }
 447
 448       virtual ~ccoef6_hypergeometric_1f2() { }
 449
 450       virtual void ccoef() const
 451       {
 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));
 454
 455         // Luke: C     ---------- START COMPUTING COEFFICIENTS USING         ----------
 456         // Luke: C     ---------- BACKWARD RECURRENCE SCHEME                 ----------
 457         // Luke: C
 458         T A4(0);
 459         T A3(0);
 460         T A2(0);
 461         T A1(boost::math::tools::root_epsilon<T>());
 462
 463         hypergeometric_pfq_base<T>::C.resize(1u, A1);
 464
 465         std::int32_t X  = N1;
 466         T            PP = X + AP;
 467
 468         const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
 469
 470         for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
 471         {
 472           --X;
 473           --PP;
 474
 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);
 479
 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);
 482
 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;
 489
 490           const T term = (((PART1 - PART2) + (PART3 - PART4)) * X_PLUS_1) / PP;
 491
 492           hypergeometric_pfq_base<T>::C.push_front(term);
 493
 494           A4 = A3;
 495           A3 = A2;
 496           A2 = A1;
 497           A1 = hypergeometric_pfq_base<T>::C.front();
 498         }
 499
 500         hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
 501       }
 502
 503     private:
 504       const T AP;
 505       const T BP;
 506       const T CP;
 507     };
 508
 509     template<typename T> class ccoef2_hypergeometric_2f1 : public hypergeometric_pfq_base<T>
 510     {
 511     public:
 512       ccoef2_hypergeometric_2f1(const T& a,
 513                                 const T& b,
 514                                 const T& c,
 515                                 const T& z,
 516                                 const T& w) : hypergeometric_pfq_base<T>(z, w),
 517                                               AP(a),
 518                                               BP(b),
 519                                               CP(c) { }
 520
 521       virtual ~ccoef2_hypergeometric_2f1() { }
 522
 523       virtual void ccoef() const
 524       {
 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));
 528
 529         // Luke: C     ---------- START COMPUTING COEFFICIENTS USING         ----------
 530         // Luke: C     ---------- BACKWARD RECURRENCE SCHEME                 ----------
 531         // Luke: C
 532         T A3(0);
 533         T A2(0);
 534         T A1(boost::math::tools::root_epsilon<T>());
 535
 536         hypergeometric_pfq_base<T>::C.resize(1u, A1);
 537
 538         std::int32_t X  = N1;
 539         std::int32_t X1 = N2;
 540         std::int32_t X3 = static_cast<std::int32_t>((X * 2) + 3);
 541
 542         T X3A = (X + 3) - AP;
 543         T X3B = (X + 3) - BP;
 544
 545         const T Z1 = T(4) / hypergeometric_pfq_base<T>::W;
 546
 547         for(std::int32_t k = static_cast<std::int32_t>(0); k < N1; ++k)
 548         {
 549           --X;
 550           --X1;
 551           --X3A;
 552           --X3B;
 553           X3 -= 2;
 554
 555           const std::int32_t X_PLUS_2 = static_cast<std::int32_t>(X + 2);
 556
 557           const T XAB = T(1) / ((X + AP) * (X + BP));
 558
 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);
 564
 565           const T term = (PART1 + PART2) - PART3;
 566
 567           hypergeometric_pfq_base<T>::C.push_front(term);
 568
 569           A3 = A2;
 570           A2 = A1;
 571           A1 = hypergeometric_pfq_base<T>::C.front();
 572         }
 573
 574         hypergeometric_pfq_base<T>::C.front() /= static_cast<std::int32_t>(2);
 575       }
 576
 577     private:
 578       const T AP;
 579       const T BP;
 580       const T CP;
 581
 582       virtual std::int32_t N() const { return static_cast<std::int32_t>(util::digit_scale<T>() * 1600.0F); }
 583     };
 584
 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);
 590   }
 591 }
 592
 593 template<class T>
 594 T examples::nr_006::luke_ccoef4_hypergeometric_0f1(const T& a, const T& x)
 595 {
 596   const ccoef4_hypergeometric_0f1<T> hypergeometric_0f1_object(a, x, T(-20));
 597
 598   hypergeometric_0f1_object.ccoef();
 599
 600   return hypergeometric_0f1_object.series();
 601 }
 602
 603 template<class T>
 604 T examples::nr_006::luke_ccoef1_hypergeometric_1f0(const T& a, const T& x)
 605 {
 606   const ccoef1_hypergeometric_1f0<T> hypergeometric_1f0_object(a, x, T(-20));
 607
 608   hypergeometric_1f0_object.ccoef();
 609
 610   return hypergeometric_1f0_object.series();
 611 }
 612
 613 template<class T>
 614 T examples::nr_006::luke_ccoef3_hypergeometric_1f1(const T& a, const T& b, const T& x)
 615 {
 616   const ccoef3_hypergeometric_1f1<T> hypergeometric_1f1_object(a, b, x, T(-20));
 617
 618   hypergeometric_1f1_object.ccoef();
 619
 620   return hypergeometric_1f1_object.series();
 621 }
 622
 623 template<class T>
 624 T examples::nr_006::luke_ccoef6_hypergeometric_1f2(const T& a, const T& b, const T& c, const T& x)
 625 {
 626   const ccoef6_hypergeometric_1f2<T> hypergeometric_1f2_object(a, b, c, x, T(-20));
 627
 628   hypergeometric_1f2_object.ccoef();
 629
 630   return hypergeometric_1f2_object.series();
 631 }
 632
 633 template<class T>
 634 T examples::nr_006::luke_ccoef2_hypergeometric_2f1(const T& a, const T& b, const T& c, const T& x)
 635 {
 636   const ccoef2_hypergeometric_2f1<T> hypergeometric_2f1_object(a, b, c, x, T(-20));
 637
 638   hypergeometric_2f1_object.ccoef();
 639
 640   return hypergeometric_2f1_object.series();
 641 }
 642
 643 int main()
 644 {
 645   stopwatch<STD_CHRONO::high_resolution_clock> my_stopwatch;
 646   float total_time = 0.0F;
 647
 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);
 653
 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);
 657
 658   std::int_least16_t i;
 659
 660   std::cout << "test hypergeometric_0f1." << std::endl;
 661   i = 1U;
 662   my_stopwatch.reset();
 663
 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)
 670                 {
 671                   const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
 672
 673                   new_value = examples::nr_006::luke_ccoef4_hypergeometric_0f1(a, x);
 674
 675                   ++i;
 676                 });
 677
 678   total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count();
 679
 680   // Print the values of Hypergeometric0F1.
 681   std::for_each(hypergeometric_0f1_results.begin(),
 682                 hypergeometric_0f1_results.end(),
 683                 [](const mp_type& h)
 684                 {
 685                   std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
 686                 });
 687
 688   std::cout << "test hypergeometric_1f0." << std::endl;
 689   i = 1U;
 690   my_stopwatch.reset();
 691
 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)
 698                 {
 699                   const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
 700
 701                   new_value = examples::nr_006::luke_ccoef1_hypergeometric_1f0(a, x);
 702
 703                   ++i;
 704                 });
 705
 706   total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count();
 707
 708   // Print the values of Hypergeometric1F0.
 709   std::for_each(hypergeometric_1f0_results.begin(),
 710                 hypergeometric_1f0_results.end(),
 711                 [](const mp_type& h)
 712                 {
 713                   std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
 714                 });
 715
 716   std::cout << "test hypergeometric_1f1." << std::endl;
 717   i = 1U;
 718   my_stopwatch.reset();
 719
 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)
 726                 {
 727                   const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
 728
 729                   new_value = examples::nr_006::luke_ccoef3_hypergeometric_1f1(a, b, x);
 730
 731                   ++i;
 732                 });
 733
 734   total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count();
 735
 736   // Print the values of Hypergeometric1F1.
 737   std::for_each(hypergeometric_1f1_results.begin(),
 738                 hypergeometric_1f1_results.end(),
 739                 [](const mp_type& h)
 740                 {
 741                   std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
 742                 });
 743
 744   std::cout << "test hypergeometric_1f2." << std::endl;
 745   i = 1U;
 746   my_stopwatch.reset();
 747
 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)
 754                 {
 755                   const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
 756
 757                   new_value = examples::nr_006::luke_ccoef6_hypergeometric_1f2(a, b, c, x);
 758
 759                   ++i;
 760                 });
 761
 762   total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count();
 763
 764   // Print the values of Hypergeometric1F2.
 765   std::for_each(hypergeometric_1f2_results.begin(),
 766                 hypergeometric_1f2_results.end(),
 767                 [](const mp_type& h)
 768                 {
 769                   std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
 770                 });
 771
 772   std::cout << "test hypergeometric_2f1." << std::endl;
 773   i = 1U;
 774   my_stopwatch.reset();
 775
 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)
 782                 {
 783                   const mp_type x(-(boost::math::constants::euler<mp_type>() * i));
 784
 785                   new_value = examples::nr_006::luke_ccoef2_hypergeometric_2f1(a, b, c, x);
 786
 787                   ++i;
 788                 });
 789
 790   total_time += STD_CHRONO::duration_cast<STD_CHRONO::duration<float> >(my_stopwatch.elapsed()).count();
 791
 792   // Print the values of Hypergeometric2F1.
 793   std::for_each(hypergeometric_2f1_results.begin(),
 794                 hypergeometric_2f1_results.end(),
 795                 [](const mp_type& h)
 796                 {
 797                   std::cout << std::setprecision(DIGIT_COUNT) << h << std::endl;
 798                 });
 799
 800   std::cout << "Total execution time = " << std::setprecision(3) << total_time << "s" << std::endl;
 801 }