Skip to content

Adding logarithm of incomplete gamma function #1338

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
JacobHass8 wants to merge 4 commits into from
+190 −10
Closed

Adding logarithm of incomplete gamma function #1338

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
da8b635
Added some context to continued_fractions example
JacobHass8 Dec 11, 2025
0a2b01b
Added log arithmetic to continued_fraction_b
JacobHass8 Dec 12, 2025
c80492f
Added log-arithmetic for continued_fraction_a
JacobHass8 Dec 12, 2025
dfb598f
Removed previous commit issues. Addded function to gamma.hpp
JacobHass8 Dec 14, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 7 additions & 9 deletions example/continued_fractions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
33 changes: 33 additions & 0 deletions include/boost/math/special_functions/gamma.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
{
Expand All @@ -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
{
Expand Down
29 changes: 28 additions & 1 deletion include/boost/math/special_functions/logaddexp.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Copy link
Collaborator

@jzmaddock jzmaddock Dec 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The trick here, is to add using std::real above, and then call real(temp), that way ADL can find boost::multiprecision::real in the multiprecision case (likewise for any other type that has their own namespace and is forbidden from putting stuff in std::).

{
return x1 + log(1.0 + exp(-temp));
Copy link
Collaborator

@jzmaddock jzmaddock Dec 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a really bad change that will kill precision when temp is small.

Copy link
Contributor Author

@JacobHass8 JacobHass8 Dec 14, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've reverted this for the logaddexp function and created a new function for adding complex numbers. I haven't been able to find an equivalent of log1p that works with complex numbers. Could I just make a new function clog1p which computes the Taylor series of $log(1+x)$ for $|x| &lt; 1$? This seems easy enough that I'm surprised I haven't found an implementation.

}

return x2 + log(1.0 + exp(temp));
}
}
}// Namespace boost::math
122 changes: 122 additions & 0 deletions include/boost/math/tools/fraction.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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{

Expand Down Expand Up @@ -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>
Expand Down Expand Up @@ -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 {

//
Expand Down Expand Up @@ -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>
Expand Down Expand Up @@ -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
Expand Down