]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
1 | |
2 | /////////////////////////////////////////////////////////////////////////////// | |
3 | // Copyright 2014 Anton Bikineev | |
4 | // Copyright 2014 Christopher Kormanyos | |
5 | // Copyright 2014 John Maddock | |
6 | // Copyright 2014 Paul Bristow | |
7 | // Distributed under the Boost | |
8 | // Software License, Version 1.0. (See accompanying file | |
9 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
10 | // | |
11 | #ifndef BOOST_MATH_HYPERGEOMETRIC_PADE_HPP | |
12 | #define BOOST_MATH_HYPERGEOMETRIC_PADE_HPP | |
13 | ||
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_USING | |
22 | ||
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 interruption | |
59 | 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& pol) | |
76 | { | |
77 | BOOST_MATH_STD_USING | |
78 | ||
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 interruption | |
117 | 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 | } } } // namespaces | |
130 | ||
131 | #endif // BOOST_MATH_HYPERGEOMETRIC_PADE_HPP |