Fix ibeta for large arguments

example/continued_fractions.cpp

@@ -122,29 +122,27 @@ inline boost::multiprecision::cpp_complex_50 gamma_Q_as_fraction(const boost::mu
    return pow(z, a) / (exp(z) * (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps)));
 }
 
-
 int main()
 {
    using namespace boost::math::tools;
 
    //[cf_gr
    golden_ratio_fraction<double> func;
-   double gr = continued_fraction_a(
+   double gr = continued_fraction_b(
       func,
-      std::numeric_limits<double>::epsilon());
+      std::numeric_limits<double>::epsilon()
+      );
    std::cout << "The golden ratio is: " << gr << std::endl;
    //]
 
-   std::cout << tan(0.5) << std::endl;
+   std::cout << "tan(0.5)=" << tan(0.5) << std::endl;
 
    std::complex<double> arg(3, 2);
-   std::cout << expint_as_fraction(5, arg) << std::endl;
+   std::cout << "E_5(3+2i)=" << expint_as_fraction(5, arg) << std::endl;
 
    std::complex<double> a(3, 3), z(3, 2);
-   std::cout << gamma_Q_as_fraction(a, z) << std::endl;
+   std::cout << "Gamma(3+3i, 3+2i)=" << gamma_Q_as_fraction(a, z) << std::endl;
 
    boost::multiprecision::cpp_complex_50 am(3, 3), zm(3, 2);
-   std::cout << gamma_Q_as_fraction(am, zm) << std::endl;
-
-   return 0;
+   std::cout << "Gamma(3+3i, 3+2i)=" << gamma_Q_as_fraction(am, zm) << std::endl;
 }

include/boost/math/special_functions/beta.hpp

@@ -1142,7 +1142,20 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no
 
    T x0 = a / (a + b);
    T y0 = b / (a + b);
-   T nu = x0 * log(x / x0) + y0 * log(y / y0);
+
+   // Expand nu about x0
+   T nu = 0;
+   for (int i=2; i<5; i++)
+   {
+      nu += pow(x-x0, i) / i * (pow(x0, -(i-1)) - pow(x0-1, -(i-1))) * pow(-1, i+1);
+   }
+   // Calculate the next term in the series
+   T remainder = pow(x-x0, 5) / 5 * (pow(x0, -4) - pow(x0-1, 4));
+
+   // If the remainder is large, then fall back to using the log formula
+   if (remainder >= tools::forth_root_epsilon<T>()){
+      nu = x0 * log(x / x0) + y0 * log(y / y0);
+   }
    //
    // Above compution is unstable, force nu to zero if
    // something went wrong:
@@ -1181,7 +1194,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no
    T mul = 1;
    if (!normalised)
       mul = boost::math::beta(a, b, pol);
-
+   T log_erf_remainder = -0.5 * log(2 * constants::pi<T>() * (a+b)) + a * log(x / x0) + b * log((1-x) / (1-x0)) + log(abs(1/nu - sqrt(x0 * (1-x0)) / (x-x0)));
    return mul * ((invert ? (1 + boost::math::erf(-nu * sqrt((a + b) / 2), pol)) / 2 : boost::math::erfc(-nu * sqrt((a + b) / 2), pol) / 2));
 }
 

include/boost/math/special_functions/gamma.hpp

@@ -35,6 +35,7 @@
 #include <boost/math/special_functions/detail/igamma_large.hpp>
 #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
 #include <boost/math/special_functions/detail/lgamma_small.hpp>
+#include <boost/math/tools/fraction.hpp>
 
 // Only needed for types larger than double
 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
@@ -393,6 +394,29 @@ struct upper_incomplete_gamma_fract
    }
 };
 
+template <class T>
+struct complex_upper_incomplete_gamma_fract
+{
+private:
+   typedef typename T::value_type scalar_type;
+   T z, a;
+   int k;
+public:
+   typedef std::pair<T, T> result_type;
+
+   complex_upper_incomplete_gamma_fract(T a1, T z1)
+      : z(z1 - a1 + scalar_type(1)), a(a1), k(0)
+   {
+   }
+
+   result_type operator()()
+   {
+      ++k;
+      z += scalar_type(2);
+      return result_type(scalar_type(k) * (a - scalar_type(k)), z);
+   }
+};
+
 template <class T>
 BOOST_MATH_GPU_ENABLED inline T upper_gamma_fraction(T a, T z, T eps)
 {
@@ -403,6 +427,15 @@ BOOST_MATH_GPU_ENABLED inline T upper_gamma_fraction(T a, T z, T eps)
    return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
 }
 
+template <class T>
+inline std::complex<T> log_upper_gamma_fraction(const std::complex<T>& a, const std::complex<T>& z)
+{
+   complex_upper_incomplete_gamma_fract<std::complex<T> > f(a, z);
+   std::complex<T> eps(std::numeric_limits<T>::epsilon());
+   boost::math::uintmax_t max_iter = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
+   return a * log(z) - z - boost::math::logaddexpcomplex(log(z - a + T(1)), boost::math::tools::log_continued_fraction_a(f, eps, max_iter));
+}
+
 template <class T>
 struct lower_incomplete_gamma_series
 {

include/boost/math/special_functions/logaddexp.hpp

@@ -38,4 +38,31 @@ Real logaddexp(Real x1, Real x2) noexcept
     return x2 + log1p(exp(temp));
 }
 
-}} // Namespace boost::math
+template <typename T>
+T logaddexpcomplex(T x1, T x2) noexcept
+{
+    using std::log;
+    using std::exp;
+    using std::abs;
+
+    // Validate inputs first
+    if (!(boost::math::isfinite)(x1))
+    {
+        return x1;
+    }
+    else if (!(boost::math::isfinite)(x2))
+    {
+        return x2;
+    }
+
+    const T temp = x1 - x2;
+
+    if (std::real(temp) > 0)
+    {
+        return x1 + log(1.0 + exp(-temp));
+    }
+
+    return x2 + log(1.0 + exp(temp));
+}
+}
+}// Namespace boost::math

include/boost/math/tools/fraction.hpp

@@ -18,6 +18,9 @@
 #include <boost/math/tools/precision.hpp>
 #include <boost/math/tools/complex.hpp>
 #include <boost/math/tools/cstdint.hpp>
+#include <boost/math/special_functions/logaddexp.hpp>
+#include <limits>
+#include <iostream>
 
 namespace boost{ namespace math{ namespace tools{
 
@@ -159,6 +162,55 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type
    return f;
 }
 
+template <typename Gen, typename U>
+BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type log_continued_fraction_b_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
+      noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 
+      #ifndef BOOST_MATH_HAS_GPU_SUPPORT
+      // SYCL can not handle this condition so we only check float on that platform
+      && noexcept(std::declval<Gen>()())
+      #endif
+      )
+{
+   BOOST_MATH_STD_USING // ADL of std names
+
+   using traits = detail::fraction_traits<Gen>;
+   using result_type = typename traits::result_type;
+   using value_type = typename traits::value_type;
+   using integer_type = typename integer_scalar_type<result_type>::type;
+   using scalar_type = typename scalar_type<result_type>::type;
+
+   integer_type const zero(0), one(1);
+
+   result_type tiny = detail::tiny_value<result_type>::get();
+   scalar_type terminator = abs(factor);
+
+   value_type v = g();
+
+   result_type f, C, D, delta;
+   f = log(traits::b(v));
+   if(f == -std::numeric_limits<result_type>::infinity())
+      f = log(tiny);
+   C = f;
+   D = log(tiny);
+
+   boost::math::uintmax_t counter(max_terms);
+   do{
+      v = g();
+      D = -logaddexp(log(traits::b(v)), log(traits::a(v)) + D);
+      if(D == -std::numeric_limits<result_type>::infinity())
+         D = log(tiny);
+      C = logaddexp(log(traits::b(v)), log(traits::a(v)) - C);
+      if(C == -std::numeric_limits<result_type>::infinity())
+         C = log(tiny);
+      delta = C + D;
+      f = f + delta;
+   }while((abs(delta) > terminator) && --counter);
+
+   max_terms = max_terms - counter;
+
+   return f;
+}
+
 } // namespace detail
 
 template <typename Gen, typename U>
@@ -219,6 +271,18 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type
    return detail::continued_fraction_b_impl(g, factor, max_terms);
 }
 
+template <typename Gen, typename U>
+BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type log_continued_fraction_b(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
+   noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 
+         #ifndef BOOST_MATH_HAS_GPU_SUPPORT
+         && noexcept(std::declval<Gen>()())
+         #endif
+         )
+{
+   return detail::log_continued_fraction_b_impl(g, factor, max_terms);
+}
+
+
 namespace detail {
 
 //
@@ -286,6 +350,53 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type
    return a0/f;
 }
 
+template <typename Gen, typename U>
+BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type log_continued_fraction_a_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
+   noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 
+   #ifndef BOOST_MATH_HAS_GPU_SUPPORT
+   && noexcept(std::declval<Gen>()())
+   #endif
+   )
+{
+   BOOST_MATH_STD_USING // ADL of std names
+
+   using traits = detail::fraction_traits<Gen>;
+   using result_type = typename traits::result_type;
+   using value_type = typename traits::value_type;
+   using integer_type = typename integer_scalar_type<result_type>::type;
+   using scalar_type = typename scalar_type<result_type>::type;
+
+   integer_type const zero(0), one(1);
+
+   result_type tiny = detail::tiny_value<result_type>::get();
+   scalar_type terminator = abs(factor);
+
+   value_type v = g();
+
+   result_type f, C, D, delta, a0;
+   f = log(traits::b(v));
+   a0 = log(traits::a(v));
+   if(f == -std::numeric_limits<result_type>::infinity())
+      f = log(tiny);
+   C = f;
+   D = log(tiny);
+
+   boost::math::uintmax_t counter(max_terms);
+   do{
+      v = g();
+      D = -logaddexpcomplex(log(traits::b(v)), log(traits::a(v)) + D);
+      if(D == -std::numeric_limits<result_type>::infinity())
+         D = log(tiny);
+      C = logaddexpcomplex(log(traits::b(v)), log(traits::a(v)) - C);
+      if(C == -std::numeric_limits<result_type>::infinity())
+         C = log(tiny);
+      delta = C + D;
+      f = f + delta;
+   }while((abs(delta) > terminator) && --counter);
+
+   return a0 - f;
+}
+
 } // namespace detail
 
 template <typename Gen, typename U>
@@ -347,6 +458,17 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type
    return detail::continued_fraction_a_impl(g, factor, max_terms);
 }
 
+template <typename Gen, typename U>
+BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type log_continued_fraction_a(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
+   noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type) 
+   #ifndef BOOST_MATH_HAS_GPU_SUPPORT
+   && noexcept(std::declval<Gen>()())
+   #endif
+   )
+{
+   return detail::log_continued_fraction_a_impl(g, factor, max_terms);
+}
+
 } // namespace tools
 } // namespace math
 } // namespace boost