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)
10 #ifndef BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP
11 #define BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP
15 namespace boost{ namespace math{ namespace detail{
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>
20 inline T hypergeometric_1F1_rational(const T& ap, const T& cp, const T& zp, const Policy& )
24 static const T zero = T(0), one = T(1), two = T(2), three = T(3);
26 // Luke: C ------------- INITIALIZATION -------------
30 T ct1 = ap * (z / cp);
31 T ct2 = z2 / (one + cp);
39 T b2 = one + ((one + ap) * (z2 / cp));
41 T b3 = one + ((two + b2) * (((two + ap) / three) * ct2));
42 T a3 = b3 - ((one + ct2) * ct1);
45 const unsigned max_iterations = boost::math::policies::get_max_series_iterations<Policy>();
47 T a4 = T(0), b4 = T(0);
48 T result = T(0), prev_result = a3 / b3;
50 for (unsigned k = 2; k < max_iterations; ++k)
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);
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);
68 // condition for interruption
69 if ((fabs(result) * boost::math::tools::epsilon<T>()) > fabs(result - prev_result) / fabs(result))
72 b1 = b2; b2 = b3; b3 = b4;
73 a1 = a2; a2 = a3; a3 = a4;
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>
88 inline T hypergeometric_2F1_rational(const T& ap, const T& bp, const T& cp, const T& zp, const unsigned n, const Policy& )
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);
95 // Luke: C ------------- INITIALIZATION -------------
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;
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)));
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;
119 std::array<T, 9u> d = {{
120 ((half_7 - ab) * z2) - sabz,
122 abz1_div_4 - (two * sabz),
125 cp * cp1_mul_cp1_inc,
132 T a4 = T(0), b4 = T(0);
133 for (unsigned k = 2; k < n; ++k)
135 // Luke: C ----- CALCULATION OF THE MULTIPLIERS -----
136 // Luke: C ----------- FOR THE RECURSION ------------
137 T g3 = (d[2] / d[7]) * (d[1] / d[5]);
141 T g1 = one + (((d[1] + d[0]) / d[6]) / d[3]);
142 T g2 = (d[1] / d[4]) / d[6];
145 g2 *= cp1 - (xi + ((d[2] + d[0]) / d[6]));
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;
156 d[5] += three * d[4];
167 #endif // BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP