1 /////////////////////////////////////////////////////////////////////////////// 2 // Copyright 2014 Anton Bikineev 3 // Copyright 2014 Christopher Kormanyos 4 // Copyright 2014 John Maddock 5 // Copyright 2014 Paul Bristow 6 // Distributed under the Boost 7 // Software License, Version 1.0. (See accompanying file 8 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 9 // 10 #ifndef BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP 11 #define BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP 12 13 #include <boost/array.hpp> 14 15 namespace boost{ namespace math{ namespace detail{ 16 17 // Luke: C ------- SUBROUTINE R1F1P(AP, CP, Z, A, B, N) --------- 18 // Luke: C --- RATIONAL APPROXIMATION OF 1F1( AP ; CP ; -Z ) ---- 19 template <class T, class Policy> hypergeometric_1F1_rational(const T & ap,const T & cp,const T & zp,const Policy &)20 inline T hypergeometric_1F1_rational(const T& ap, const T& cp, const T& zp, const Policy& ) 21 { 22 BOOST_MATH_STD_USING 23 24 static const T zero = T(0), one = T(1), two = T(2), three = T(3); 25 26 // Luke: C ------------- INITIALIZATION ------------- 27 const T z = -zp; 28 const T z2 = z / two; 29 30 T ct1 = ap * (z / cp); 31 T ct2 = z2 / (one + cp); 32 T xn3 = zero; 33 T xn2 = one; 34 T xn1 = two; 35 T xn0 = three; 36 37 T b1 = one; 38 T a1 = one; 39 T b2 = one + ((one + ap) * (z2 / cp)); 40 T a2 = b2 - ct1; 41 T b3 = one + ((two + b2) * (((two + ap) / three) * ct2)); 42 T a3 = b3 - ((one + ct2) * ct1); 43 ct1 = three; 44 45 const unsigned max_iterations = boost::math::policies::get_max_series_iterations<Policy>(); 46 47 T a4 = T(0), b4 = T(0); 48 T result = T(0), prev_result = a3 / b3; 49 50 for (unsigned k = 2; k < max_iterations; ++k) 51 { 52 // Luke: C ----- CALCULATION OF THE MULTIPLIERS ----- 53 // Luke: C ----------- FOR THE RECURSION ------------ 54 ct2 = (z2 / ct1) / (cp + xn1); 55 const T g1 = one + (ct2 * (xn2 - ap)); 56 ct2 *= ((ap + xn1) / (cp + xn2)); 57 const T g2 = ct2 * ((cp - xn1) + (((ap + xn0) / (ct1 + two)) * z2)); 58 const T g3 = ((ct2 * z2) * (((z2 / ct1) / (ct1 - two)) * ((ap + xn2)) / (cp + xn3))) * (ap - xn2); 59 60 // Luke: C ------- THE RECURRENCE RELATIONS --------- 61 // Luke: C ------------ ARE AS FOLLOWS -------------- 62 b4 = (g1 * b3) + (g2 * b2) + (g3 * b1); 63 a4 = (g1 * a3) + (g2 * a2) + (g3 * a1); 64 65 prev_result = result; 66 result = a4 / b4; 67 68 // condition for interruption 69 if ((fabs(result) * boost::math::tools::epsilon<T>()) > fabs(result - prev_result) / fabs(result)) 70 break; 71 72 b1 = b2; b2 = b3; b3 = b4; 73 a1 = a2; a2 = a3; a3 = a4; 74 75 xn3 = xn2; 76 xn2 = xn1; 77 xn1 = xn0; 78 xn0 += 1; 79 ct1 += two; 80 } 81 82 return result; 83 } 84 85 // Luke: C ----- SUBROUTINE R2F1P(AB, BP, CP, Z, A, B, N) ------- 86 // Luke: C -- RATIONAL APPROXIMATION OF 2F1( AB , BP; CP ; -Z ) - 87 template <class T, class Policy> hypergeometric_2F1_rational(const T & ap,const T & bp,const T & cp,const T & zp,const unsigned n,const Policy &)88 inline T hypergeometric_2F1_rational(const T& ap, const T& bp, const T& cp, const T& zp, const unsigned n, const Policy& ) 89 { 90 BOOST_MATH_STD_USING 91 92 static const T one = T(1), two = T(2), three = T(3), four = T(4), 93 six = T(6), half_7 = T(3.5), half_3 = T(1.5), forth_3 = T(0.75); 94 95 // Luke: C ------------- INITIALIZATION ------------- 96 const T z = -zp; 97 const T z2 = z / two; 98 99 T sabz = (ap + bp) * z; 100 const T ab = ap * bp; 101 const T abz = ab * z; 102 const T abz1 = z + (abz + sabz); 103 const T abz2 = abz1 + (sabz + (three * z)); 104 const T cp1 = cp + one; 105 const T ct1 = cp1 + cp1; 106 107 T b1 = one; 108 T a1 = one; 109 T b2 = one + (abz1 / (cp + cp)); 110 T a2 = b2 - (abz / cp); 111 T b3 = one + ((abz2 / ct1) * (one + (abz1 / ((-six) + (three * ct1))))); 112 T a3 = b3 - ((abz / cp) * (one + ((abz2 - abz1) / ct1))); 113 sabz /= four; 114 115 const T abz1_div_4 = abz1 / four; 116 const T cp1_inc = cp1 + one; 117 const T cp1_mul_cp1_inc = cp1 * cp1_inc; 118 119 boost::array<T, 9u> d = {{ 120 ((half_7 - ab) * z2) - sabz, 121 abz1_div_4, 122 abz1_div_4 - (two * sabz), 123 cp1_inc, 124 cp1_mul_cp1_inc, 125 cp * cp1_mul_cp1_inc, 126 half_3, 127 forth_3, 128 forth_3 * z 129 }}; 130 131 T xi = three; 132 T a4 = T(0), b4 = T(0); 133 for (unsigned k = 2; k < n; ++k) 134 { 135 // Luke: C ----- CALCULATION OF THE MULTIPLIERS ----- 136 // Luke: C ----------- FOR THE RECURSION ------------ 137 T g3 = (d[2] / d[7]) * (d[1] / d[5]); 138 d[1] += d[8] + sabz; 139 d[2] += d[8] - sabz; 140 g3 *= d[1] / d[6]; 141 T g1 = one + (((d[1] + d[0]) / d[6]) / d[3]); 142 T g2 = (d[1] / d[4]) / d[6]; 143 d[7] += two * d[6]; 144 ++d[6]; 145 g2 *= cp1 - (xi + ((d[2] + d[0]) / d[6])); 146 147 // Luke: C ------- THE RECURRENCE RELATIONS --------- 148 // Luke: C ------------ ARE AS FOLLOWS -------------- 149 b4 = (g1 * b3) + (g2 * b2) + (g3 * b1); 150 a4 = (g1 * a3) + (g2 * a2) + (g3 * a1); 151 b1 = b2; b2 = b3; b3 = b4; 152 a1 = a2; a2 = a3; a3 = a4; 153 154 d[8] += z2; 155 d[0] += two * d[8]; 156 d[5] += three * d[4]; 157 d[4] += two * d[3]; 158 ++d[3]; 159 ++xi; 160 } 161 162 return a4 / b4; 163 } 164 165 } } } // namespaces 166 167 #endif // BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP 168