Merge pull request #1025 from boostorg/pr1007

Improve ibeta power term calculation in the Sterling case

Author: jzmaddock

include/boost/math/distributions/detail/inv_discrete_quantile.hpp

@@ -147,7 +147,7 @@ typename Dist::value_type
    // we're assuming that "guess" is likely to be accurate
    // to the nearest int or so:
    //
-   else if(adder != 0)
+   else if((adder != 0) && (a + adder != a))
    {
       //
       // If we're looking for a large result, then bump "adder" up

include/boost/math/special_functions/beta.hpp

@@ -473,20 +473,108 @@ T ibeta_power_terms(T a,
    if ((shift_a == 0) && (shift_b == 0))
    {
       T power1, power2;
+      bool need_logs = false;
       if (a < b)
       {
-         power1 = pow((x * y * c * c) / (a * b), a);
-         power2 = pow((y * c) / b, b - a);
+         BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
+         {
+            power1 = pow((x * y * c * c) / (a * b), a);
+            power2 = pow((y * c) / b, b - a);
+         }
+         else
+         {
+            // We calculate these logs purely so we can check for overflow in the power functions
+            T l1 = log((x * y * c * c) / (a * b));
+            T l2 = log((y * c) / b);
+            if ((l1 * a > tools::log_min_value<T>()) && (l1 * a < tools::log_max_value<T>()) && (l2 * (b - a) < tools::log_max_value<T>()) && (l2 * (b - a) > tools::log_min_value<T>()))
+            {
+               power1 = pow((x * y * c * c) / (a * b), a);
+               power2 = pow((y * c) / b, b - a);
+            }
+            else
+            {
+               need_logs = true;
+            }
+         }
       }
       else
       {
-         power1 = pow((x * y * c * c) / (a * b), b);
-         power2 = pow((x * c) / a, a - b);
+         BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
+         {
+            power1 = pow((x * y * c * c) / (a * b), b);
+            power2 = pow((x * c) / a, a - b);
+         }
+         else
+         {
+            // We calculate these logs purely so we can check for overflow in the power functions
+            T l1 = log((x * y * c * c) / (a * b)) * b;
+            T l2 = log((x * c) / a) * (a - b);
+            if ((l1 * a > tools::log_min_value<T>()) && (l1 * a < tools::log_max_value<T>()) && (l2 * (b - a) < tools::log_max_value<T>()) && (l2 * (b - a) > tools::log_min_value<T>()))
+            {
+               power1 = pow((x * y * c * c) / (a * b), b);
+               power2 = pow((x * c) / a, a - b);
+            }
+            else
+               need_logs = true;
+         }
       }
-      if (!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2))
+      BOOST_IF_CONSTEXPR(std::numeric_limits<T>::has_infinity)
       {
-         // We have to use logs :(
-         return prefix * exp(a * log(x * c / a) + b * log(y * c / b)) * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol));
+         if (!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2))
+         {
+            need_logs = true;
+         }
+      }
+      if (need_logs)
+      {
+         //
+         // We want:
+         //
+         // (xc / a)^a (yc / b)^b
+         //
+         // But we know that one or other term will over / underflow and combining the logs will be next to useless as that will cause significant cancellation.
+         // If we assume b > a and express z ^ b as(z ^ b / a) ^ a with z = (yc / b) then we can move one power term inside the other :
+         //
+         // ((xc / a) * (yc / b)^(b / a))^a
+         //
+         // However, we're not quite there yet, as the term being exponentiated is quite likely to be close to unity, so let:
+         //
+         // xc / a = 1 + (xb - ya) / a
+         //
+         // analogously let :
+         //
+         // 1 + p = (yc / b) ^ (b / a) = 1 + expm1((b / a) * log1p((ya - xb) / b))
+         //
+         // so putting the two together we have :
+         //
+         // exp(a * log1p((xb - ya) / a + p + p(xb - ya) / a))
+         //
+         // Analogously, when a > b we can just swap all the terms around.
+         // 
+         // Finally, there are a few cases (x or y is unity) when the above logic can't be used
+         // or where there is no logarithmic cancellation and accuracy is better just using
+         // the regular formula:
+         //
+         T xc_a = x * c / a;
+         T yc_b = y * c / b;
+         if ((x == 1) || (y == 1) || (fabs(xc_a - 1) > 0.25) || (fabs(yc_b - 1) > 0.25))
+         {
+            // The above logic fails, the result is almost certainly zero:
+            power1 = exp(log(xc_a) * a + log(yc_b) * b);
+            power2 = 1;
+         }
+         else if (b > a)
+         {
+            T p = boost::math::expm1((b / a) * boost::math::log1p((y * a - x * b) / b));
+            power1 = exp(a * log1p((x * b - y * a) / a + p * (x * c / a)));
+            power2 = 1;
+         }
+         else
+         {
+            T p = boost::math::expm1((a / b) * boost::math::log1p((x * b - y * a) / a));
+            power1 = exp(b * log1p((y * a - x * b) / b + p * (y * c / b)));
+            power2 = 1;
+         }
       }
       return prefix * power1 * power2 * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol));
    }

test/test_binomial.cpp

@@ -716,10 +716,7 @@ void test_spots(RealType T)
 
    check_out_of_range<boost::math::binomial_distribution<RealType> >(1, 1); // (All) valid constructor parameter values.
 
-   // TODO: Generic ibeta_power_terms has accuracy issue when long
-   // double is not precise enough, causing overflow in this case.
-   if(!std::is_same<RealType, real_concept>::value ||
-      sizeof(boost::math::concepts::real_concept_base_type) > 8)
+   // See bug reported here: https://github.com/boostorg/math/pull/1007
    {
       using namespace boost::math::policies;
       typedef policy<discrete_quantile<integer_round_outwards> > Policy;