Add lbeta

doc/fp_utilities/ulps_plots.qbk

@@ -20,7 +20,7 @@
                   size_t samples = 10000, bool perturb_abscissas = false, int random_seed = -1);
         //
         // Set a clip value to restrict the range of ULP's shown, values outside [-clip,clip]
-        // will be shown at the clip boundary and in a different color (red by default):
+        // will be shown at the clip boundary as crosses, in the color of their function by default:
         //
         ulps_plot& clip(PreciseReal clip);
         //
@@ -44,8 +44,8 @@
         //
         ulps_plot& font_color(std::string const & font_color);
         //
-        // Sets the color of points which are cropped, or set to an empty
-        // string to hide them completely (default is "red"):
+        // Sets one color for all points which are cropped, or set to an empty
+        // string to hide them completely (default is the color of each point's function):
         //
         ulps_plot& crop_color(std::string const & font_color);
         //

doc/sf/beta.qbk

@@ -13,6 +13,12 @@
 
    template <class T1, class T2, class ``__Policy``>
    BOOST_MATH_GPU_ENABLED ``__sf_result`` beta(T1 a, T2 b, const ``__Policy``&);
+
+   template <class T1, class T2>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` lbeta(T1 a, T2 b);
+
+   template <class T1, class T2, class ``__Policy``>
+   BOOST_MATH_GPU_ENABLED ``__sf_result`` lbeta(T1 a, T2 b, const ``__Policy``&);
 
    }} // namespaces
 
@@ -29,6 +35,11 @@ The beta function is defined by:
 The return type of these functions is computed using the __arg_promotion_rules
 when T1 and T2 are different types.
 
+`lbeta` returns the natural logarithm of the beta function, log B(a, b).
+It stays finite and accurate where `beta` underflows, for example when both /a/ and /b/ are large,
+and where lgamma(a) + lgamma(b) - lgamma(a + b) would lose every digit to cancellation,
+when one argument is much larger than the other.
+
 [h4 Accuracy]
 
 The following table shows peak errors for various domains of input arguments,
@@ -69,9 +80,21 @@ immediately cuts the maximum error in half:
 This may not be the final solution, but it is very competitive compared to
 other implementation methods.
 
+`lbeta` takes the logarithm of this form term by term. With /a/ >= /b/, /c = a + b/ and /xgh = x + g - 1\/2/,
+and since /agh = cgh - b/, the two power terms merge:
+
+log B(a, b) = log(L(a) L(b) / L(c)) + (a - 1\/2) log1p(-b \/ cgh) + b log(bgh \/ cgh) + (1 - log bgh) \/ 2
+
+where /L/ is the scaled Lanczos sum. None of these terms cancel when one argument is large, and none underflow
+when both are. In double precision the result is within 2 ulps where |log B| >= 1, and within 1 epsilon
+absolute error elsewhere, for arguments across the whole exponent range.
+Without a __lanczos approximation, `lbeta` takes logarithms in the generic method below, except that when one
+argument is small it uses lgamma and __tgamma_delta_ratio directly, and when both lie between 1 and the
+threshold for Stirling's approximation it takes the logarithm of `beta` itself, which cannot overflow or underflow there.
+
 The generic implementation - where no __lanczos approximation is available - is
 implemented in a very similar way to the generic version of the gamma function
-by means of Sterling's approximation.
+by means of Stirling's approximation.
 Again in order to avoid numerical overflow the power terms that prefix the series
 are collected together
 

include/boost/math/special_functions/beta.hpp

@@ -205,6 +205,161 @@ BOOST_MATH_GPU_ENABLED T beta_imp(T a, T b, const lanczos::undefined_lanczos& l,
 } // template <class T>T beta_imp(T a, T b, const lanczos::undefined_lanczos& l)
 #endif
 
+//
+// Logarithm of Beta(a,b), from the same Lanczos approximation as beta_imp.
+// With a >= b, c = a + b and xgh = x + g - 0.5, the log of beta_imp's product is
+//
+//    log(S(a) S(b) / S(c)) + (a - 0.5) log1p(-b / cgh) + b log(bgh / cgh) + (1 - log(bgh)) / 2,
+//
+// where S is the scaled Lanczos sum, using agh = cgh - b to merge beta_imp's two powers.
+// No large terms cancel, unlike lgamma(a) + lgamma(b) - lgamma(a + b) when one argument is large,
+// and nothing underflows, unlike log(beta(a, b)) when both are.
+//
+template <class T, class Lanczos, class Policy>
+BOOST_MATH_GPU_ENABLED T lbeta_imp(T a, T b, const Lanczos&, const Policy& pol)
+{
+   BOOST_MATH_STD_USING  // for ADL of std names
+
+   // Reject NaN before any comparison sends it down a branch meant for ordered arguments:
+   if((BOOST_MATH_NAMESPACE::isnan)(a) || (BOOST_MATH_NAMESPACE::isnan)(b))
+      return policies::raise_domain_error<T>("boost::math::lbeta<%1%>(%1%,%1%)", "The arguments to the beta function must not be NaN (got %1%).", (BOOST_MATH_NAMESPACE::isnan)(a) ? a : b, pol);
+   if(a <= 0)
+      return policies::raise_domain_error<T>("boost::math::lbeta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol);
+   if(b <= 0)
+      return policies::raise_domain_error<T>("boost::math::lbeta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol);
+   // log(beta(a, b)) tends to -infinity when either argument does:
+   if((BOOST_MATH_NAMESPACE::isinf)(a) || (BOOST_MATH_NAMESPACE::isinf)(b))
+      return -policies::raise_overflow_error<T>("boost::math::lbeta<%1%>(%1%,%1%)", nullptr, pol);
+
+   T c = a + b;
+
+   // Special cases, as in beta_imp. For tiny b, log(beta) = -log(b) - b digamma(a) + O(b); unlike in
+   // beta_imp, the b term is not negligible next to -log(b) when a is huge, and log(a) is close enough to
+   // digamma(a) wherever c == a:
+   if((c == a) && (b < tools::epsilon<T>()))
+      return -log(b) - b * log(a);
+   else if((c == b) && (a < tools::epsilon<T>()))
+      return -log(a) - a * log(b);
+   if(b == 1)
+      return -log(a);
+   else if(a == 1)
+      return -log(b);
+   else if(c < tools::epsilon<T>())
+      return log(c) - log(a) - log(b);
+
+   if(a < b)
+   {
+      BOOST_MATH_GPU_SAFE_SWAP(a, b);
+   }
+
+   T bgh = static_cast<T>(b + Lanczos::g() - 0.5f);
+   T cgh = static_cast<T>(c + Lanczos::g() - 0.5f);
+   T result = log(Lanczos::lanczos_sum_expG_scaled(a) * (Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c)));
+   result += (a - T(0.5)) * BOOST_MATH_NAMESPACE::log1p(-b / cgh, pol);
+   result += b * log(bgh / cgh);
+   result += (1 - log(bgh)) / 2;
+   return result;
+} // template <class T, class Lanczos> lbeta_imp(T a, T b, const Lanczos&)
+
+//
+// Logarithm of the generic Beta(a,b) without Lanczos approximation support:
+//
+#ifndef BOOST_MATH_HAS_GPU_SUPPORT
+template <class T, class Policy>
+BOOST_MATH_GPU_ENABLED T lbeta_imp(T a, T b, const lanczos::undefined_lanczos& l, const Policy& pol)
+{
+   BOOST_MATH_STD_USING
+
+   // Reject NaN before any comparison sends it down a branch meant for ordered arguments:
+   if((BOOST_MATH_NAMESPACE::isnan)(a) || (BOOST_MATH_NAMESPACE::isnan)(b))
+      return policies::raise_domain_error<T>("boost::math::lbeta<%1%>(%1%,%1%)", "The arguments to the beta function must not be NaN (got %1%).", (BOOST_MATH_NAMESPACE::isnan)(a) ? a : b, pol);
+   if(a <= 0)
+      return policies::raise_domain_error<T>("boost::math::lbeta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol);
+   if(b <= 0)
+      return policies::raise_domain_error<T>("boost::math::lbeta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol);
+   // log(beta(a, b)) tends to -infinity when either argument does:
+   if((BOOST_MATH_NAMESPACE::isinf)(a) || (BOOST_MATH_NAMESPACE::isinf)(b))
+      return -policies::raise_overflow_error<T>("boost::math::lbeta<%1%>(%1%,%1%)", nullptr, pol);
+
+   // Fix the argument order so that lbeta(a, b) == lbeta(b, a) exactly:
+   if (a < b)
+   {
+      BOOST_MATH_GPU_SAFE_SWAP(a, b);
+   }
+   const T c = a + b;
+
+   // Special cases, as in the Lanczos version:
+   if ((c == a) && (b < tools::epsilon<T>()))
+      return -log(b) - b * log(a);
+   else if ((c == b) && (a < tools::epsilon<T>()))
+      return -log(a) - a * log(b);
+   if (b == 1)
+      return -log(a);
+   else if (a == 1)
+      return -log(b);
+   else if (c < tools::epsilon<T>())
+      return log(c) - log(a) - log(b);
+
+   // Regular cases start here:
+   const T min_sterling = minimum_argument_for_bernoulli_recursion<T>();
+
+   long shift_a = 0;
+   long shift_b = 0;
+
+   if(a < min_sterling)
+      shift_a = 1 + ltrunc(min_sterling - a);
+   if(b < min_sterling)
+      shift_b = 1 + ltrunc(min_sterling - b);
+   long shift_c = shift_a + shift_b;
+
+   if ((shift_a == 0) && (shift_b == 0))
+   {
+      // log of pow(a / c, a) * pow(b / c, b), with log1p for the larger argument, whose ratio is near one:
+      T power_terms = a < b ? a * log(a / c) + b * BOOST_MATH_NAMESPACE::log1p(-a / c, pol) : a * BOOST_MATH_NAMESPACE::log1p(-b / c, pol) + b * log(b / c);
+      return power_terms + log(scaled_tgamma_no_lanczos(a, pol) * (scaled_tgamma_no_lanczos(b, pol) / scaled_tgamma_no_lanczos(c, pol)));
+   }
+   else if ((a < 1) && (b < 1))
+   {
+      return BOOST_MATH_NAMESPACE::lgamma(a, pol) + (BOOST_MATH_NAMESPACE::lgamma(b, pol) - BOOST_MATH_NAMESPACE::lgamma(c, pol));
+   }
+   else if((a < 1) || (b < 1) || (shift_a == 0) || (shift_b == 0))
+   {
+      // One argument, say a, is small: log(gamma(a) gamma(b) / gamma(a + b)) directly, since shifting
+      // a up and back would cancel terms of size shift_a * log(c) when b is large. The ratio of gammas
+      // is about b^-a, which only a narrow exponent range can underflow; then fall back to recursion.
+      T small_arg = a < b ? a : b;
+      T large_arg = a < b ? b : a;
+      T ratio = BOOST_MATH_NAMESPACE::tgamma_delta_ratio(large_arg, small_arg, pol);
+      if ((ratio >= tools::min_value<T>()) || (small_arg < 1))
+         return BOOST_MATH_NAMESPACE::lgamma(small_arg, pol) + log(ratio);
+   }
+   else
+   {
+      // Both arguments lie in [1, min_sterling), so beta itself can neither overflow nor underflow,
+      // and its recursion multiplies rather than cancels: taking the log only reduces its relative error.
+      return log(beta_imp(a, b, l, pol));
+   }
+   {
+      // Only reached when one argument is huge and the gamma ratio above underflows.
+      T result = lbeta_imp(T(a + shift_a), T(b + shift_b), l, pol);
+      //
+      // Recursion, summed as logs: as a product it overflows when one argument is huge:
+      //
+      for (long i = 0; i < shift_c; ++i)
+      {
+         T ratio = c + i;
+         if (i < shift_a)
+            ratio /= a + i;
+         if (i < shift_b)
+            ratio /= b + i;
+         result += log(ratio);
+      }
+      return result;
+   }
+
+} // template <class T>T lbeta_imp(T a, T b, const lanczos::undefined_lanczos& l)
+#endif
+
 //
 // Compute the leading power terms in the incomplete Beta:
 //
@@ -1761,6 +1916,31 @@ BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2>::type
    return BOOST_MATH_NAMESPACE::beta(a, b, policies::policy<>());
 }
 
+BOOST_MATH_EXPORT template <class RT1, class RT2, class Policy>
+BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2>::type
+   lbeta(RT1 a, RT2 b, const Policy&)
+{
+   BOOST_FPU_EXCEPTION_GUARD
+   typedef typename tools::promote_args<RT1, RT2>::type result_type;
+   typedef typename policies::evaluation<result_type, Policy>::type value_type;
+   typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
+   typedef typename policies::normalise<
+      Policy,
+      policies::promote_float<false>,
+      policies::promote_double<false>,
+      policies::discrete_quantile<>,
+      policies::assert_undefined<> >::type forwarding_policy;
+
+   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::lbeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), evaluation_type(), forwarding_policy()), "boost::math::lbeta<%1%>(%1%,%1%)");
+}
+
+BOOST_MATH_EXPORT template <class RT1, class RT2>
+BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2>::type
+   lbeta(RT1 a, RT2 b)
+{
+   return BOOST_MATH_NAMESPACE::lbeta(a, b, policies::policy<>());
+}
+
 BOOST_MATH_EXPORT template <class RT1, class RT2, class RT3, class Policy>
 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
    beta(RT1 a, RT2 b, RT3 x, const Policy&)

include/boost/math/special_functions/math_fwd.hpp

@@ -128,6 +128,14 @@ BOOST_MATH_NAMESPACE_BEGIN // Math functions (in roughly alphabetic order).
    BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2, RT3>
          beta(RT1 a, RT2 b, RT3 x, const Policy& pol); // Beta function (3 arguments).
 
+   BOOST_MATH_EXPORT template <class RT1, class RT2>
+   BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2>
+         lbeta(RT1 a, RT2 b); // Log of the beta function.
+
+   BOOST_MATH_EXPORT template <class RT1, class RT2, class Policy>
+   BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2>
+         lbeta(RT1 a, RT2 b, const Policy& pol); // Log of the beta function.
+
    BOOST_MATH_EXPORT template <class RT1, class RT2, class RT3>
    BOOST_MATH_GPU_ENABLED tools::promote_args_t<RT1, RT2, RT3>
          betac(RT1 a, RT2 b, RT3 x);
@@ -1315,6 +1323,10 @@ BOOST_MATH_NAMESPACE_BEGIN // Math functions (in roughly alphabetic order).
    template <class RT1, class RT2, class A>\
    BOOST_MATH_GPU_ENABLED inline BOOST_MATH_NAMESPACE::tools::promote_args_t<RT1, RT2, A> \
    beta(RT1 a, RT2 b, A x){ return ::BOOST_MATH_NAMESPACE::beta(a, b, x, Policy()); }\
+\
+   template <class RT1, class RT2>\
+   BOOST_MATH_GPU_ENABLED inline BOOST_MATH_NAMESPACE::tools::promote_args_t<RT1, RT2> \
+   lbeta(RT1 a, RT2 b) { return ::BOOST_MATH_NAMESPACE::lbeta(a, b, Policy()); }\
 \
    template <class RT1, class RT2, class RT3>\
    BOOST_MATH_GPU_ENABLED inline BOOST_MATH_NAMESPACE::tools::promote_args_t<RT1, RT2, RT3> \

include/boost/math/tools/ulps_plot.hpp

@@ -291,14 +291,18 @@ class ulps_plot {
                     fs << "<circle cx='" << x << "' cy='" << y << "' r='1' fill='" << plot.nan_color_ << "'/>\n";
                     y = y_scale(static_cast<PreciseReal>(-plot.clip_));
                     fs << "<circle cx='" << x << "' cy='" << y << "' r='1' fill='" << plot.nan_color_ << "'/>\n";
+                    continue;
                 }
                 if (plot.clip_ > 0 && static_cast<PreciseReal>(abs(ulp[j])) > plot.clip_)
                 {
-                   if (plot.crop_color_ == "")
+                   // Mark cropped points with a cross, in the color of their function unless one was set,
+                   // so that they stay distinguishable when several functions are plotted:
+                   std::string const & crop_color = plot.crop_color_set_ ? plot.crop_color_ : color;
+                   if (crop_color == "")
                       continue;
                    CoarseReal x = x_scale(plot.coarse_abscissas_[j]);
                    PreciseReal y = y_scale(static_cast<PreciseReal>(ulp[j] < 0 ? -plot.clip_ : plot.clip_));
-                   fs << "<circle cx='" << x << "' cy='" << y << "' r='1' fill='" << plot.crop_color_ << "'/>\n";
+                   fs << "<path d='M" << x - 2.5 << " " << y - 2.5 << " l5 5 m0 -5 l-5 5' stroke='" << crop_color << "' stroke-width='1'/>\n";
                 }
                 else
                 {
@@ -407,6 +411,7 @@ class ulps_plot {
     std::string background_color_;
     std::string font_color_;
     std::string crop_color_;
+    bool crop_color_set_ = false;
     std::string nan_color_;
 };
 
@@ -470,6 +475,7 @@ template<class F, typename PreciseReal, typename CoarseReal>
 ulps_plot<F, PreciseReal, CoarseReal>& ulps_plot<F, PreciseReal, CoarseReal>::crop_color(std::string const & color)
 {
     crop_color_ = color;
+    crop_color_set_ = true;
     return *this;
 }