Fixes ibeta for large arguments by JacobHass8 · Pull Request #1363 · boostorg/math

JacobHass8 · 2026-02-17T18:47:17Z

Only the expansion for $\eta$ has been added. When to use the erf expansion and relevant tests need to be added.

JacobHass8 · 2026-02-17T19:45:42Z

@jzmaddock is there a way to see if the continued fraction is converging quickly or not? Maybe a ratio of subsequent terms or a change in the value of the continued fraction?

In this spirit, the easiest think would just be to add a try/catch statement as done here

try
{
   fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
}
// If series converges slowly, fall back to erf approximation
catch (boost::math::evaluation_error& e)
{
   fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
}

Maybe this will have unintended consequences though?

mborland · 2026-02-17T20:02:13Z

@jzmaddock is there a way to see if the continued fraction is converging quickly or not? Maybe a ratio of subsequent terms or a change in the value of the continued fraction?

Non-programatically you can define BOOST_MATH_INSTRUMENT and you'll get all that information dumped into your terminal.

In this spirit, the easiest think would just be to add a try/catch statement as done here
try
{
   fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
}
// If series converges slowly, fall back to erf approximation
catch (boost::math::evaluation_error& e)
{
   fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
}
Maybe this will have unintended consequences though?

We test and run on a number of noexcept as you can probably see from the CI. Right now in a noexcept environment it looks like you can't really tell what happened. Perhaps changing the signature to:

BOOST_MATH_GPU_ENABLED inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative, boost::math::uintmax_t iters = nullptr);

Since max_terms is written to by reference in the continued fraction impl. Then you could do something a la

uintmax_t iters;
fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative, &iters);
If (too many iters)
{
   fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
}

Defaulting to, and checking the nullptr before writing after the throwing location would probably keep this change from affecting other code.

JacobHass8 · 2026-02-17T20:26:54Z

Since max_terms is written to by reference in the continued fraction impl. Then you could do something a la
uintmax_t iters;
fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative, &iters);
If (too many iters)
{
   fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
}
Defaulting to, and checking the nullptr before writing after the throwing location would probably keep this change from affecting other code.

The issue is that before returning, ibeta_fraction2 runs

boost::math::policies::check_series_iterations<T>("boost::math::ibeta", max_terms, pol);

which throws the evaluation error. I think that your solution would work, but then we would have to get rid of the check above in ibeta_fraction2 and run it after every call of ibeta_fraction2. Does that make sense? I'm not sure if I explained that well.

Maybe there's an easier solution I'm not seeing though. Here's the function in question:

template <class T, class Policy>
BOOST_MATH_GPU_ENABLED inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative)
{
   typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
   BOOST_MATH_STD_USING
   T result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
   if(p_derivative)
   {
      *p_derivative = result;
      BOOST_MATH_ASSERT(*p_derivative >= 0);
   }
   if(result == 0)
      return result;

   ibeta_fraction2_t<T> f(a, b, x, y);
   boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations<Policy>();
   T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>(), max_terms);
   boost::math::policies::check_series_iterations<T>("boost::math::ibeta", max_terms, pol); // Want to catch this error! 
   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
   BOOST_MATH_INSTRUMENT_VARIABLE(result);
   return result / fract;
}

mborland · 2026-02-17T20:43:11Z

In this case non-convergence isn't exceptional, it's somewhat expected right? Passing a modified policy to ibeta_fraction2 like ignore_error would keep you from throwing at that check. Then this works:

template <class T, class Policy>
BOOST_MATH_GPU_ENABLED inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative, uintmax_t* iters = nullptr)
{
   typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
   BOOST_MATH_STD_USING
   T result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
   if(p_derivative)
   {
      *p_derivative = result;
      BOOST_MATH_ASSERT(*p_derivative >= 0);
   }
   if(result == 0)
      return result;

   ibeta_fraction2_t<T> f(a, b, x, y);
   boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations<Policy>();
   T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>(), max_terms);
   boost::math::policies::check_series_iterations<T>("boost::math::ibeta", max_terms, pol); // Nothing happens
   if(iters && max_iter >= policies::get_max_series_iterations<Policy>()) *iters = max_iter;
   BOOST_MATH_INSTRUMENT_VARIABLE(fract);
   BOOST_MATH_INSTRUMENT_VARIABLE(result);
   return result / fract;
}

JacobHass8 · 2026-02-17T22:04:33Z

In this case non-convergence isn't exceptional, it's somewhat expected right? Passing a modified policy to ibeta_fraction2 like ignore_error would keep you from throwing at that check. Then this works:

I think that should work! I'm not very familiar with policies though. How would I modify the policy to ignore_error on within ibeta_imp?

mborland · 2026-02-17T23:05:50Z

In this case non-convergence isn't exceptional, it's somewhat expected right? Passing a modified policy to ibeta_fraction2 like ignore_error would keep you from throwing at that check. Then this works:

I think that should work! I'm not very familiar with policies though. How would I modify the policy to ignore_error on within ibeta_imp?

using ignore_iterations_policy = typename boost::math::policies::normalise<
        Policy,
        boost::math::policies::max_series_iterations_error<boost::math::policies::ignore_error>
    >::type;

    ignore_iterations_policy new_policy;

And then use the new_policy object to call ibeta_fraction2

JacobHass8 · 2026-02-18T04:00:13Z

using ignore_iterations_policy = typename boost::math::policies::normalise<
        Policy,
        boost::math::policies::max_series_iterations_error<boost::math::policies::ignore_error>
    >::type;

    ignore_iterations_policy new_policy;

I implemented this, but am getting a ton of errors when running all the ibeta tests. The issue is that there are a number of tests where iter >= policies::get_max_series_iterations<Policy>(). In these cases, boost::math::policies::check_series_iterations<T>("boost::math::ibeta", max_terms, pol) doesn't throw an error, but the series has converged and the answer should be accepted (at least according to the test data). Just checking iter >= policies::get_max_series_iterations<Policy>() will cause the erf approximation to be used in these test cases which is less accurate. I think that I had a bug which was causing the number of iterations to be artificially inflated for the tests. I'm finding the number of iterations are actually much more reasonable now.

I'm not sure what the solution here should be. However, the try/catch statement seems the simplest and all of the test cases pass with this implementation.

jzmaddock · 2026-02-18T11:04:37Z

I'll try and look into the policy thing later.

jzmaddock · 2026-02-18T13:45:30Z

It's a bit of a cop-out, but this patch passes tests locally for me:

diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp
index a277c65e8..569a8fb23 100644
--- a/include/boost/math/special_functions/beta.hpp
+++ b/include/boost/math/special_functions/beta.hpp
@@ -1584,46 +1584,39 @@ BOOST_MATH_GPU_ENABLED T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, b
       else
       {
          // a and b both large:
-         bool use_asym = false;
          T ma = BOOST_MATH_GPU_SAFE_MAX(a, b);
-         T xa = ma == a ? x : y;
          T saddle = ma / (a + b);
-         T powers = 0;
-         if ((ma > 1e-5f / tools::epsilon<T>()) && (ma / BOOST_MATH_GPU_SAFE_MIN(a, b) < (xa < saddle ? 2 : 15)))
-         {
-            if (a == b)
-               use_asym = true;
-            else
-            {
-               powers = exp(log(x / (a / (a + b))) * a + log(y / (b / (a + b))) * b);
-               if (powers < tools::epsilon<T>())
-                  use_asym = true;
-            }
-         }
-         if(use_asym)
+         if ((ma > 1e-5f / tools::epsilon<T>()) && (fabs(saddle - x) < 1e-12) && (1 - saddle > 0.125))
          {
             fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
-            if (fract * tools::epsilon<T>() < powers)
-            {
-               // Erf approximation failed, correction term is too large, fall back:
-               fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
-            }
-            else
-               invert = false;
+            invert = false;
          }
          else
          {
-            try
+            typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
+            T local_result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol);
+            if (p_derivative)
             {
-               fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
+               *p_derivative = local_result;
+               BOOST_MATH_ASSERT(*p_derivative >= 0);
             }
-            // If series converges slowly, fall back to erf approximation
-            catch (boost::math::evaluation_error& e)
+            if (local_result != 0)
             {
-               fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
+               ibeta_fraction2_t<T> f(a, b, x, y);
+               boost::math::uintmax_t max_terms = boost::math::policies::get_max_series_iterations<Policy>();
+               T local_fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>(), max_terms);
+               if (max_terms >= boost::math::policies::get_max_series_iterations<Policy>())
+               {
+                  // Continued fraction failed, fall back to asymptotic expansion:
+                  fract = ibeta_large_ab(a, b, x, y, invert, normalised, pol);
+                  invert = false;
+               }
+               else
+                  fract = local_result / local_fract;
             }
-         }  
-         BOOST_MATH_INSTRUMENT_VARIABLE(fract);
+            else
+               fract = 0;
+         }
       }
    }
    if(p_derivative)

JacobHass8 · 2026-02-18T23:38:42Z

It's a bit of a cop-out, but this patch passes tests locally for me:

If I'm understanding your code correctly, you've simply just written out all the logic from ibeta_fraction2 without the series convergence test. Instead, you check if the maximum number of iterations has been reached and use the ibeta_large_ab approximation in this case.

This is kind of what I had in mind, but wasn't implementing well. I think this is a reasonable solution! We've copied the code from ibeta_fraction2, but this makes more sense to me than defining a new policy.

JacobHass8 · 2026-02-18T23:55:39Z

As far as test data goes, how should I go about testing this region since Mathematica hasn't been able to produce results. The simplest test is to make sure that ibeta is always increasing as $x$ increases. I was thinking test that ibeta increases in the range $x=\mathrm{saddle} \pm 10^{-10,\ldots, -15}$. This would ensure that we traverse the intersection of the continued fraction and erf approximation.

I could also implement the erf approximation in Mathematica and test that we're close to these values. We wouldn't be comparing to "true" data, but maybe this is the closest that anyone can get?

codecov · 2026-02-19T00:57:25Z

Codecov Report

❌ Patch coverage is 96.87500% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 95.43%. Comparing base (c42193d) to head (1611bae).
⚠️ Report is 27 commits behind head on develop.

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

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1363      +/-   ##
===========================================
+ Coverage    95.32%   95.43%   +0.11%     
===========================================
  Files          825      761      -64     
  Lines        67994    62344    -5650     
===========================================
- Hits         64815    59498    -5317     
+ Misses        3179     2846     -333

Files with missing lines	Coverage Δ
test/test_ibeta.hpp	`99.09% <100.00%> (+0.08%)`	⬆️
include/boost/math/special_functions/beta.hpp	`99.26% <95.45%> (-0.74%)`	⬇️

... and 82 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...1611bae. 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-19T08:31:55Z

Hey everyone,
just wanted to mention that not every stress test SciPy users come up with warrant a definite "fix" in boost. I do not want to leave boost maintainers with the impression that they get bombarded with edge case handling by us. I also forwarded a few of such issues to you guys in the past but please, do not feel pressurised to respond and fix all of them if they cause more trouble than they are worth.

If you get the feeling that we start playing whack-a-mole in the sense that one fix causes trouble in other cases in this regime, it should be fine if you declare certain inputs infeasible @jzmaddock .

jzmaddock · 2026-02-19T16:50:00Z

If I'm understanding your code correctly, you've simply just written out all the logic from ibeta_fraction2 without the series convergence test. Instead, you check if the maximum number of iterations has been reached and use the ibeta_large_ab approximation in this case.

Correct, we can stop a call to ibeta_fraction2 from throwing at this location, but we don't get back the iteration count only the result so far, so a bit of cut and paste seemed the best solution.

As far as test data goes, how should I go about testing this region since Mathematica hasn't been able to produce results. The simplest test is to make sure that ibeta is always increasing as x increases. I was thinking test that ibeta increases in the range x = saddle ± 10 − 10 , … , − 15 . This would ensure that we traverse the intersection of the continued fraction and erf approximation.

That's sensible, at least we should do a bit of due diligence. I suspect the error rates are so large in this area that we'll get a lot of noise, and a monotonically increasing result is probably too much to ask for. What we could do is generate a scatter plot for a handful of values of a and b, and check things look sensible as x goes through the saddle point and we switch methods.

I could also implement the erf approximation in Mathematica and test that we're close to these values. We wouldn't be comparing to "true" data, but maybe this is the closest that anyone can get?

Hmm, the problem with that is we don't necessarily have a good handle on what the error rate inherent to the approximation is?

Hey everyone,
just wanted to mention that not every stress test SciPy users come up with warrant a definite "fix" in boost. I do not want to leave boost maintainers with the impression that they get bombarded with edge case handling by us. I also forwarded a few of such issues to you guys in the past but please, do not feel pressurised to respond and fix all of them if they cause more trouble than they are worth.

If you get the feeling that we start playing whack-a-mole in the sense that one fix causes trouble in other cases in this regime, it should be fine if you declare certain inputs infeasible @jzmaddock .

For sure, and we're definitely on the limit with this one, certainly any results in this domain are likely to be very approximate, but I also think we can do a little bit better than at present.

jzmaddock · 2026-02-19T17:35:12Z

I tried a quick scatter chart (google gemini AI which I don't necessarily trust!!) for a = 10000000272564224, b = 9965820922822656 and distance from saddle point:

Looks broadly OK in this case?

JacobHass8 · 2026-02-19T18:30:57Z

I tried a quick scatter chart (google gemini AI which I don't necessarily trust!!) for a = 10000000272564224, b = 9965820922822656 and distance from saddle point:

Was this made with the files in the current PR? I'd say that looks pretty good too! I can try to write up a test to ensure the values are first off not nan and are increasing.

jzmaddock · 2026-02-19T19:28:28Z

Was this made with the files in the current PR?

This PR with the patch I suggested applied. It's just one data point of course, but it's encouraging at least. The values do seem to be monotonically increasing too, though I suspect if you zoomed in on a small enough interval that might change as you see more "noise". This is probably true of all the special functions though.

JacobHass8 · 2026-02-19T20:15:35Z

This PR with the patch I suggested applied. It's just one data point of course, but it's encouraging at least. The values do seem to be monotonically increasing too, though I suspect if you zoomed in on a small enough interval that might change as you see more "noise". This is probably true of all the special functions though.

I added some very basic testing to see if ibeta is monotonically increasing. It's more of a sanity check, but maybe that is the best that we can do? Note, the tests flat out fail for floats but that's probably acceptable.

// Bug testing for large a,b and x close to a / (a+b). See PR 1363.
// The values for a,b are just too large for floats to handle.
if (!std::is_same<T, float>::value)
{
   T a_values[2] = {10000000272564224, 3.1622776601699636e+16};
   T b_values[2] = {9965820922822656, 3.130654883566682e+18};
   T delta[8] = {0, 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10, 1e-9};
   T a, b, x;
   for (unsigned int i=0; i<2; i++){
      a = a_values[i];
      b = b_values[i];
      x = a / (a+b); // roughly the median of ibeta

      
      for (unsigned int j=0; j < 7; j++)
      {
         BOOST_CHECK(boost::math::ibeta(a, b, x + delta[j+1]) > boost::math::ibeta(a, b, x + delta[j]));
         BOOST_CHECK(boost::math::ibeta(a, b, x - delta[j]) > boost::math::ibeta(a, b, x - delta[j+1]));  
      }
   }
}

JacobHass8 added 2 commits February 17, 2026 10:44

Added nu expansion

638bfc4

Added try/catch statement

3c176f1

Accepted review comments to copy code from ibeta_fraction2

8b5b1fd

masoudahg00 mentioned this pull request Feb 19, 2026

Fix type passed into BOOST_MATH_NOEXCEPT() in specialization of epsilon(). #14

Merged

Added some sparse testing

12047b6

Relaxed tests

1611bae

Conversation

JacobHass8 commented Feb 17, 2026

Uh oh!

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

Uh oh!

Uh oh!

mborland commented Feb 17, 2026

Uh oh!

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

Uh oh!

Uh oh!

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

Uh oh!

Uh oh!

JacobHass8 commented Feb 17, 2026

Uh oh!

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

Uh oh!

Uh oh!

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

Uh oh!

Uh oh!

jzmaddock commented Feb 18, 2026

Uh oh!

jzmaddock commented Feb 18, 2026

Uh oh!

JacobHass8 commented Feb 18, 2026

Uh oh!

JacobHass8 commented Feb 18, 2026

Uh oh!

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

Uh oh!

Codecov Report

Uh oh!

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

Uh oh!

Uh oh!

jzmaddock commented Feb 19, 2026

Uh oh!

jzmaddock commented Feb 19, 2026

Uh oh!

JacobHass8 commented Feb 19, 2026

Uh oh!

jzmaddock commented Feb 19, 2026

Uh oh!

JacobHass8 commented Feb 19, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Comments

JacobHass8 commented Feb 17, 2026 •

edited

Loading

JacobHass8 commented Feb 17, 2026 •

edited

Loading

mborland commented Feb 17, 2026 •

edited

Loading

mborland commented Feb 17, 2026 •

edited

Loading

JacobHass8 commented Feb 18, 2026 •

edited

Loading

codecov bot commented Feb 19, 2026 •

edited

Loading

dschmitz89 commented Feb 19, 2026 •

edited

Loading