Skip to content

Implement the log of the incomplete beta function #1359

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

Open
JacobHass8 wants to merge 5 commits into
base: develop
Choose a base branch
from
+138 −41
Open

Implement the log of the incomplete beta function #1359

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5e32630
Implemented keyword for ibeta_imp
JacobHass8 Feb 10, 2026
d05919a
Fixed bug with type real_concept
JacobHass8 Feb 11, 2026
9b29552
Added lightweight test_ibeta to Jamfile.v2
JacobHass8 Feb 11, 2026
4b33948
Added grid testing for log_ibeta
JacobHass8 Feb 11, 2026
f1c2e6f
Changed return infinity to throwing overflow error
JacobHass8 Feb 11, 2026
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
98 changes: 89 additions & 9 deletions include/boost/math/special_functions/beta.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1195,7 +1195,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_large_ab(T a, T b, T x, T y, bool invert, bool no
//

template <class T, class Policy>
BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative)
BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative, bool logarithm=false)
{
constexpr auto function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)";
typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
Expand Down Expand Up @@ -1235,12 +1235,32 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b
if(b == 0)
return policies::raise_domain_error<T>(function, "The arguments a and b to the incomplete beta function cannot both be zero, with x=%1%.", x, pol);
if(b > 0)
return static_cast<T>(inv ? 0 : 1);
{
if (!logarithm)
{
return static_cast<T>(inv ? 0 : 1);
}
else
{
return inv ? -policies::raise_overflow_error<T>(function, nullptr, pol) : T(0);
}

}

}
else if(b == 0)
{
if(a > 0)
return static_cast<T>(inv ? 1 : 0);
{
if (!logarithm)
{
return static_cast<T>(inv ? 1 : 0);
}
else
{
return inv ? T(0) : -policies::raise_overflow_error<T>(function, nullptr, pol);
}
}
}
}
else
Expand All @@ -1257,15 +1277,30 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b
{
*p_derivative = (a == 1) ? (T)1 : (a < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
}
return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0));
if (!logarithm)
{
return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0));
}
else
{
return (invert ? (normalised ? 0 : log(boost::math::beta(a, b, pol))) : -policies::raise_overflow_error<T>(function, nullptr, pol));
}
}
if(x == 1)
{
if(p_derivative)
{
*p_derivative = (b == 1) ? T(1) : (b < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2);
}
return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
if (!logarithm)
{
return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
}
else
{
return (invert == 0 ? (normalised ? 0 : log(boost::math::beta(a, b, pol))) : -policies::raise_overflow_error<T>(function, nullptr, pol));
}

}
if((a == 0.5f) && (b == 0.5f))
{
Expand All @@ -1277,7 +1312,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b
T p = invert ? asin(sqrt(y)) / constants::half_pi<T>() : asin(sqrt(x)) / constants::half_pi<T>();
if(!normalised)
p *= constants::pi<T>();
return p;
return logarithm ? log(p) : p;
}
if(a == 1)
{
Expand All @@ -1294,7 +1329,15 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b
{
if(p_derivative)
*p_derivative = 1;
return invert ? y : x;
if (!logarithm)
{
return invert ? y : x;
}
else
{
return invert ? log(y) : log(x);
Copy link
Collaborator

@jzmaddock jzmaddock Feb 11, 2026

Choose a reason for hiding this comment

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

Just a heads up that these conditionals can trip up expression-template-enabled multiprecision types. You might be safer with log(invert ? y : x) but we can fix this up later.

}

}

if(p_derivative)
Expand All @@ -1308,7 +1351,7 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b
p = invert ? T(-boost::math::powm1(x, a, pol)) : T(pow(x, a));
if(!normalised)
p /= a;
return p;
return logarithm ? log(p) : p;
}

if(BOOST_MATH_GPU_SAFE_MIN(a, b) <= 1)
Expand Down Expand Up @@ -1625,7 +1668,15 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b
}
}
}
return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract;
if (!logarithm)
{
return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract;
}
else
{
return invert ? (normalised ? boost::math::log1p(-fract) : log(boost::math::beta(a, b, pol) - fract)) : log(fract);
}

} // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised)

template <class T, class Policy>
Expand All @@ -1634,6 +1685,12 @@ BOOST_MATH_GPU_ENABLED inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool
return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(nullptr));
}

template <class T, class Policy>
BOOST_MATH_GPU_ENABLED inline T log_ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised)
{
return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(nullptr), true);
}

template <class T, class Policy>
BOOST_MATH_GPU_ENABLED T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
{
Expand Down Expand Up @@ -1794,6 +1851,29 @@ BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
return boost::math::ibeta(a, b, x, policies::policy<>());
}

template <class RT1, class RT2, class RT3, class Policy>
BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
log_ibeta(RT1 a, RT2 b, RT3 x, const Policy&)
Copy link
Collaborator

@jzmaddock jzmaddock Feb 11, 2026

Choose a reason for hiding this comment

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

We need to think about consistency here, based on gamma_p/gamma_q/gamma we should just add a single "l" prefix to make libeta, though I admit it doesn't read too well! Thoughts anyone?

Copy link
Member

@mborland mborland Feb 11, 2026

Choose a reason for hiding this comment

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

I'd agree that precedent says we should go with libeta even though it's kind of ugly.

Copy link
Contributor Author

@JacobHass8 JacobHass8 Feb 11, 2026

Choose a reason for hiding this comment

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

I could go either way with libeta or log_ibeta. The reason I didn't go with libeta is it looks like a library named eta. Using libeta would be consistent with lgamma which might be preferable.

In regards to gamma_p/gamma_q, the incomplete beta function already has the complement implemented as ibetac. To be consistent with this, I would suggest using libeta and libetac.

{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<RT1, RT2, RT3>::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, forwarding_policy>(detail::log_ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, true), "boost::math::ibeta<%1%>(%1%,%1%,%1%)");
}
template <class RT1, class RT2, class RT3>
BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
log_ibeta(RT1 a, RT2 b, RT3 x)
{
return boost::math::log_ibeta(a, b, x, policies::policy<>());
}

template <class RT1, class RT2, class RT3, class Policy>
BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
ibetac(RT1 a, RT2 b, RT3 x, const Policy&)
Expand Down
1 change: 1 addition & 0 deletions test/Jamfile.v2
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -263,6 +263,7 @@ test-suite special_fun :
[ run test_gamma_mp.cpp /boost/test//boost_unit_test_framework : : : release <define>TEST=3 [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : <build>no ] : test_gamma_mp_3 ]
[ run test_hankel.cpp /boost/test//boost_unit_test_framework ]
[ run test_hermite.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ]
[ run test_ibeta.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework ]
Copy link
Collaborator

@jzmaddock jzmaddock Feb 11, 2026

Choose a reason for hiding this comment

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

This should be covered by the following lines - we've deliberately split testing ibeta up into smaller chunks because if we test too many types in one go, we get massive object files which fail on cygwin/mingw and run the system out of memory on other CI machines. The "test everything at once" behaviour that is the default is really just there for the convenience of local testing.

Copy link
Contributor Author

@JacobHass8 JacobHass8 Feb 11, 2026

Choose a reason for hiding this comment

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

Okay! I'll go ahead and delete this line then. I added it because when I ran ../../../b2 test_ibeta from the test folder, test_ibeta wasn't found. It output a ton of linking errors.

[ run test_ibeta.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework
: # command line
: # input files
Expand Down
80 changes: 48 additions & 32 deletions test/test_ibeta.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -453,39 +453,55 @@ void test_spots(T)
static_cast<T>(4),
ldexp(static_cast<T>(1 + static_cast<T>(1.0) / 1024), -351)),
static_cast<T>(2.381008060978474962211278613067275529112106932635520021e-212L), tolerance);
BOOST_CHECK_CLOSE(
::boost::math::beta(
static_cast<T>(2),
static_cast<T>(4),
ldexp(static_cast<T>(1 + static_cast<T>(1.0) / 2048), -351)),
static_cast<T>(2.378685692854274898232669682422430136513931911501225435e-212L), tolerance);
BOOST_CHECK_CLOSE(
::boost::math::ibeta(
static_cast<T>(3),
static_cast<T>(5),
ldexp(static_cast<T>(1 + static_cast<T>(15) / 16), -268)),
static_cast<T>(2.386034198603463687323052353589201848077110231388968865e-240L), tolerance);
BOOST_CHECK_CLOSE(
::boost::math::ibeta_derivative(
static_cast<T>(2),
static_cast<T>(4),
ldexp(static_cast<T>(1), -557)),
static_cast<T>(4.23957586190238472641508753637420672781472122471791800210e-167L), tolerance * 4);
BOOST_CHECK_CLOSE(
::boost::math::ibeta_derivative(
static_cast<T>(2),
static_cast<T>(4.5),
ldexp(static_cast<T>(1), -557)),
static_cast<T>(5.24647512910420109893867082626308082567071751558842352760e-167L), tolerance * 20);


T tiny = boost::math::tools::min_value<T>() / 2;
T small = boost::math::tools::epsilon<T>();
if (tiny != 0)
BOOST_CHECK_CLOSE(
::boost::math::beta(
static_cast<T>(2),
static_cast<T>(4),
ldexp(static_cast<T>(1 + static_cast<T>(1.0) / 2048), -351)),
static_cast<T>(2.378685692854274898232669682422430136513931911501225435e-212L), tolerance);
BOOST_CHECK_CLOSE(
::boost::math::ibeta(
static_cast<T>(3),
static_cast<T>(5),
ldexp(static_cast<T>(1 + static_cast<T>(15) / 16), -268)),
static_cast<T>(2.386034198603463687323052353589201848077110231388968865e-240L), tolerance);
BOOST_CHECK_CLOSE(
::boost::math::ibeta_derivative(
static_cast<T>(2),
static_cast<T>(4),
ldexp(static_cast<T>(1), -557)),
static_cast<T>(4.23957586190238472641508753637420672781472122471791800210e-167L), tolerance * 4);
BOOST_CHECK_CLOSE(
::boost::math::ibeta_derivative(
static_cast<T>(2),
static_cast<T>(4.5),
ldexp(static_cast<T>(1), -557)),
static_cast<T>(5.24647512910420109893867082626308082567071751558842352760e-167L), tolerance * 20);

// Test log_ibeta special function
// First check log_ibeta matches ibeta on a grid

T a_vals[] = { 0.25, 1., 5., 10., 50.};
T x_vals[] = { 0.0078125, 0.0625, 0.25, 0.75};

for (T a : a_vals)
{
for (T b : a_vals)
{
BOOST_CHECK_EQUAL(boost::math::ibeta(tiny, small, small), 1);
for (T x : x_vals)
{
BOOST_CHECK_CLOSE(exp(::boost::math::log_ibeta(a, b, x)), ::boost::math::ibeta(a, b, x), tolerance);
}
}
BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast<T>(2), static_cast<T>(1), static_cast<T>(0)), 0);
BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast<T>(1), static_cast<T>(2), static_cast<T>(0)), 0);
}

T tiny = boost::math::tools::min_value<T>() / 2;
T small = boost::math::tools::epsilon<T>();
if (tiny != 0)
{
BOOST_CHECK_EQUAL(boost::math::ibeta(tiny, small, small), 1);
}
BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast<T>(2), static_cast<T>(1), static_cast<T>(0)), 0);
BOOST_CHECK_EQUAL(boost::math::ibeta(static_cast<T>(1), static_cast<T>(2), static_cast<T>(0)), 0);
}

Loading