boost/math/special_functions/detail/hypergeometric_rational.hpp

   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>
  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>
  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