1d newton

include/boost/math/distributions/kolmogorov_smirnov.hpp

@@ -382,7 +382,7 @@ inline RealType quantile(const kolmogorov_smirnov_distribution<RealType, Policy>
    std::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
 
    return tools::newton_raphson_iterate(detail::kolmogorov_smirnov_quantile_functor<RealType, Policy>(dist, p),
-           k, RealType(0), boost::math::tools::max_value<RealType>(), get_digits, m);
+           k, RealType(0), RealType(1), get_digits, m);
 } // quantile
 
 template <class RealType, class Policy>
@@ -407,7 +407,7 @@ inline RealType quantile(const complemented2_type<kolmogorov_smirnov_distributio
 
    return tools::newton_raphson_iterate(
            detail::kolmogorov_smirnov_complementary_quantile_functor<RealType, Policy>(dist, p),
-           k, RealType(0), boost::math::tools::max_value<RealType>(), get_digits, m);
+           k, RealType(0), RealType(1), get_digits, m);
 } // quantile (complemented)
 
 template <class RealType, class Policy>

include/boost/math/special_functions/beta.hpp

@@ -1021,18 +1021,15 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de
    BOOST_MATH_ASSERT((p_derivative == 0) || normalised);
 
    if(!(boost::math::isfinite)(a))
-      return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
+      return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol);
    if(!(boost::math::isfinite)(b))
-      return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
-   if(!(boost::math::isfinite)(x))
+      return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol);
+   if (!(0 <= x && x <= 1))
       return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);
 
    if(p_derivative)
       *p_derivative = -1; // value not set.
 
-   if((x < 0) || (x > 1))
-      return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
-
    if(normalised)
    {
       if(a < 0)
@@ -1422,18 +1419,16 @@ T ibeta_derivative_imp(T a, T b, T x, const Policy& pol)
    // start with the usual error checks:
    //
    if (!(boost::math::isfinite)(a))
-      return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol);
+      return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be finite (got a=%1%).", a, pol);
    if (!(boost::math::isfinite)(b))
-      return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol);
-   if (!(boost::math::isfinite)(x))
+      return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be finite (got b=%1%).", b, pol);
+   if (!(0 <= x && x <= 1))
       return policies::raise_domain_error<T>(function, "The argument x to the incomplete beta function must be in [0,1] (got x=%1%).", x, pol);
 
    if(a <= 0)
       return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol);
    if(b <= 0)
       return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol);
-   if((x < 0) || (x > 1))
-      return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol);
    //
    // Now the corner cases:
    //

include/boost/math/special_functions/detail/ibeta_inverse.hpp

@@ -25,23 +25,16 @@ namespace boost{ namespace math{ namespace detail{
 template <class T>
 struct temme_root_finder
 {
-   temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {}
+   temme_root_finder(const T t_, const T a_) : t(t_), a(a_) {
+      const T x_extrema = 1 / (1 + a);
+      BOOST_MATH_ASSERT(0 < x_extrema && x_extrema < 1);
+   }
 
    boost::math::tuple<T, T> operator()(T x)
    {
       BOOST_MATH_STD_USING // ADL of std names
 
       T y = 1 - x;
-      if(y == 0)
-      {
-         T big = tools::max_value<T>() / 4;
-         return boost::math::make_tuple(static_cast<T>(-big), static_cast<T>(-big));
-      }
-      if(x == 0)
-      {
-         T big = tools::max_value<T>() / 4;
-         return boost::math::make_tuple(static_cast<T>(-big), big);
-      }
       T f = log(x) + a * log(y) + t;
       T f1 = (1 / x) - (a / (y));
       return boost::math::make_tuple(f, f1);
@@ -410,6 +403,14 @@ T temme_method_3_ibeta_inverse(T a, T b, T p, T q, const Policy& pol)
    T lower = eta < mu ? cross : 0;
    T upper = eta < mu ? 1 : cross;
    T x = (lower + upper) / 2;
+
+   // Special cases due to lack of numerical precision
+   if (cross == 1) {
+      return 1;
+   } else if (cross == 0) {
+      return 0;
+   }
+
    x = tools::newton_raphson_iterate(
       temme_root_finder<T>(u, mu), x, lower, upper, policies::digits<T, Policy>() / 2);
 #ifdef BOOST_INSTRUMENT