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
Closed

Adding logarithm of incomplete gamma function #1338

Show file tree
Hide file tree
Changes from 3 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
41 changes: 28 additions & 13 deletions example/continued_fractions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -115,36 +115,51 @@ inline std::complex<T> gamma_Q_as_fraction(const std::complex<T>& a, const std::
return pow(z, a) / (exp(z) *(z - a + T(1) + boost::math::tools::continued_fraction_a(f, eps)));
}
//]
inline boost::multiprecision::cpp_complex_50 gamma_Q_as_fraction(const boost::multiprecision::cpp_complex_50& a, const boost::multiprecision::cpp_complex_50& z)
// inline boost::multiprecision::cpp_complex_50 gamma_Q_as_fraction(const boost::multiprecision::cpp_complex_50& a, const boost::multiprecision::cpp_complex_50& z)
// {
// upper_incomplete_gamma_fract<boost::multiprecision::cpp_complex_50> f(a, z);
// boost::multiprecision::cpp_complex_50 eps(std::numeric_limits<boost::multiprecision::cpp_complex_50::value_type>::epsilon());
// return pow(z, a) / (exp(z) * (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps)));
// }
template <class T>
inline std::complex<T> log_gamma_Q_as_fraction(const std::complex<T>& a, const std::complex<T>& z)
{
upper_incomplete_gamma_fract<boost::multiprecision::cpp_complex_50> f(a, z);
boost::multiprecision::cpp_complex_50 eps(std::numeric_limits<boost::multiprecision::cpp_complex_50::value_type>::epsilon());
return pow(z, a) / (exp(z) * (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps)));
upper_incomplete_gamma_fract<std::complex<T> > f(a, z);
std::complex<T> eps(std::numeric_limits<T>::epsilon());
return a * log(z) - z - boost::math::logaddexp(log(z - a + T(1)), boost::math::tools::continued_fraction_a(f, eps, true));
}


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::cout << "The golden ratio is: " << gr << std::endl;
std::numeric_limits<double>::epsilon(),
true);
std::cout << "The golden ratio is: " << exp(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;

// Off in the imaginary part by a factor of 2*pi
std::cout << "Without log-arithmetic: ln(Gamma(3+3i, 3+2i))=" << log(gamma_Q_as_fraction(a, z)) << std::endl;
std::cout << "With log-arithmetic: ln(Gamma(3+3i, 3+2i))=" << log_gamma_Q_as_fraction(a, z) << std::endl;

std::complex<double> s(10, 0);
std::complex<double> x(1000, 0);
std::cout << "Gamma(10, 1000)=" << log_gamma_Q_as_fraction(s, x) << std::endl;

boost::multiprecision::cpp_complex_50 am(3, 3), zm(3, 2);
std::cout << gamma_Q_as_fraction(am, zm) << std::endl;
// boost::multiprecision::cpp_complex_50 am(3, 3), zm(3, 2);
// std::cout << "Gamma(3+3i, 3+2i)=" << gamma_Q_as_fraction(am, zm) << std::endl;

return 0;
}
9 changes: 5 additions & 4 deletions include/boost/math/special_functions/logaddexp.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,14 +7,15 @@
#include <limits>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/constants/constants.hpp>
#include <complex>

namespace boost { namespace math {

// Calculates log(exp(x1) + exp(x2))
template <typename Real>
Real logaddexp(Real x1, Real x2) noexcept
{
using std::log1p;
using std::log;
using std::exp;
using std::abs;

Expand All @@ -30,12 +31,12 @@ Real logaddexp(Real x1, Real x2) noexcept

const Real temp = x1 - x2;

if (temp > 0)
if (std::real(temp) > 0)

@jzmaddock jzmaddock Dec 13, 2025

Copy link
Copy Markdown
Collaborator

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 + log1p(exp(-temp));
return x1 + log(1.0 + exp(-temp));

@jzmaddock jzmaddock Dec 13, 2025

Copy link
Copy Markdown
Collaborator

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.

@JacobHass8 JacobHass8 Dec 14, 2025

Copy link
Copy Markdown
Contributor Author

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 + log1p(exp(temp));
return x2 + log(1.0 + exp(temp));
}

}} // Namespace boost::math
136 changes: 92 additions & 44 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,8 @@
#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>

namespace boost{ namespace math{ namespace tools{

Expand Down Expand Up @@ -110,7 +112,7 @@ namespace detail {
//

template <typename Gen, typename U>
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms, bool logArithmetic)

@jzmaddock jzmaddock Dec 13, 2025

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Since the log and non-log versions basically have completely different code, I see this as two different functions?

@JacobHass8 JacobHass8 Dec 14, 2025

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

I've separated this into a new function in gamma.hpp. I wasn't sure what the correct implementation should have been.

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
Expand Down Expand Up @@ -140,19 +142,40 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type
C = f;
D = 0;

if (logArithmetic){
f = log(f);
C = log(C);
D = log(tiny);
}

boost::math::uintmax_t counter(max_terms);
do{
v = g();
D = traits::b(v) + traits::a(v) * D;
if(D == result_type(0))
D = tiny;
C = traits::b(v) + traits::a(v) / C;
if(C == zero)
C = tiny;
D = one/D;
delta = C*D;
f = f * delta;
}while((abs(delta - one) > terminator) && --counter);
if(!logArithmetic){
do{
v = g();
D = traits::b(v) + traits::a(v) * D;
if(D == result_type(0))
D = tiny;
C = traits::b(v) + traits::a(v) / C;
if(C == zero)
C = tiny;
D = one/D;
delta = C*D;
f = f * delta;
}while((abs(delta - one) > terminator) && --counter);
}
else{
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;

Expand All @@ -162,30 +185,30 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type
} // namespace detail

template <typename Gen, typename U>
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, boost::math::uintmax_t& max_terms, bool logArithmetic=false)
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::continued_fraction_b_impl(g, factor, max_terms);
return detail::continued_fraction_b_impl(g, factor, max_terms, logArithmetic);
}

template <typename Gen, typename U>
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor)
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, const U& factor, bool logArithmetic=false)
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::uintmax_t max_terms = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
return detail::continued_fraction_b_impl(g, factor, max_terms);
return detail::continued_fraction_b_impl(g, factor, max_terms, logArithmetic);
}

template <typename Gen>
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, bool logArithmetic=false)
noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
#ifndef BOOST_MATH_HAS_GPU_SUPPORT
&& noexcept(std::declval<Gen>()())
Expand All @@ -199,11 +222,11 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type

result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
boost::math::uintmax_t max_terms = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
return detail::continued_fraction_b_impl(g, factor, max_terms);
return detail::continued_fraction_b_impl(g, factor, max_terms, logArithmetic);
}

template <typename Gen>
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::math::uintmax_t& max_terms)
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::math::uintmax_t& max_terms, bool logArithmetic=false)
noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
#ifndef BOOST_MATH_HAS_GPU_SUPPORT
&& noexcept(std::declval<Gen>()())
Expand All @@ -216,7 +239,7 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type
using result_type = typename traits::result_type;

result_type factor = ldexp(1.0f, 1 - bits); // 1 / pow(result_type(2), bits);
return detail::continued_fraction_b_impl(g, factor, max_terms);
return detail::continued_fraction_b_impl(g, factor, max_terms, logArithmetic);
}

namespace detail {
Expand All @@ -236,7 +259,7 @@ namespace detail {
// Note that the first a1 and b1 returned by generator Gen are both used.
//
template <typename Gen, typename U>
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a_impl(Gen& g, const U& factor, boost::math::uintmax_t& max_terms, bool logArithmetic=false)

@jzmaddock jzmaddock Dec 13, 2025

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Again this looks like two different functions to me.

noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
#ifndef BOOST_MATH_HAS_GPU_SUPPORT
&& noexcept(std::declval<Gen>()())
Expand Down Expand Up @@ -266,53 +289,78 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type
C = f;
D = 0;

if (logArithmetic){
f = log(f);
C = log(C);
D = log(tiny);
}

boost::math::uintmax_t counter(max_terms);

do{
v = g();
D = traits::b(v) + traits::a(v) * D;
if(D == zero)
D = tiny;
C = traits::b(v) + traits::a(v) / C;
if(C == zero)
C = tiny;
D = one/D;
delta = C*D;
f = f * delta;
}while((abs(delta - one) > terminator) && --counter);
if (!logArithmetic){
do{
v = g();
D = traits::b(v) + traits::a(v) * D;
if(D == zero)
D = tiny;
C = traits::b(v) + traits::a(v) / C;
if(C == zero)
C = tiny;
D = one/D;
delta = C*D;
f = f * delta;
}while((abs(delta - one) > terminator) && --counter);
}
else{
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 a0/f;
if (!logArithmetic){
return a0/f;
}
else{
return log(a0) - f;
}
}

} // namespace detail

template <typename Gen, typename U>
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, boost::math::uintmax_t& max_terms)
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, boost::math::uintmax_t& max_terms, bool logArithmetic=false)
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::continued_fraction_a_impl(g, factor, max_terms);
return detail::continued_fraction_a_impl(g, factor, max_terms, logArithmetic);
}

template <typename Gen, typename U>
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor)
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, const U& factor, bool logArithmetic=false)
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::uintmax_t max_iter = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();
return detail::continued_fraction_a_impl(g, factor, max_iter);
return detail::continued_fraction_a_impl(g, factor, max_iter, logArithmetic);
}

template <typename Gen>
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, bool logArithmetic=false)
noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
#ifndef BOOST_MATH_HAS_GPU_SUPPORT
&& noexcept(std::declval<Gen>()())
Expand All @@ -327,11 +375,11 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type
result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
boost::math::uintmax_t max_iter = (boost::math::numeric_limits<boost::math::uintmax_t>::max)();

return detail::continued_fraction_a_impl(g, factor, max_iter);
return detail::continued_fraction_a_impl(g, factor, max_iter, logArithmetic);
}

template <typename Gen>
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::math::uintmax_t& max_terms)
BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::math::uintmax_t& max_terms, bool logArithmetic=false)
noexcept(BOOST_MATH_IS_FLOAT(typename detail::fraction_traits<Gen>::result_type)
#ifndef BOOST_MATH_HAS_GPU_SUPPORT
&& noexcept(std::declval<Gen>()())
Expand All @@ -344,7 +392,7 @@ BOOST_MATH_GPU_ENABLED inline typename detail::fraction_traits<Gen>::result_type
using result_type = typename traits::result_type;

result_type factor = ldexp(1.0f, 1-bits); // 1 / pow(result_type(2), bits);
return detail::continued_fraction_a_impl(g, factor, max_terms);
return detail::continued_fraction_a_impl(g, factor, max_terms, logArithmetic);
}

} // namespace tools
Expand Down
Loading