Fixes ibeta for large arguments

[JacobHass8]

Closes #1361 and scipy/scipy#24566.

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

[JacobHass8]

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

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

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]

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]

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]

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