ceph/src/boost/boost/math/special_functions/detail/hypergeometric_pade.hpp

   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