Implement the log of the incomplete beta function #1359

JacobHass8 · 2026-02-10T21:06:52Z

codecov · 2026-02-11T03:57:05Z

Codecov Report

❌ Patch coverage is 52.38095% with 10 lines in your changes missing coverage. Please review.
✅ Project coverage is 92.01%. Comparing base (c42193d) to head (f1c2e6f).

Files with missing lines	Patch %	Lines
include/boost/math/special_functions/beta.hpp	52.38%	10 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1359      +/-   ##
===========================================
- Coverage    95.32%   92.01%   -3.32%     
===========================================
  Files          825      598     -227     
  Lines        67994    56050   -11944     
===========================================
- Hits         64815    51574   -13241     
- Misses        3179     4476    +1297

Files with missing lines	Coverage Δ
include/boost/math/special_functions/beta.hpp	`88.32% <52.38%> (-11.68%)`	⬇️

... and 322 files with indirect coverage changes

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update c42193d...f1c2e6f. Read the comment docs.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

dschmitz89 · 2026-02-11T12:25:37Z

I'd like to ask to reconsider exposing two separate public functions, once ibeta and its log libeta. This would be more consistent with the API style you used so far as you have lgamma, lgamma_p and lgamma_q.

You could have one private function which includes the boolean parameter log carrying out the heavy lifting under the hood (pretty much what you build here). This would then by called by the public functions ibeta and lbeta.

WDYT?

include/boost/math/special_functions/beta.hpp

jzmaddock · 2026-02-11T16:14:18Z

include/boost/math/special_functions/beta.hpp

+      }
+      else
+      {
+         return (invert ? (normalised ? 0 : log(boost::math::beta(a, b, pol))) : -std::numeric_limits<T>::infinity());


overflow_error again, and this opens a can of worms because the regular beta can suffer underflow too, so that probably needs to be logarithmed too.

I noticed there wasn't a lbeta function too. For now, it seems like taking log(beta(a, b)) works pretty well. I can open another PR after to implement lbeta though.

include/boost/math/special_functions/beta.hpp

jzmaddock · 2026-02-11T16:17:12Z

include/boost/math/special_functions/beta.hpp

+         }
+         else
+         {
+            return invert ? log(y) : log(x);


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.

jzmaddock · 2026-02-11T16:23:31Z

test/Jamfile.v2

   [ 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  ]


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.

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.

jzmaddock · 2026-02-11T16:25:28Z

include/boost/math/special_functions/beta.hpp


+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&)


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?

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

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.

jzmaddock · 2026-02-11T16:29:33Z

I'd like to ask to reconsider exposing two separate public functions, once ibeta and its log libeta. This would be more consistent with the API style you used so far as you have lgamma, lgamma_p and lgamma_q.

That's what we have here, but the log functions are called log_ibeta: see my comment on consistency above.

jzmaddock · 2026-02-11T16:50:48Z

This is a good first start, but isn't currently much improvement on calling log(ibeta(a, b, x)), I started looking at what it would take to improve things, and boy, it's going to get complicated :(

Take for example all the places where we invoke the series evaluation, for example:

math/include/boost/math/special_functions/beta.hpp

Line 1330 in c42193d

    
           fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);

and subsequent cases.

That leads us into several failure modes:

Here:

math/include/boost/math/special_functions/beta.hpp

Line 685 in c42193d

    
           result = 0;  // denorms cause overflow in the Lanzos series, result will be zero anyway

a or b is a denorm and the following code blows up unless we back out. This one might not be the end of the world.

This one is probably ignorable too:

math/include/boost/math/special_functions/beta.hpp

Line 697 in c42193d

if (!(boost::math::isfinite)(result))

Then we hit a whole world of badness here:

math/include/boost/math/special_functions/beta.hpp

Line 705 in c42193d

if((l1 > tools::log_min_value<T>())

we want to take the compute by logs branch but not exponentiate the result (which means we need to pass the logarithm parameter down to this function too.

This

math/include/boost/math/special_functions/beta.hpp

Line 744 in c42193d

ibeta_series_t<T> s(a, b, x, result);

is not quite as bad as it looks, the result so far is used as a multiplier to the first value of the series (so that if it's zero or small we terminate straight away). In the log world, we could use 1 as the multiplier, then take the log and add to the result so far.

But now we have a problem - the caller isn't expecting us to return the log of the result - the best I can think of at present is to add a bool result_is_log = false; to the start of ibeta_imp and set it true when the implementation method we've called has returned a log value already... and then fix up all the complicated logic you've just already done!

And then work through the other methods of course...

JacobHass8 · 2026-02-11T17:12:57Z

This is a good first start, but isn't currently much improvement on calling log(ibeta(a, b, x)), I started looking at what it would take to improve things, and boy, it's going to get complicated :(

I'm not sure I agree with this. There are definitely areas that we can improve ibeta_imp to account for logs. However, just taking log(ibeta(a, b, x)) can lead to overflow errors where log_ibeta(a,b,x), as implemented here doesn't. For example, take

typedef double T;

T a = 300.25; 
T b = 2.75;
T x = static_cast<T>(1) / 2048;

T val = boost::math::log_ibeta(a, b, x);
T prev = log(boost::math::ibeta(a, b, x));
std::cout << "log_ibeta(" << a << "," << b << ")=" <<  val << std::endl; // log_ibeta(300.25,2.75)=-2279.78
std::cout << "log(ibeta(" << a << "," << b << "))=" << prev << std::endl; // log(ibeta(300.25,2.75))=-inf

and

typedef double T;

T a = 20; 
T b = 40;
T x = 0.9;

T val = boost::math::log_ibeta(a, b, x);
T prev = log(boost::math::ibeta(a, b, x));
std::cout << "log_ibeta(" << a << "," << b << ")=" <<  val << std::endl; // log_ibeta(20,40)=-1.98955e-26
std::cout << "log(ibeta(" << a << "," << b << "))=" << prev << std::endl; // log(ibeta(300.25,2.75))=0

and lastly,

typedef double T;

T a = 500; 
T b = 500;
T x = 0.9;

T val = boost::math::log_ibeta(a, b, x);
T prev = log(boost::math::ibeta(a, b, x));
std::cout << "log_ibeta(" << a << "," << b << ")=" <<  val << std::endl; // log_ibeta(500,500)=-2.23212e-224
std::cout << "log(ibeta(" << a << "," << b << "))=" << prev << std::endl; // log(ibeta(300.25,2.75))=0

In my testing, log(boost::math::ibeta(a, b, x) under/overflows prematurely in both tails. The implementation in this PR certainly isn't the best solution, but does provide better coverage in both tails.

Maybe I'm missing something easy though. I'll wait to see what you think before pushing ahead any farther on this PR.

dschmitz89 · 2026-02-11T18:22:36Z

I'd like to ask to reconsider exposing two separate public functions, once ibeta and its log libeta. This would be more consistent with the API style you used so far as you have lgamma, lgamma_p and lgamma_q.

That's what we have here, but the log functions are called log_ibeta: see my comment on consistency above.

Ah ok, then nevermind and sorry for the noise!

jzmaddock · 2026-02-12T10:54:53Z

I think you're code succeeds only because you have internal promotion to long double turn ON, if you turn it off (or indeed build with msvc) then your log_ibeta raises an overflow error from log(0).

JacobHass8 added 2 commits February 10, 2026 13:05

Implemented keyword for ibeta_imp

5e32630

Fixed bug with type real_concept

d05919a

JacobHass8 added 2 commits February 10, 2026 20:33

Added lightweight test_ibeta to Jamfile.v2

9b29552

Added grid testing for log_ibeta

4b33948

jzmaddock reviewed Feb 11, 2026

View reviewed changes

include/boost/math/special_functions/beta.hpp Outdated Show resolved Hide resolved

jzmaddock reviewed Feb 11, 2026

View reviewed changes

include/boost/math/special_functions/beta.hpp Outdated Show resolved Hide resolved

jzmaddock reviewed Feb 11, 2026

View reviewed changes

include/boost/math/special_functions/beta.hpp Outdated Show resolved Hide resolved

jzmaddock reviewed Feb 11, 2026

View reviewed changes

Changed return infinity to throwing overflow error

f1c2e6f

Implement the log of the incomplete beta function #1359

Are you sure you want to change the base?

Implement the log of the incomplete beta function #1359

Conversation

JacobHass8 commented Feb 10, 2026

Uh oh!

codecov bot commented Feb 11, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

dschmitz89 commented Feb 11, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jzmaddock commented Feb 11, 2026

Uh oh!

jzmaddock commented Feb 11, 2026

Uh oh!

JacobHass8 commented Feb 11, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dschmitz89 commented Feb 11, 2026

Uh oh!

jzmaddock commented Feb 12, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

codecov bot commented Feb 11, 2026 •

edited

Loading

dschmitz89 commented Feb 11, 2026 •

edited

Loading

JacobHass8 commented Feb 11, 2026 •

edited

Loading