Import SOC headers, Test and fix up 1F0 and 0F1.

include/boost/math/special_functions/detail/hypergeometric_0f1_bessel.hpp

@@ -0,0 +1,45 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+#ifndef BOOST_MATH_HYPERGEOMETRIC_0F1_BESSEL_HPP
+  #define BOOST_MATH_HYPERGEOMETRIC_0F1_BESSEL_HPP
+
+  namespace boost { namespace math { namespace detail {
+
+  template <class T, class Policy>
+  inline T hypergeometric_0f1_bessel(const T& b, const T& z, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+
+    const bool is_z_nonpositive = z <= 0;
+
+    const T sqrt_z = is_z_nonpositive ? T(sqrt(-z)) : T(sqrt(z));
+    const T bessel_mult = is_z_nonpositive ?
+      boost::math::cyl_bessel_j(b - 1, 2 * sqrt_z, pol) :
+      boost::math::cyl_bessel_i(b - 1, 2 * sqrt_z, pol) ;
+
+    if (b > boost::math::max_factorial<T>::value)
+    {
+       const T lsqrt_z = log(sqrt_z);
+       const T lsqrt_z_pow_b = (b - 1) * lsqrt_z;
+       T lg = (boost::math::lgamma(b, pol) - lsqrt_z_pow_b);
+       lg = exp(lg);
+       return lg * bessel_mult;
+    }
+    else
+    {
+       const T sqrt_z_pow_b = pow(sqrt_z, b - 1);
+       return (boost::math::tgamma(b, pol) / sqrt_z_pow_b) * bessel_mult;
+    }
+  }
+
+  } } } // namespaces
+
+#endif // BOOST_MATH_HYPERGEOMETRIC_0F1_BESSEL_HPP

include/boost/math/special_functions/detail/hypergeometric_1f1_bessel.hpp

@@ -0,0 +1,312 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+#ifndef BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP
+  #define BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP
+
+  #include <boost/math/tools/series.hpp>
+  #include <boost/math/special_functions/bessel.hpp>
+  #include <boost/math/special_functions/laguerre.hpp>
+
+  namespace boost { namespace math { namespace detail {
+
+  // declarations of helpers for 13_3_7 and 13_3_8:
+  // bessel_j recurrence relation based on finite continued fractions;
+  // it is stable forward recurrence algorithm, see William J. Lentz
+  // "Generating Bessel functions in Mie scattering calculations
+  // using continued fractions" for details
+  template <class T>
+  inline T hypergeometric_bessel_j_recurrence_next(const T& jvm1, const T& v, const T& z);
+
+  // swap values of Bessel functions:
+  template <class T>
+  inline void hypergeometric_bessel_j_recurrence_iterate(T& jvm1, T& jv, const T& v, const T& z);
+
+  // next coefficient for 13_3_7
+  template <class T>
+  inline T hypergeometric_13_3_7_coefficient_next(const T& anm3, const T& anm2, const T& a, const T& b, const unsigned n);
+
+  // next coefficient for 13_3_8
+  template <class T>
+  inline T hypergeometric_13_3_8_coefficient_next(const T& cnm3, const T& cnm2, const T& cnm1, const T& a, const T& b, const unsigned n);
+
+  // swap values of coefficients:
+  template <class T>
+  inline void hypergeometric_coefficient_13_3_7_iterate(T& anm3, T& anm2, T& anm1, T& an, const T& a, const T& b, const unsigned n);
+
+  template <class T>
+  inline void hypergeometric_coefficient_13_3_8_iterate(T& cnm3, T& cnm2, T& cnm1, T& cn, const T& a, const T& b, const unsigned n);
+
+  // term class of Abramowitz & Stegun 13_3_7 formula
+  template <class T>
+  struct hypergeometric_1f1_13_3_7_series_term
+  {
+    typedef T result_type;
+
+    hypergeometric_1f1_13_3_7_series_term(const T& a, const T& b, const T& z):
+      a(a), b(b), z(z), n(0u)
+    {
+      BOOST_MATH_STD_USING
+
+      sqrt_z_pow_n = sqrt_z = sqrt(z);
+      sqrt_2b_minus_4a_pow_n = sqrt_2b_minus_4a = sqrt(2 * (b - (2 * a)));
+      sqrt_2zb_minus_4za = sqrt_z * sqrt_2b_minus_4a;
+
+      v_current = b - 1;
+
+      anm3 = 1;
+      anm2 = 0;
+      anm1 = b / 2;
+      an = detail::hypergeometric_13_3_7_coefficient_next(anm3, anm2, a, b, 3u);
+
+      jvm1 = boost::math::cyl_bessel_j(v_current, sqrt_2zb_minus_4za);
+      v_current++;
+      jv = detail::hypergeometric_bessel_j_recurrence_next(jvm1, v_current, sqrt_2zb_minus_4za);
+
+      term = jvm1;
+      ++n;
+    }
+
+    T operator()()
+    {
+      const T result = term;
+
+      iterate();
+      term = ((anm2 * sqrt_z_pow_n) / sqrt_2b_minus_4a_pow_n) * jv;
+
+      return result;
+    }
+
+  private:
+    void iterate()
+    {
+      ++n; ++v_current;
+      sqrt_z_pow_n *= sqrt_z;
+      sqrt_2b_minus_4a_pow_n *= sqrt_2b_minus_4a;
+
+      detail::hypergeometric_coefficient_13_3_7_iterate(anm3, anm2, anm1, an, a, b, n + 2);
+      detail::hypergeometric_bessel_j_recurrence_iterate(jvm1, jv, v_current, sqrt_2zb_minus_4za);
+    }
+
+    const T a, b, z;
+    unsigned n;
+    T sqrt_z, sqrt_2b_minus_4a, sqrt_2zb_minus_4za;
+    T sqrt_z_pow_n, sqrt_2b_minus_4a_pow_n;
+    T v_current;
+    T anm3, anm2, anm1, an;
+    T jvm1, jv;
+    T term;
+  };
+
+  // function for 13_3_7 evaluation
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_13_3_7_series(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+
+    const T sqrt_bz_div_2_minus_az = sqrt(((b * z) / 2) - (a * z));
+    const T prefix = ((boost::math::tgamma(b, pol) * sqrt_bz_div_2_minus_az) /
+        pow(sqrt_bz_div_2_minus_az, b)) * exp(z / 2);
+
+    detail::hypergeometric_1f1_13_3_7_series_term<T> s(a, b, z);
+    boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
+#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+    T zero = 0;
+    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
+#else
+    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
+#endif
+    boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_13_3_7_series<%1%>(%1%,%1%,%1%)", max_iter, pol);
+    return prefix * result;
+  }
+
+  // term class of Abramowitz & Stegun 13_3_8 formula
+  template <class T>
+  struct hypergeometric_1f1_13_3_8_series_term
+  {
+    typedef T result_type;
+
+    static const T h;
+
+    hypergeometric_1f1_13_3_8_series_term(const T& a, const T& b, const T& z):
+      a(a), b(b), z(z), n(0u)
+    {
+      BOOST_MATH_STD_USING
+
+      sqrt_minus_az = sqrt(-a * z);
+      const T double_sqrt_minus_az = 2 * sqrt_minus_az;
+
+      v_current = b - 1;
+      z_pow_n = z; sqrt_minus_az_pow_n = sqrt_minus_az;
+
+      cnm3 = 1;
+      cnm2 = -b * h;
+      cnm1 = (((1 - (2 * h)) * a) + ((b * (b + 1)) * (h * h))) / 2;
+      cn = detail::hypergeometric_13_3_8_coefficient_next(cnm3, cnm2, cnm1, a, b, 3u);
+
+      jvm1 = boost::math::cyl_bessel_j(v_current, double_sqrt_minus_az);
+      ++v_current;
+      jv = detail::hypergeometric_bessel_j_recurrence_next(jvm1, b, double_sqrt_minus_az);
+    }
+
+    T operator()()
+    {
+      ++n;
+      switch (n - 1)
+      {
+        case 0u:
+          return jvm1;
+        case 1u:
+          detail::hypergeometric_bessel_j_recurrence_iterate(jvm1, jv, v_current, T(2 * sqrt_minus_az));
+          return ((cnm2 * z) / sqrt_minus_az) * jvm1;
+      }
+
+      ++v_current;
+      z_pow_n *= z;
+      sqrt_minus_az_pow_n *= sqrt_minus_az;
+      const T result = ((cnm1 * z_pow_n) / sqrt_minus_az_pow_n) * jv;
+
+      detail::hypergeometric_coefficient_13_3_8_iterate(cnm3, cnm2, cnm1, cn, a, b, n + 1);
+      detail::hypergeometric_bessel_j_recurrence_iterate(jvm1, jv, v_current, T(2 * sqrt_minus_az));
+      return result;
+    }
+
+  private:
+    const T a, b, z;
+    unsigned n;
+    T sqrt_minus_az;
+    T z_pow_n, sqrt_minus_az_pow_n;
+    T v_current;
+    T cnm3, cnm2, cnm1, cn;
+    T jvm1, jv;
+  };
+
+  template <class T>
+  const T hypergeometric_1f1_13_3_8_series_term<T>::h = -boost::math::constants::pi<T>() / T(10.);
+
+  // function for 13_3_8 evaluation
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_13_3_8_series(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+
+    const T sqrt_minus_az_pow_b_minus_one = pow(sqrt(-a * z), b - 1);
+    const T prefix = (boost::math::tgamma(b, pol) / sqrt_minus_az_pow_b_minus_one) *
+                      exp(detail::hypergeometric_1f1_13_3_8_series_term<T>::h * z);
+
+    detail::hypergeometric_1f1_13_3_8_series_term<T> s(a, b, z);
+    boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
+#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+    T zero = 0;
+    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
+#else
+    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
+#endif
+    boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_13_3_8_series<%1%>(%1%,%1%,%1%)", max_iter, pol);
+    return prefix * result;
+  }
+
+  // definitions of helpers:
+
+  template <class T>
+  inline T hypergeometric_bessel_j_recurrence_next(const T& jvm1, const T& v, const T& z)
+  {
+    BOOST_MATH_STD_USING
+
+    // TODO: we can estimate a possible capacity for these vectors
+    std::vector<T> finite_cf_nums;
+    std::vector<T> finite_cf_denoms;
+
+    const T a_coef_when_i_is_one = (2 * v) / z;
+    const T a_coef_when_i_is_two = (-2 * (v + 1)) / z;
+
+    finite_cf_nums.push_back(a_coef_when_i_is_one);
+    finite_cf_nums.push_back(a_coef_when_i_is_two + (1. / a_coef_when_i_is_one));
+
+    finite_cf_denoms.push_back(a_coef_when_i_is_two);
+
+    unsigned i = 3u;
+
+    do {
+      const T temp = (2 * (v + (i - 1))) / z;
+      const T next_coef = i & 1 ? temp : T(-temp);
+
+      finite_cf_nums.push_back(next_coef + (1. / finite_cf_nums.back()));
+      finite_cf_denoms.push_back(next_coef + (1. / finite_cf_denoms.back()));
+
+      ++i;
+    } while (fabs(finite_cf_nums.back() - finite_cf_denoms.back()) > boost::math::tools::epsilon<T>());
+
+    T ratio = finite_cf_nums.front();
+
+    for (i = 0u; i < finite_cf_denoms.size(); ++i)
+      ratio *= finite_cf_nums[i+1] / finite_cf_denoms[i];
+
+    return jvm1 / ratio;
+  }
+
+  template <class T>
+  inline void hypergeometric_bessel_j_recurrence_iterate(T& jvm1, T& jv, const T& v, const T& z)
+  {
+    using std::swap;
+
+    swap(jvm1, jv);
+    jv = detail::hypergeometric_bessel_j_recurrence_next(jvm1, v, z);
+  }
+
+  template <class T>
+  inline T hypergeometric_13_3_7_coefficient_next(const T& anm3, const T& anm2, const T& a, const T& b, const unsigned n)
+  {
+    const T term_anm2 = ((n + b) - 2) * anm2;
+    const T term_anm3 = ((2 * a) - b) * anm3;
+
+    return (term_anm2 + term_anm3) / n;
+
+  }
+
+  template <class T>
+  inline T hypergeometric_13_3_8_coefficient_next(const T& cnm3, const T& cnm2, const T& cnm1, const T& a, const T& b, const unsigned n)
+  {
+    static const T& h = detail::hypergeometric_1f1_13_3_8_series_term<T>::h;
+    const T one_minus_two_h = 1 - (2 * h);
+    const T h_minus_one = h - 1;
+
+    const T term_cnm1  = ((one_minus_two_h * n) - (b * h)) * cnm1;
+    const T term_cnm2 = ((one_minus_two_h * a) - ((h * h_minus_one) * (b + (n - 1)))) * cnm2;
+    const T term_cnm3 = ((-h * h_minus_one) * a) * cnm3;
+
+    return ((term_cnm1 + term_cnm2) + term_cnm3) / (n + 1);
+  }
+
+  template <class T>
+  inline void hypergeometric_coefficient_13_3_7_iterate(T& anm3, T& anm2, T& anm1, T& an, const T& a, const T& b, const unsigned n)
+  {
+    using std::swap;
+
+    swap(anm3, anm2);
+    swap(anm2, anm1);
+    swap(anm1, an);
+    an = detail::hypergeometric_13_3_7_coefficient_next(anm3, anm2, a, b, n);
+  }
+
+  template <class T>
+  inline void hypergeometric_coefficient_13_3_8_iterate(T& cnm3, T& cnm2, T& cnm1, T& cn, const T& a, const T& b, const unsigned n)
+  {
+    using std::swap;
+
+    swap(cnm3, cnm2);
+    swap(cnm2, cnm1);
+    swap(cnm1, cn);
+    cn = detail::hypergeometric_13_3_8_coefficient_next(cnm3, cnm2, cnm1, a, b, n);
+  }
+
+  } } } // namespaces
+
+#endif // BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP

include/boost/math/special_functions/detail/hypergeometric_1f1_recurrence.hpp

@@ -0,0 +1,200 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_HYPERGEOMETRIC_1F1_RECURRENCE_HPP_
+  #define BOOST_HYPERGEOMETRIC_1F1_RECURRENCE_HPP_
+
+  #include <boost/math/special_functions/modf.hpp>
+  #include <boost/math/special_functions/next.hpp>
+
+  #include <boost/math/tools/recurrence.hpp>
+
+  namespace boost { namespace math { namespace detail {
+
+  template <class T>
+  struct hypergeometric_1f1_recurrence_a_coefficients
+  {
+    typedef boost::math::tuple<T, T, T> result_type;
+
+    hypergeometric_1f1_recurrence_a_coefficients(const T& a, const T& b, const T& z):
+    a(a), b(b), z(z)
+    {
+    }
+
+    result_type operator()(boost::intmax_t i) const
+    {
+      const T ai = a + i;
+
+      const T an = -ai;
+      const T bn = (b - (2 * ai)) - z;
+      const T cn = b - ai;
+
+      return boost::math::make_tuple(an, bn, cn);
+    }
+
+  private:
+    const T a, b, z;
+  };
+
+  template <class T>
+  struct hypergeometric_1f1_recurrence_b_coefficients
+  {
+    typedef boost::math::tuple<T, T, T> result_type;
+
+    hypergeometric_1f1_recurrence_b_coefficients(const T& a, const T& b, const T& z):
+    a(a), b(b), z(z)
+    {
+    }
+
+    result_type operator()(boost::intmax_t i) const
+    {
+      const T bi = b + i;
+
+      const T an = z * (bi - a);
+      const T bn = bi * ((z + bi) - 1);
+      const T cn = bi * (bi - 1);
+
+      return boost::math::make_tuple(an, bn, cn);
+    }
+
+  private:
+    const T a, b, z;
+  };
+
+  template <class T>
+  struct hypergeometric_1f1_recurrence_a_and_b_coefficients
+  {
+    typedef boost::math::tuple<T, T, T> result_type;
+
+    hypergeometric_1f1_recurrence_a_and_b_coefficients(const T& a, const T& b, const T& z):
+    a(a), b(b), z(z)
+    {
+    }
+
+    result_type operator()(boost::intmax_t i) const
+    {
+      const T ai = a + i;
+      const T bi = b + i;
+
+      const T an = ai * z;
+      const T bn = bi * ((1 - bi) + z);
+      const T cn = bi * (1 - bi);
+
+      return boost::math::make_tuple(an, bn, cn);
+    }
+
+  private:
+    const T a, b, z;
+  };
+
+  // forward declaration for initial values
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_imp(const T& a, const T& b, const T& z, const Policy& pol);
+
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_backward_recurrence_for_negative_a(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING // modf, frexp, fabs, pow
+
+    boost::intmax_t integer_part = 0;
+    const T bk = modf(b, &integer_part);
+    T ak = modf(a, &integer_part);
+
+    int exp_of_a = 0; frexp(a, &exp_of_a);
+    int exp_of_b = 0; frexp(b, &exp_of_b);
+
+    const bool are_fractional_parts_close_enough =
+      fabs(boost::math::float_distance(ak, bk)) <= pow(2, (std::max)(exp_of_a, exp_of_b));
+
+    if ((a < b) && (b < 0) && (are_fractional_parts_close_enough)) // TODO: has to be researched deeper
+    {
+      ak = b - 1;
+      integer_part -= (boost::math::lltrunc(ceil(b)) - 1);
+    }
+
+    T first = detail::hypergeometric_1f1_imp(ak, b, z, pol);
+    --ak;
+    T second = detail::hypergeometric_1f1_imp(ak, b, z, pol);
+
+    detail::hypergeometric_1f1_recurrence_a_coefficients<T> s(ak, b, z);
+
+    return tools::solve_recurrence_relation_backward(s,
+                                                     static_cast<unsigned int>(std::abs(integer_part)),
+                                                     first,
+                                                     second);
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_forward_recurrence_for_positive_a(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING // modf, fabs
+
+    boost::intmax_t integer_part = 0;
+    T ak = modf(a, &integer_part);
+
+    T first = detail::hypergeometric_1f1_imp(ak, b, z, pol);
+    ++ak;
+    T second = detail::hypergeometric_1f1_imp(ak, b, z, pol);
+
+    detail::hypergeometric_1f1_recurrence_a_coefficients<T> s(ak, b, z);
+
+    return tools::solve_recurrence_relation_forward(s, integer_part, first, second);
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_backward_recurrence_for_negative_b(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING // modf, fabs
+
+    boost::intmax_t integer_part = 0;
+    T bk = modf(b, &integer_part);
+
+    T first = detail::hypergeometric_1f1_imp(a, bk, z, pol);
+    --bk;
+    T second = detail::hypergeometric_1f1_imp(a, bk, z, pol);
+
+    detail::hypergeometric_1f1_recurrence_b_coefficients<T> s(a, bk, z);
+
+    return tools::solve_recurrence_relation_backward(s,
+                                                     static_cast<unsigned int>(std::abs(integer_part)),
+                                                     first,
+                                                     second);
+  }
+
+  // this method works provided that integer part of a is the same as integer part of b
+  // we don't make this check inside it.
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_backward_recurrence_for_negative_a_and_b(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING // modf, fabs
+
+    boost::intmax_t integer_part = 0;
+    T ak = modf(a, &integer_part);
+    T bk = modf(b, &integer_part);
+
+    T first = detail::hypergeometric_1f1_imp(ak, bk, z, pol);
+    --ak; --bk;
+    T second = detail::hypergeometric_1f1_imp(ak, bk, z, pol);
+
+    detail::hypergeometric_1f1_recurrence_a_and_b_coefficients<T> s(ak, bk, z);
+
+    return tools::solve_recurrence_relation_backward(s, fabs(integer_part), first, second);
+  }
+
+  // ranges
+  template <class T>
+  inline bool hypergeometric_1f1_is_a_small_enough(const T& a)
+  {
+    return a < -10; // TODO: make dependent on precision
+  }
+
+  } } } // namespaces
+
+#endif // BOOST_HYPERGEOMETRIC_1F1_RECURRENCE_HPP_

include/boost/math/special_functions/detail/hypergeometric_asym.hpp

@@ -0,0 +1,87 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+#ifndef BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP
+  #define BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP
+
+  #include <boost/math/special_functions/gamma.hpp>
+
+  namespace boost { namespace math {
+
+  // forward declaration of used function 2f0
+  template <class T1, class T2, class T3, class Policy>
+  inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_2f0(T1 a1, T2 a2, T3 z, const Policy& /* pol */);
+
+  namespace detail {
+
+  // assumes a and b are not non-positive integers
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_asym_positive_series(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+    static const char* const function = "boost::math::hypergeometric_1f1_asym_positive_series<%1%>(%1%,%1%,%1%)";
+
+    const bool is_a_integer = (a == floor(a));
+    const bool is_b_integer = (b == floor(b));
+
+    BOOST_ASSERT(!(is_a_integer && a <= 0) && !(is_b_integer && b <= 0));
+
+    const T gamma_ratio = (b > 0 && a > 0) ?
+      boost::math::tgamma_ratio(b, a, pol) :
+      T(boost::math::tgamma(b, pol) / boost::math::tgamma(a, pol));
+    const T prefix_a = (exp(z) * gamma_ratio) * pow(z, (a - b));
+
+    return prefix_a * boost::math::hypergeometric_2f0(b - a, 1 - a, 1 / z, pol);
+  }
+
+  // assumes b and (b - a) are not non-positive integers
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_asym_negative_series(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+    static const char* const function = "boost::math::hypergeometric_1f1_asym_negative_series<%1%>(%1%,%1%,%1%)";
+
+    const T b_minus_a = b - a;
+
+    const bool is_a_integer = (a == floor(a));
+    const bool is_b_minus_a_integer = (b_minus_a == floor(b_minus_a));
+
+    BOOST_ASSERT(!(is_a_integer && a <= 0) && !(is_b_minus_a_integer && b_minus_a <= 0));
+
+    const T gamma_ratio = (b > 0 && b_minus_a > 0) ?
+      boost::math::tgamma_ratio(b, b_minus_a, pol) :
+      boost::math::tgamma(b) / boost::math::tgamma((b_minus_a), pol);
+    const T prefix_b = gamma_ratio / pow(-z, a);
+
+    return prefix_b * boost::math::hypergeometric_2f0(a, 1 - b_minus_a, -1 / z, pol);
+  }
+
+  // experimental range
+  template <class T>
+  inline bool hypergeometric_1f1_asym_region(const T& a, const T& b, const T& z)
+  {
+    BOOST_MATH_STD_USING
+
+    const T the_max_of_one_and_b_minus_a  ((std::max)(T(1), fabs(b - a)));
+    const T the_max_of_one_and_one_minus_a((std::max)(T(1), fabs(1 - a)));
+
+    const T the_product_of_these_maxima(the_max_of_one_and_b_minus_a * the_max_of_one_and_one_minus_a);
+
+    const T abs_of_z(fabs(z));
+
+    if ((abs_of_z > 100) && (the_product_of_these_maxima < (abs_of_z / 2))) // TODO: not a good way
+      return true;
+
+    return false;
+  }
+
+  } } } // namespaces
+
+#endif // BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP

include/boost/math/special_functions/detail/hypergeometric_cf.hpp

@@ -0,0 +1,228 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+#ifndef BOOST_MATH_HYPERGEOMETRIC_CF_HPP
+  #define BOOST_MATH_HYPERGEOMETRIC_CF_HPP
+
+  namespace boost { namespace math { namespace detail {
+
+  // primary template for term of continued fraction
+  template <class T, unsigned p, unsigned q>
+  struct hypergeometric_pfq_cf_term;
+
+  // partial specialization for 0f1
+  template <class T>
+  struct hypergeometric_pfq_cf_term<T, 0u, 1u>
+  {
+    typedef std::pair<T,T> result_type;
+
+    hypergeometric_pfq_cf_term(const T& b, const T& z):
+      n(1), b(b), z(z),
+      term(std::make_pair(T(0), T(1)))
+    {
+    }
+
+    result_type operator()()
+    {
+      const result_type result = term;
+      ++b; ++n;
+      numer = -(z / (b * n));
+      term = std::make_pair(numer, 1 - numer);
+      return result;
+    }
+
+  private:
+    unsigned n;
+    T b;
+    const T z;
+    T numer;
+    result_type term;
+  };
+
+  // partial specialization for 1f0
+  template <class T>
+  struct hypergeometric_pfq_cf_term<T, 1u, 0u>
+  {
+    typedef std::pair<T,T> result_type;
+
+    hypergeometric_pfq_cf_term(const T& a, const T& z):
+      n(1), a(a), z(z),
+      term(std::make_pair(T(0), T(1)))
+    {
+    }
+
+    result_type operator()()
+    {
+      const result_type result = term;
+      ++a; ++n;
+      numer = -((a * z) / n);
+      term = std::make_pair(numer, 1 - numer);
+      return result;
+    }
+
+  private:
+    unsigned n;
+    T a;
+    const T z;
+    T numer;
+    result_type term;
+  };
+
+  // partial specialization for 1f1
+  template <class T>
+  struct hypergeometric_pfq_cf_term<T, 1u, 1u>
+  {
+    typedef std::pair<T,T> result_type;
+
+    hypergeometric_pfq_cf_term(const T& a, const T& b, const T& z):
+      n(1), a(a), b(b), z(z),
+      term(std::make_pair(T(0), T(1)))
+    {
+    }
+
+    result_type operator()()
+    {
+      const result_type result = term;
+      ++a; ++b; ++n;
+      numer = -((a * z) / (b * n));
+      term = std::make_pair(numer, 1 - numer);
+      return result;
+    }
+
+  private:
+    unsigned n;
+    T a, b;
+    const T z;
+    T numer;
+    result_type term;
+  };
+
+  // partial specialization for 1f2
+  template <class T>
+  struct hypergeometric_pfq_cf_term<T, 1u, 2u>
+  {
+    typedef std::pair<T,T> result_type;
+
+    hypergeometric_pfq_cf_term(const T& a, const T& b, const T& c, const T& z):
+      n(1), a(a), b(b), c(c), z(z),
+      term(std::make_pair(T(0), T(1)))
+    {
+    }
+
+    result_type operator()()
+    {
+      const result_type result = term;
+      ++a; ++b; ++c; ++n;
+      numer = -((a * z) / ((b * c) * n));
+      term = std::make_pair(numer, 1 - numer);
+      return result;
+    }
+
+  private:
+    unsigned n;
+    T a, b, c;
+    const T z;
+    T numer;
+    result_type term;
+  };
+
+  // partial specialization for 2f1
+  template <class T>
+  struct hypergeometric_pfq_cf_term<T, 2u, 1u>
+  {
+    typedef std::pair<T,T> result_type;
+
+    hypergeometric_pfq_cf_term(const T& a, const T& b, const T& c, const T& z):
+      n(1), a(a), b(b), c(c), z(z),
+      term(std::make_pair(T(0), T(1)))
+    {
+    }
+
+    result_type operator()()
+    {
+      const result_type result = term;
+      ++a; ++b; ++c; ++n;
+      numer = -(((a * b) * z) / (c * n));
+      term = std::make_pair(numer, 1 - numer);
+      return result;
+    }
+
+  private:
+    unsigned n;
+    T a, b, c;
+    const T z;
+    T numer;
+    result_type term;
+  };
+
+  template <class T, unsigned p, unsigned q, class Policy>
+  inline T compute_cf_pfq(detail::hypergeometric_pfq_cf_term<T, p, q>& term, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+    const T result = tools::continued_fraction_b(
+      term,
+      boost::math::policies::get_epsilon<T, Policy>(),
+      max_iter);
+    boost::math::policies::check_series_iterations<T>(
+      "boost::math::hypergeometric_pfq_cf<%1%>(%1%,%1%,%1%)",
+      max_iter,
+      pol);
+    return result;
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_0f1_cf(const T& b, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_cf_term<T, 0u, 1u> f(b, z);
+    T result = detail::compute_cf_pfq(f, pol);
+    result = ((z / b) / result) + 1;
+    return result;
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_1f0_cf(const T& a, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_cf_term<T, 1u, 0u> f(a, z);
+    T result = detail::compute_cf_pfq(f, pol);
+    result = ((a * z) / result) + 1;
+    return result;
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_cf(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_cf_term<T, 1u, 1u> f(a, b, z);
+    T result = detail::compute_cf_pfq(f, pol);
+    result = (((a * z) / b) / result) + 1;
+    return result;
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_1f2_cf(const T& a, const T& b, const T& c, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_cf_term<T, 1u, 2u> f(a, b, c, z);
+    T result = detail::compute_cf_pfq(f, pol);
+    result = (((a * z) / (b * c)) / result) + 1;
+    return result;
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_2f1_cf(const T& a, const T& b, const T& c, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_cf_term<T, 2u, 1u> f(a, b, c, z);
+    T result = detail::compute_cf_pfq(f, pol);
+    result = ((((a * b) * z) / c) / result) + 1;
+    return result;
+  }
+
+  } } } // namespaces
+
+#endif // BOOST_MATH_HYPERGEOMETRIC_CF_HPP

include/boost/math/special_functions/detail/hypergeometric_pade.hpp

@@ -0,0 +1,131 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+#ifndef BOOST_MATH_HYPERGEOMETRIC_PADE_HPP
+  #define BOOST_MATH_HYPERGEOMETRIC_PADE_HPP
+
+  namespace boost{ namespace math{ namespace detail{
+
+  // Luke: C ---------- SUBROUTINE R1F1P(CP, Z, A, B, N) ----------
+  // Luke: C ----- PADE APPROXIMATION OF 1F1( 1 ; CP ; -Z ) -------
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_pade(const T& cp, const T& zp, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+
+    static const T one = T(1);
+
+    // Luke: C ------------- INITIALIZATION -------------
+    const T z = -zp;
+    const T zz = z * z;
+    T b0 = one;
+    T a0 = one;
+    T xi1 = one;
+    T ct1 = cp + one;
+    T cp1 = cp - one;
+
+    T b1 = one + (z / ct1);
+    T a1 = b1 - (z / cp);
+
+    const unsigned max_iterations = boost::math::policies::get_max_series_iterations<Policy>();
+
+    T b2 = T(0), a2 = T(0);
+    T result = T(0), prev_result = a1 / b1;
+
+    for (unsigned k = 1; k < max_iterations; ++k)
+    {
+      // Luke: C ----- CALCULATION OF THE MULTIPLIERS -----
+      // Luke: C ----------- FOR THE RECURSION ------------
+      const T ct2 = ct1 * ct1;
+      const T g1 = one + ((cp1 / (ct2 + ct1 + ct1)) * z);
+      const T g2 = ((xi1 / (ct2 - one)) * ((xi1 + cp1) / ct2)) * zz;
+
+      // Luke: C ------- THE RECURRENCE RELATIONS ---------
+      // Luke: C ------------ ARE AS FOLLOWS --------------
+      b2 = (g1 * b1) + (g2 * b0);
+      a2 = (g1 * a1) + (g2 * a0);
+
+      prev_result = result;
+      result = a2 / b2;
+
+      // condition for interruption
+      if ((fabs(result) * boost::math::tools::epsilon<T>()) > fabs(result - prev_result))
+        break;
+
+      b0 = b1; b1 = b2;
+      a0 = a1; a1 = a2;
+
+      ct1 += 2;
+      ++xi1;
+    }
+
+    return a2 / b2;
+  }
+
+  // Luke: C -------- SUBROUTINE R2F1P(BP, CP, Z, A, B, N) --------
+  // Luke: C ---- PADE APPROXIMATION OF 2F1( 1 , BP; CP ; -Z ) ----
+  template <class T, class Policy>
+  inline T hypergeometric_2f1_pade(const T& bp, const T& cp, const T& zp, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+
+    static const T one = T(1);
+
+    // Luke: C ---------- INITIALIZATION -----------
+    const T z = -zp;
+    const T zz = z * z;
+    T b0 = one;
+    T a0 = one;
+    T xi1 = one;
+    T ct1 = cp;
+    const T b1c1 = (cp - one) * (bp - one);
+
+    T b1 = one + ((z / (cp + one)) * (bp + one));
+    T a1 = b1 - ((bp / cp) * z);
+
+    const unsigned max_iterations = boost::math::policies::get_max_series_iterations<Policy>();
+
+    T b2 = T(0), a2 = T(0);
+    T result = T(0), prev_result = a1 / b1;
+
+    for (unsigned k = 1; k < max_iterations; ++k)
+    {
+      // Luke: C ----- CALCULATION OF THE MULTIPLIERS -----
+      // Luke: C ----------- FOR THE RECURSION ------------
+      const T ct2 = ct1 + xi1;
+      const T ct3 = ct2 * ct2;
+      const T g2 = (((((ct1 / ct3) * (bp - ct1)) / (ct3 - one)) * xi1) * (bp + xi1)) * zz;
+      ++xi1;
+      const T g1 = one + (((((xi1 + xi1) * ct1) + b1c1) / (ct3 + ct2 + ct2)) * z);
+
+      // Luke: C ------- THE RECURRENCE RELATIONS ---------
+      // Luke: C ------------ ARE AS FOLLOWS --------------
+      b2 = (g1 * b1) + (g2 * b0);
+      a2 = (g1 * a1) + (g2 * a0);
+
+      prev_result = result;
+      result = a2 / b2;
+
+      // condition for interruption
+      if ((fabs(result) * boost::math::tools::epsilon<T>()) > fabs(result - prev_result))
+        break;
+
+      b0 = b1; b1 = b2;
+      a0 = a1; a1 = a2;
+
+      ++ct1;
+    }
+
+    return a2 / b2;
+  }
+
+  } } } // namespaces
+
+#endif // BOOST_MATH_HYPERGEOMETRIC_PADE_HPP

include/boost/math/special_functions/detail/hypergeometric_rational.hpp

@@ -0,0 +1,165 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+#ifndef BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP
+  #define BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP
+
+  #include <boost/array.hpp>
+
+  namespace boost{ namespace math{ namespace detail{
+
+  // Luke: C ------- SUBROUTINE R1F1P(AP, CP, Z, A, B, N) ---------
+  // Luke: C --- RATIONAL APPROXIMATION OF 1F1( AP ; CP ; -Z ) ----
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_rational(const T& ap, const T& cp, const T& zp, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+
+    static const T zero = T(0), one = T(1), two = T(2), three = T(3);
+
+    // Luke: C ------------- INITIALIZATION -------------
+    const T z = -zp;
+    const T z2 = z / two;
+
+    T ct1 = ap * (z / cp);
+    T ct2 = z2 / (one + cp);
+    T xn3 = zero;
+    T xn2 = one;
+    T xn1 = two;
+    T xn0 = three;
+
+    T b1 = one;
+    T a1 = one;
+    T b2 = one + ((one + ap) * (z2 / cp));
+    T a2 = b2 - ct1;
+    T b3 = one + ((two + b2) * (((two + ap) / three) * ct2));
+    T a3 = b3 - ((one + ct2) * ct1);
+    ct1 = three;
+
+    const unsigned max_iterations = boost::math::policies::get_max_series_iterations<Policy>();
+
+    T a4 = T(0), b4 = T(0);
+    T result = T(0), prev_result = a3 / b3;
+
+    for (unsigned k = 2; k < max_iterations; ++k)
+    {
+      // Luke: C ----- CALCULATION OF THE MULTIPLIERS -----
+      // Luke: C ----------- FOR THE RECURSION ------------
+      ct2 = (z2 / ct1) / (cp + xn1);
+      const T g1 = one + (ct2 * (xn2 - ap));
+      ct2 *= ((ap + xn1) / (cp + xn2));
+      const T g2 = ct2 * ((cp - xn1) + (((ap + xn0) / (ct1 + two)) * z2));
+      const T g3 = ((ct2 * z2) * (((z2 / ct1) / (ct1 - two)) * ((ap + xn2)) / (cp + xn3))) * (ap - xn2);
+
+      // Luke: C ------- THE RECURRENCE RELATIONS ---------
+      // Luke: C ------------ ARE AS FOLLOWS --------------
+      b4 = (g1 * b3) + (g2 * b2) + (g3 * b1);
+      a4 = (g1 * a3) + (g2 * a2) + (g3 * a1);
+
+      prev_result = result;
+      result = a4 / b4;
+
+      // condition for interruption
+      if ((fabs(result) * boost::math::tools::epsilon<T>()) > fabs(result - prev_result))
+        break;
+
+      b1 = b2; b2 = b3; b3 = b4;
+      a1 = a2; a2 = a3; a3 = a4;
+
+      xn3 = xn2; xn2 = xn1; xn1 = xn0; ++xn0;
+      ct1 += two;
+    }
+
+    return result;
+  }
+
+  // Luke: C ----- SUBROUTINE R2F1P(AB, BP, CP, Z, A, B, N) -------
+  // Luke: C -- RATIONAL APPROXIMATION OF 2F1( AB , BP; CP ; -Z ) -
+  template <class T, class Policy>
+  inline T hypergeometric_2f1_rational(const T& ap, const T& bp, const T& cp, const T& zp, const unsigned n, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+
+    static const T one = T(1), two = T(2), three = T(3), four = T(4),
+                   six = T(6), half_7 = T(3.5), half_3 = T(1.5), forth_3 = T(0.75);
+
+    // Luke: C ------------- INITIALIZATION -------------
+    const T z = -zp;
+    const T z2 = z / two;
+
+    T sabz = (ap + bp) * z;
+    const T ab = ap * bp;
+    const T abz = ab * z;
+    const T abz1 = z + (abz + sabz);
+    const T abz2 = abz1 + (sabz + (three * z));
+    const T cp1 = cp + one;
+    const T ct1 = cp1 + cp1;
+
+    T b1 = one;
+    T a1 = one;
+    T b2 = one + (abz1 / (cp + cp));
+    T a2 = b2 - (abz / cp);
+    T b3 = one + ((abz2 / ct1) * (one + (abz1 / ((-six) + (three * ct1)))));
+    T a3 = b3 - ((abz / cp) * (one + ((abz2 - abz1) / ct1)));
+    sabz /= four;
+
+    const T abz1_div_4 = abz1 / four;
+    const T cp1_inc = cp1 + one;
+    const T cp1_mul_cp1_inc = cp1 * cp1_inc;
+
+    boost::array<T, 9u> d = {{
+      ((half_7 - ab) * z2) - sabz,
+      abz1_div_4,
+      abz1_div_4 - (two * sabz),
+      cp1_inc,
+      cp1_mul_cp1_inc,
+      cp * cp1_mul_cp1_inc,
+      half_3,
+      forth_3,
+      forth_3 * z
+    }};
+
+    T xi = three;
+    T a4 = T(0), b4 = T(0);
+    for (unsigned k = 2; k < n; ++k)
+    {
+      // Luke: C ----- CALCULATION OF THE MULTIPLIERS -----
+      // Luke: C ----------- FOR THE RECURSION ------------
+      T g3 = (d[2] / d[7]) * (d[1] / d[5]);
+      d[1] += d[8] + sabz;
+      d[2] += d[8] - sabz;
+      g3 *= d[1] / d[6];
+      T g1 = one + (((d[1] + d[0]) / d[6]) / d[3]);
+      T g2 = (d[1] / d[4]) / d[6];
+      d[7] += two * d[6];
+      ++d[6];
+      g2 *= cp1 - (xi + ((d[2] + d[0]) / d[6]));
+
+      // Luke: C ------- THE RECURRENCE RELATIONS ---------
+      // Luke: C ------------ ARE AS FOLLOWS --------------
+      b4 = (g1 * b3) + (g2 * b2) + (g3 * b1);
+      a4 = (g1 * a3) + (g2 * a2) + (g3 * a1);
+      b1 = b2; b2 = b3; b3 = b4;
+      a1 = a2; a2 = a3; a3 = a4;
+
+      d[8] += z2;
+      d[0] += two * d[8];
+      d[5] += three * d[4];
+      d[4] += two * d[3];
+      ++d[3];
+      ++xi;
+    }
+
+    return a4 / b4;
+  }
+
+  } } } // namespaces
+
+#endif // BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP

include/boost/math/special_functions/detail/hypergeometric_separated_series.hpp

@@ -0,0 +1,51 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+#ifndef BOOST_MATH_HYPERGEOMETRIC_SEPARATED_SERIES_HPP
+  #define BOOST_MATH_HYPERGEOMETRIC_SEPARATED_SERIES_HPP
+
+  namespace boost { namespace math { namespace detail {
+
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_separated_series(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+
+    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+    const T factor = policies::get_epsilon<T, Policy>();
+
+    T denom = 1, numer = 1;
+    T intermediate_result = 1, result = 1;
+    T a_pochhammer = a, z_pow = z;
+    unsigned N = 0;
+    while (--max_iter)
+    {
+      ++N;
+      const T mult = (((b + N) - 1) * N);
+      denom *= mult; numer *= mult;
+      numer += a_pochhammer * z_pow;
+
+      result = numer / denom;
+
+      if (fabs(factor * result) > fabs(result - intermediate_result))
+        break;
+
+      intermediate_result = result;
+
+      a_pochhammer *= (a + N);
+      z_pow *= z;
+    }
+
+    return result;
+  }
+
+  } } } // namespaces
+
+#endif // BOOST_MATH_HYPERGEOMETRIC_SEPARATED_SERIES_HPP

include/boost/math/special_functions/detail/hypergeometric_series.hpp

@@ -0,0 +1,242 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+//
+#ifndef BOOST_MATH_HYPERGEOMETRIC_SERIES_HPP
+  #define BOOST_MATH_HYPERGEOMETRIC_SERIES_HPP
+
+#include <boost/math/tools/series.hpp>
+#include <boost/math/policies/error_handling.hpp>
+
+  namespace boost { namespace math { namespace detail {
+
+  // primary template for term of Taylor series
+  template <class T, unsigned p, unsigned q>
+  struct hypergeometric_pfq_generic_series_term;
+
+  // partial specialization for 0F1
+  template <class T>
+  struct hypergeometric_pfq_generic_series_term<T, 0u, 1u>
+  {
+    typedef T result_type;
+
+    hypergeometric_pfq_generic_series_term(const T& b, const T& z)
+       : n(0), term(1), b(b), z(z)
+    {
+    }
+
+    T operator()()
+    {
+      BOOST_MATH_STD_USING
+      const T r = term;
+      term *= ((1 / ((b + n) * (n + 1))) * z);
+      ++n;
+      return r;
+    }
+
+  private:
+    unsigned n;
+    T term;
+    const T b, z;
+  };
+
+  // partial specialization for 1F0
+  template <class T>
+  struct hypergeometric_pfq_generic_series_term<T, 1u, 0u>
+  {
+    typedef T result_type;
+
+    hypergeometric_pfq_generic_series_term(const T& a, const T& z)
+       : n(0), term(1), a(a), z(z)
+    {
+    }
+
+    T operator()()
+    {
+      BOOST_MATH_STD_USING
+      const T r = term;
+      term *= (((a + n) / (n + 1)) * z);
+      ++n;
+      return r;
+    }
+
+  private:
+    unsigned n;
+    T term;
+    const T a, z;
+  };
+
+  // partial specialization for 1F1
+  template <class T>
+  struct hypergeometric_pfq_generic_series_term<T, 1u, 1u>
+  {
+    typedef T result_type;
+
+    hypergeometric_pfq_generic_series_term(const T& a, const T& b, const T& z)
+       : n(0), term(1), a(a), b(b), z(z)
+    {
+    }
+
+    T operator()()
+    {
+      BOOST_MATH_STD_USING
+      const T r = term;
+      term *= (((a + n) / ((b + n) * (n + 1))) * z);
+      ++n;
+      return r;
+    }
+
+  private:
+    unsigned n;
+    T term;
+    const T a, b, z;
+  };
+
+  // partial specialization for 1F2
+  template <class T>
+  struct hypergeometric_pfq_generic_series_term<T, 1u, 2u>
+  {
+    typedef T result_type;
+
+    hypergeometric_pfq_generic_series_term(const T& a, const T& b1, const T& b2, const T& z)
+       : n(0), term(1), a(a), b1(b1), b2(b2), z(z)
+    {
+    }
+
+    T operator()()
+    {
+      BOOST_MATH_STD_USING
+      const T r = term;
+      term *= (((a + n) / ((b1 + n) * (b2 + n) * (n + 1))) * z);
+      ++n;
+      return r;
+    }
+
+  private:
+    unsigned n;
+    T term;
+    const T a, b1, b2, z;
+  };
+
+  // partial specialization for 2F0
+  template <class T>
+  struct hypergeometric_pfq_generic_series_term<T, 2u, 0u>
+  {
+    typedef T result_type;
+
+    hypergeometric_pfq_generic_series_term(const T& a1, const T& a2, const T& z)
+       : n(0), term(1), a1(a1), a2(a2), z(z)
+    {
+    }
+
+    T operator()()
+    {
+      BOOST_MATH_STD_USING
+      const T r = term;
+      term *= (((a1 + n) * (a2 + n) / (n + 1)) * z);
+      ++n;
+      return r;
+    }
+
+  private:
+    unsigned n;
+    T term;
+    const T a1, a2, z;
+  };
+
+  // partial specialization for 2F1
+  template <class T>
+  struct hypergeometric_pfq_generic_series_term<T, 2u, 1u>
+  {
+    typedef T result_type;
+
+    hypergeometric_pfq_generic_series_term(const T& a1, const T& a2, const T& b, const T& z)
+       : n(0), term(1), a1(a1), a2(a2), b(b), z(z)
+    {
+    }
+
+    T operator()()
+    {
+      BOOST_MATH_STD_USING
+      const T r = term;
+      term *= (((a1 + n) * (a2 + n) / ((b + n) * (n + 1))) * z);
+      ++n;
+      return r;
+    }
+
+  private:
+    unsigned n;
+    T term;
+    const T a1, a2, b, z;
+  };
+
+  // we don't need to define extra check and make a polinom from
+  // series, when p(i) and q(i) are negative integers and p(i) >= q(i)
+  // as described in functions.wolfram.alpha, because we always
+  // stop summation when result (in this case numerator) is zero.
+  template <class T, unsigned p, unsigned q, class Policy>
+  inline T sum_pfq_series(detail::hypergeometric_pfq_generic_series_term<T, p, q>& term, const Policy& pol)
+  {
+    BOOST_MATH_STD_USING
+    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+    const T zero = 0;
+    const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
+#else
+    const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
+#endif
+    policies::check_series_iterations<T>("boost::math::hypergeometric_pfq_generic_series<%1%>(%1%,%1%,%1%)", max_iter, pol);
+    return result;
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_0f1_generic_series(const T& b, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_generic_series_term<T, 0u, 1u> s(b, z);
+    return detail::sum_pfq_series(s, pol);
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_1f0_generic_series(const T& a, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_generic_series_term<T, 1u, 0u> s(a, z);
+    return detail::sum_pfq_series(s, pol);
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_1f1_generic_series(const T& a, const T& b, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_generic_series_term<T, 1u, 1u> s(a, b, z);
+    return detail::sum_pfq_series(s, pol);
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_1f2_generic_series(const T& a, const T& b1, const T& b2, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_generic_series_term<T, 1u, 2u> s(a, b1, b2, z);
+    return detail::sum_pfq_series(s, pol);
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_2f0_generic_series(const T& a1, const T& a2, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_generic_series_term<T, 2u, 0u> s(a1, a2, z);
+    return detail::sum_pfq_series(s, pol);
+  }
+
+  template <class T, class Policy>
+  inline T hypergeometric_2f1_generic_series(const T& a1, const T& a2, const T& b, const T& z, const Policy& pol)
+  {
+    detail::hypergeometric_pfq_generic_series_term<T, 2u, 1u> s(a1, a2, b, z);
+    return detail::sum_pfq_series(s, pol);
+  }
+
+  } } } // namespaces
+
+#endif // BOOST_MATH_HYPERGEOMETRIC_SERIES_HPP

include/boost/math/special_functions/hypergeometric_0F1.hpp

@@ -0,0 +1,110 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef _BOOST_HYPERGEOMETRIC_0F1_HPP_
+  #define _BOOST_HYPERGEOMETRIC_0F1_HPP_
+
+#include <boost/math/policies/policy.hpp>
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/special_functions/detail/hypergeometric_series.hpp>
+#include <boost/math/special_functions/detail/hypergeometric_0f1_bessel.hpp>
+
+namespace boost { namespace math { namespace detail {
+
+
+   template <class T>
+   struct hypergeometric_0F1_cf
+   {
+      T b, z;
+      int k;
+      hypergeometric_0F1_cf(T b_, T z_) : b(b_), z(z_), k(0) {}
+      typedef std::pair<T, T> result_type;
+
+      result_type operator()()
+      {
+         ++k;
+         return std::make_pair(-z / ((k + 1) * (b + k)), 1 + z / ((k + 1) * (b + k)));
+      }
+   };
+
+   template <class T, class Policy>
+   T hypergeometric_0F1_cf_imp(T b, T z, const Policy& pol, const char* function)
+   {
+      hypergeometric_0F1_cf<T> evaluator(b, z);
+      boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+      T cf = tools::continued_fraction_a(evaluator, policies::get_epsilon<T, Policy>(), max_iter);
+      policies::check_series_iterations<T>(function, max_iter, pol);
+      return 1 + z / (b * (1 + cf));
+   }
+
+
+   template <class T, class Policy>
+   inline T hypergeometric_0f1_imp(const T& b, const T& z, const Policy& pol)
+   {
+      const char* function = "boost::math::hypergeometric_0f1<%1%,%1%>(%1%, %1%)";
+      BOOST_MATH_STD_USING
+
+         // some special cases
+         if (z == 0)
+            return T(1);
+
+      if ((b <= 0) && (b == floor(b)))
+         return policies::raise_pole_error<T>(
+            function,
+            "Evaluation of 0f1 with nonpositive integer b = %1%.", b, pol);
+
+      if (z < -5 && b > -5)
+      {
+         // Series is alternating and divergent, need to do something else here,
+         // Bessel function relation is much more accurate, unless |b| is similarly
+         // large to |z|, otherwise the CF formula suffers from cancellation when 
+         // the result would be very small.
+         if (fabs(z / b) > 4)
+            return hypergeometric_0f1_bessel(b, z, pol);
+         return hypergeometric_0F1_cf_imp(b, z, pol, function);
+      }
+      // evaluation through Taylor series looks
+      // more precisious than Bessel relation:
+      // detail::hypergeometric_0f1_bessel(b, z, pol);
+      return detail::hypergeometric_0f1_generic_series(b, z, pol);
+   }
+
+} // namespace detail
+
+template <class T1, class T2, class Policy>
+inline typename tools::promote_args<T1, T2>::type hypergeometric_0F1(T1 b, T2 z, const Policy& /* pol */)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+      typedef typename tools::promote_args<T1, T2>::type result_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename policies::normalise<
+      Policy,
+      policies::promote_float<false>,
+      policies::promote_double<false>,
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+   return policies::checked_narrowing_cast<result_type, Policy>(
+      detail::hypergeometric_0f1_imp<value_type>(
+         static_cast<value_type>(b),
+         static_cast<value_type>(z),
+         forwarding_policy()),
+      "boost::math::hypergeometric_0f1<%1%>(%1%,%1%)");
+}
+
+template <class T1, class T2>
+inline typename tools::promote_args<T1, T2>::type hypergeometric_0F1(T1 b, T2 z)
+{
+   return hypergeometric_0F1(b, z, policies::policy<>());
+}
+
+
+} } // namespace boost::math
+
+#endif // _BOOST_HYPERGEOMETRIC_2014_04_07_HPP_

include/boost/math/special_functions/hypergeometric_1F0.hpp

@@ -0,0 +1,73 @@
+
+///////////////////////////////////////////////////////////////////////////////
+//  Copyright 2014 Anton Bikineev
+//  Copyright 2014 Christopher Kormanyos
+//  Copyright 2014 John Maddock
+//  Copyright 2014 Paul Bristow
+//  Distributed under the Boost
+//  Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef _BOOST_HYPERGEOMETRIC_1F0_HPP_
+  #define _BOOST_HYPERGEOMETRIC_1F0_HPP_
+
+#include <boost/math/policies/policy.hpp>
+#include <boost/math/policies/error_handling.hpp>
+
+
+namespace boost { namespace math { namespace detail {
+
+template <class T, class Policy>
+inline T hypergeometric_1f0_imp(const T& a, const T& z, const Policy& pol)
+{
+   static const char* function = "boost::math::hypergeometric_1f0<%1%,%1%>(%1%, %1%)";
+   BOOST_MATH_STD_USING // pow
+
+   if (z == 1)
+      return policies::raise_pole_error<T>(
+         function,
+         "Evaluation of 1F0 with z = %1%.",
+         z,
+         pol);
+   if (1 - z < 0)
+   {
+      if (floor(a) != a)
+         return policies::raise_domain_error<T>(function,
+            "Result is complex when a is non-integral and z > 1, but got z = %1%", z, pol);
+   }
+   // more naive and convergent method than series
+   return pow(T(1 - z), T(-a));
+}
+
+} // namespace detail
+
+template <class T1, class T2, class Policy>
+inline typename tools::promote_args<T1, T2>::type hypergeometric_1f0(T1 a, T2 z, const Policy&)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename tools::promote_args<T1, T2>::type result_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename policies::normalise<
+      Policy,
+      policies::promote_float<false>,
+      policies::promote_double<false>,
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+   return policies::checked_narrowing_cast<result_type, Policy>(
+      detail::hypergeometric_1f0_imp<value_type>(
+         static_cast<value_type>(a),
+         static_cast<value_type>(z),
+         forwarding_policy()),
+      "boost::math::hypergeometric_1f0<%1%>(%1%,%1%)");
+}
+
+template <class T1, class T2>
+inline typename tools::promote_args<T1, T2>::type hypergeometric_1f0(T1 a, T2 z)
+{
+   return hypergeometric_1f0(a, z, policies::policy<>());
+}
+
+
+  } } // namespace boost::math
+
+#endif // _BOOST_HYPERGEOMETRIC_2014_04_07_HPP_