special_functions/detail/hypergeometric_pade.hpp

0001
0002 ///////////////////////////////////////////////////////////////////////////////
0003 //  Copyright 2014 Anton Bikineev
0004 //  Copyright 2014 Christopher Kormanyos
0005 //  Copyright 2014 John Maddock
0006 //  Copyright 2014 Paul Bristow
0007 //  Distributed under the Boost
0008 //  Software License, Version 1.0. (See accompanying file
0009 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0010 //
0011 #ifndef BOOST_MATH_HYPERGEOMETRIC_PADE_HPP
0012 #define BOOST_MATH_HYPERGEOMETRIC_PADE_HPP
0013
0014   namespace boost{ namespace math{ namespace detail{
0015
0016   // Luke: C ---------- SUBROUTINE R1F1P(CP, Z, A, B, N) ----------
0017   // Luke: C ----- PADE APPROXIMATION OF 1F1( 1 ; CP ; -Z ) -------
0018   template <class T, class Policy>
0019   inline T hypergeometric_1F1_pade(const T& cp, const T& zp, const Policy& )
0020   {
0021     BOOST_MATH_STD_USING
0022
0023     static const T one = T(1);
0024
0025     // Luke: C ------------- INITIALIZATION -------------
0026     const T z = -zp;
0027     const T zz = z * z;
0028     T b0 = one;
0029     T a0 = one;
0030     T xi1 = one;
0031     T ct1 = cp + one;
0032     T cp1 = cp - one;
0033
0034     T b1 = one + (z / ct1);
0035     T a1 = b1 - (z / cp);
0036
0037     const unsigned max_iterations = boost::math::policies::get_max_series_iterations<Policy>();
0038
0039     T b2 = T(0), a2 = T(0);
0040     T result = T(0), prev_result;
0041
0042     for (unsigned k = 1; k < max_iterations; ++k)
0043     {
0044       // Luke: C ----- CALCULATION OF THE MULTIPLIERS -----
0045       // Luke: C ----------- FOR THE RECURSION ------------
0046       const T ct2 = ct1 * ct1;
0047       const T g1 = one + ((cp1 / (ct2 + ct1 + ct1)) * z);
0048       const T g2 = ((xi1 / (ct2 - one)) * ((xi1 + cp1) / ct2)) * zz;
0049
0050       // Luke: C ------- THE RECURRENCE RELATIONS ---------
0051       // Luke: C ------------ ARE AS FOLLOWS --------------
0052       b2 = (g1 * b1) + (g2 * b0);
0053       a2 = (g1 * a1) + (g2 * a0);
0054
0055       prev_result = result;
0056       result = a2 / b2;
0057
0058       // condition for interruption
0059       if ((fabs(result) * boost::math::tools::epsilon<T>()) > fabs(result - prev_result))
0060         break;
0061
0062       b0 = b1; b1 = b2;
0063       a0 = a1; a1 = a2;
0064
0065       ct1 += 2;
0066       xi1 += 1;
0067     }
0068
0069     return a2 / b2;
0070   }
0071
0072   // Luke: C -------- SUBROUTINE R2F1P(BP, CP, Z, A, B, N) --------
0073   // Luke: C ---- PADE APPROXIMATION OF 2F1( 1 , BP; CP ; -Z ) ----
0074   template <class T, class Policy>
0075   inline T hypergeometric_2F1_pade(const T& bp, const T& cp, const T& zp, const Policy&)
0076   {
0077     BOOST_MATH_STD_USING
0078
0079     static const T one = T(1);
0080
0081     // Luke: C ---------- INITIALIZATION -----------
0082     const T z = -zp;
0083     const T zz = z * z;
0084     T b0 = one;
0085     T a0 = one;
0086     T xi1 = one;
0087     T ct1 = cp;
0088     const T b1c1 = (cp - one) * (bp - one);
0089
0090     T b1 = one + ((z / (cp + one)) * (bp + one));
0091     T a1 = b1 - ((bp / cp) * z);
0092
0093     const unsigned max_iterations = boost::math::policies::get_max_series_iterations<Policy>();
0094
0095     T b2 = T(0), a2 = T(0);
0096     T result = T(0), prev_result = a1 / b1;
0097
0098     for (unsigned k = 1; k < max_iterations; ++k)
0099     {
0100       // Luke: C ----- CALCULATION OF THE MULTIPLIERS -----
0101       // Luke: C ----------- FOR THE RECURSION ------------
0102       const T ct2 = ct1 + xi1;
0103       const T ct3 = ct2 * ct2;
0104       const T g2 = (((((ct1 / ct3) * (bp - ct1)) / (ct3 - one)) * xi1) * (bp + xi1)) * zz;
0105       ++xi1;
0106       const T g1 = one + (((((xi1 + xi1) * ct1) + b1c1) / (ct3 + ct2 + ct2)) * z);
0107
0108       // Luke: C ------- THE RECURRENCE RELATIONS ---------
0109       // Luke: C ------------ ARE AS FOLLOWS --------------
0110       b2 = (g1 * b1) + (g2 * b0);
0111       a2 = (g1 * a1) + (g2 * a0);
0112
0113       prev_result = result;
0114       result = a2 / b2;
0115
0116       // condition for interruption
0117       if ((fabs(result) * boost::math::tools::epsilon<T>()) > fabs(result - prev_result))
0118         break;
0119
0120       b0 = b1; b1 = b2;
0121       a0 = a1; a1 = a2;
0122
0123       ++ct1;
0124     }
0125
0126     return a2 / b2;
0127   }
0128
0129   } } } // namespaces
0130
0131 #endif // BOOST_MATH_HYPERGEOMETRIC_PADE_HPP