132 lines · plain
1 2///////////////////////////////////////////////////////////////////////////////3// Copyright 2014 Anton Bikineev4// Copyright 2014 Christopher Kormanyos5// Copyright 2014 John Maddock6// Copyright 2014 Paul Bristow7// Distributed under the Boost8// Software License, Version 1.0. (See accompanying file9// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)10//11#ifndef BOOST_MATH_HYPERGEOMETRIC_PADE_HPP12#define BOOST_MATH_HYPERGEOMETRIC_PADE_HPP13 14 namespace boost{ namespace math{ namespace detail{15 16 // Luke: C ---------- SUBROUTINE R1F1P(CP, Z, A, B, N) ----------17 // Luke: C ----- PADE APPROXIMATION OF 1F1( 1 ; CP ; -Z ) -------18 template <class T, class Policy>19 inline T hypergeometric_1F1_pade(const T& cp, const T& zp, const Policy& )20 {21 BOOST_MATH_STD_USING22 23 static const T one = T(1);24 25 // Luke: C ------------- INITIALIZATION -------------26 const T z = -zp;27 const T zz = z * z;28 T b0 = one;29 T a0 = one;30 T xi1 = one;31 T ct1 = cp + one;32 T cp1 = cp - one;33 34 T b1 = one + (z / ct1);35 T a1 = b1 - (z / cp);36 37 const unsigned max_iterations = boost::math::policies::get_max_series_iterations<Policy>();38 39 T b2 = T(0), a2 = T(0);40 T result = T(0), prev_result;41 42 for (unsigned k = 1; k < max_iterations; ++k)43 {44 // Luke: C ----- CALCULATION OF THE MULTIPLIERS -----45 // Luke: C ----------- FOR THE RECURSION ------------46 const T ct2 = ct1 * ct1;47 const T g1 = one + ((cp1 / (ct2 + ct1 + ct1)) * z);48 const T g2 = ((xi1 / (ct2 - one)) * ((xi1 + cp1) / ct2)) * zz;49 50 // Luke: C ------- THE RECURRENCE RELATIONS ---------51 // Luke: C ------------ ARE AS FOLLOWS --------------52 b2 = (g1 * b1) + (g2 * b0);53 a2 = (g1 * a1) + (g2 * a0);54 55 prev_result = result;56 result = a2 / b2;57 58 // condition for interruption59 if ((fabs(result) * boost::math::tools::epsilon<T>()) > fabs(result - prev_result))60 break;61 62 b0 = b1; b1 = b2;63 a0 = a1; a1 = a2;64 65 ct1 += 2;66 xi1 += 1;67 }68 69 return a2 / b2;70 }71 72 // Luke: C -------- SUBROUTINE R2F1P(BP, CP, Z, A, B, N) --------73 // Luke: C ---- PADE APPROXIMATION OF 2F1( 1 , BP; CP ; -Z ) ----74 template <class T, class Policy>75 inline T hypergeometric_2F1_pade(const T& bp, const T& cp, const T& zp, const Policy&)76 {77 BOOST_MATH_STD_USING78 79 static const T one = T(1);80 81 // Luke: C ---------- INITIALIZATION -----------82 const T z = -zp;83 const T zz = z * z;84 T b0 = one;85 T a0 = one;86 T xi1 = one;87 T ct1 = cp;88 const T b1c1 = (cp - one) * (bp - one);89 90 T b1 = one + ((z / (cp + one)) * (bp + one));91 T a1 = b1 - ((bp / cp) * z);92 93 const unsigned max_iterations = boost::math::policies::get_max_series_iterations<Policy>();94 95 T b2 = T(0), a2 = T(0);96 T result = T(0), prev_result = a1 / b1;97 98 for (unsigned k = 1; k < max_iterations; ++k)99 {100 // Luke: C ----- CALCULATION OF THE MULTIPLIERS -----101 // Luke: C ----------- FOR THE RECURSION ------------102 const T ct2 = ct1 + xi1;103 const T ct3 = ct2 * ct2;104 const T g2 = (((((ct1 / ct3) * (bp - ct1)) / (ct3 - one)) * xi1) * (bp + xi1)) * zz;105 ++xi1;106 const T g1 = one + (((((xi1 + xi1) * ct1) + b1c1) / (ct3 + ct2 + ct2)) * z);107 108 // Luke: C ------- THE RECURRENCE RELATIONS ---------109 // Luke: C ------------ ARE AS FOLLOWS --------------110 b2 = (g1 * b1) + (g2 * b0);111 a2 = (g1 * a1) + (g2 * a0);112 113 prev_result = result;114 result = a2 / b2;115 116 // condition for interruption117 if ((fabs(result) * boost::math::tools::epsilon<T>()) > fabs(result - prev_result))118 break;119 120 b0 = b1; b1 = b2;121 a0 = a1; a1 = a2;122 123 ++ct1;124 }125 126 return a2 / b2;127 }128 129 } } } // namespaces130 131#endif // BOOST_MATH_HYPERGEOMETRIC_PADE_HPP132