@@ -25,23 +25,16 @@ namespace boost{ namespace math{ namespace detail{
25
25
template <class T >
26
26
struct temme_root_finder
27
27
{
28
- temme_root_finder (const T t_, const T a_) : t(t_), a(a_) {}
28
+ temme_root_finder (const T t_, const T a_) : t(t_), a(a_) {
29
+ const T x_extrema = 1 / (1 + a);
30
+ BOOST_MATH_ASSERT (0 < x_extrema && x_extrema < 1 );
31
+ }
29
32
30
33
boost::math::tuple<T, T> operator ()(T x)
31
34
{
32
35
BOOST_MATH_STD_USING // ADL of std names
33
36
34
37
T y = 1 - x;
35
- if (y == 0 )
36
- {
37
- T big = tools::max_value<T>() / 4 ;
38
- return boost::math::make_tuple (static_cast <T>(-big), static_cast <T>(-big));
39
- }
40
- if (x == 0 )
41
- {
42
- T big = tools::max_value<T>() / 4 ;
43
- return boost::math::make_tuple (static_cast <T>(-big), big);
44
- }
45
38
T f = log (x) + a * log (y) + t;
46
39
T f1 = (1 / x) - (a / (y));
47
40
return boost::math::make_tuple (f, f1);
@@ -410,6 +403,10 @@ T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol)
410
403
T lower = eta < mu ? cross : 0 ;
411
404
T upper = eta < mu ? 1 : cross;
412
405
T x = (lower + upper) / 2 ;
406
+
407
+ // Early exit for cases with numerical precision issues.
408
+ if (cross == 0 || cross == 1 ) { return cross; }
409
+
413
410
x = tools::newton_raphson_iterate (
414
411
temme_root_finder<T>(u, mu), x, lower, upper, policies::digits<T, Policy>() / 2 );
415
412
#ifdef BOOST_INSTRUMENT
0 commit comments