diff --git a/doc/fp_utilities/ulps_plots.qbk b/doc/fp_utilities/ulps_plots.qbk
index 3a6eca5bda..15a70c9fac 100644
--- a/doc/fp_utilities/ulps_plots.qbk
+++ b/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);
         //
diff --git a/doc/sf/beta.qbk b/doc/sf/beta.qbk
index 7e1904c254..4e6a3a24f8 100644
--- a/doc/sf/beta.qbk
+++ b/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
 
diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp
index 8d4ca5bfc5..a89205f4c6 100644
--- a/include/boost/math/special_functions/beta.hpp
+++ b/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&)
diff --git a/include/boost/math/special_functions/math_fwd.hpp b/include/boost/math/special_functions/math_fwd.hpp
index 04e0dad201..6185eda140 100644
--- a/include/boost/math/special_functions/math_fwd.hpp
+++ b/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> \
diff --git a/include/boost/math/tools/ulps_plot.hpp b/include/boost/math/tools/ulps_plot.hpp
index df71510183..2936e5496f 100644
--- a/include/boost/math/tools/ulps_plot.hpp
+++ b/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;
 }
 
diff --git a/reporting/performance/lbeta_performance.cpp b/reporting/performance/lbeta_performance.cpp
new file mode 100644
index 0000000000..d17da2dc17
--- /dev/null
+++ b/reporting/performance/lbeta_performance.cpp
@@ -0,0 +1,86 @@
+//  (C) Copyright Nick Thompson 2026.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+//
+// Compares lbeta(a, b) with the two obvious alternatives, log(beta(a, b)) and
+// lgamma(a) + lgamma(b) - lgamma(a + b), with a and b log-uniform in [1e-3, 1e6]:
+// the range where neither alternative fails outright for most pairs.
+//
+// Promotion is disabled so that each type is timed in its own arithmetic: by default,
+// float is evaluated in double, and double in long double where that is wider.
+//
+
+#include <cmath>
+#include <random>
+#include <utility>
+#include <vector>
+#include <benchmark/benchmark.h>
+#include <boost/math/special_functions/beta.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+
+using no_promotion = boost::math::policies::policy<boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>>;
+
+template <typename Real>
+std::vector<std::pair<Real, Real>> arguments()
+{
+    std::mt19937_64 mt {42};
+    std::uniform_real_distribution<double> dist(-3, 6);
+    std::vector<std::pair<Real, Real>> v(1024);
+    for (auto& p : v)
+    {
+        p = {static_cast<Real>(std::pow(10.0, dist(mt))), static_cast<Real>(std::pow(10.0, dist(mt)))};
+    }
+    return v;
+}
+
+template <typename Real>
+void lbeta_performance(benchmark::State& state)
+{
+    auto v = arguments<Real>();
+    std::size_t i = 0;
+    for (auto _ : state)
+    {
+        auto const & p = v[i++ & 1023];
+        benchmark::DoNotOptimize(boost::math::lbeta(p.first, p.second, no_promotion()));
+    }
+}
+
+template <typename Real>
+void log_beta_performance(benchmark::State& state)
+{
+    using std::log;
+    auto v = arguments<Real>();
+    std::size_t i = 0;
+    for (auto _ : state)
+    {
+        auto const & p = v[i++ & 1023];
+        benchmark::DoNotOptimize(log(boost::math::beta(p.first, p.second, no_promotion())));
+    }
+}
+
+template <typename Real>
+void lgamma_sum_performance(benchmark::State& state)
+{
+    using boost::math::lgamma;
+    auto v = arguments<Real>();
+    std::size_t i = 0;
+    for (auto _ : state)
+    {
+        auto const & p = v[i++ & 1023];
+        benchmark::DoNotOptimize(lgamma(p.first, no_promotion()) + lgamma(p.second, no_promotion()) - lgamma(p.first + p.second, no_promotion()));
+    }
+}
+
+BENCHMARK_TEMPLATE(lbeta_performance, float);
+BENCHMARK_TEMPLATE(lbeta_performance, double);
+BENCHMARK_TEMPLATE(lbeta_performance, long double);
+BENCHMARK_TEMPLATE(log_beta_performance, float);
+BENCHMARK_TEMPLATE(log_beta_performance, double);
+BENCHMARK_TEMPLATE(log_beta_performance, long double);
+BENCHMARK_TEMPLATE(lgamma_sum_performance, float);
+BENCHMARK_TEMPLATE(lgamma_sum_performance, double);
+BENCHMARK_TEMPLATE(lgamma_sum_performance, long double);
+
+BENCHMARK_MAIN();
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index b13093a5b1..e5de8e296d 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -217,6 +217,7 @@ test-suite special_fun :
    [ run test_bessel_k_prime.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework  ]
    [ run bessel_iterator_test.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework  ]
    [ run test_beta.cpp test_instances//test_instances pch_light /boost/test//boost_unit_test_framework  ]
+   [ run test_lbeta.cpp ]
    [ run test_bessel_airy_zeros.cpp /boost/test//boost_unit_test_framework  ]
    [ run test_bernoulli_constants.cpp /boost/test//boost_unit_test_framework  ]
    [ run test_binomial_coeff.cpp pch /boost/test//boost_unit_test_framework  ]
diff --git a/test/compile_test/instantiate.hpp b/test/compile_test/instantiate.hpp
index 96c8964866..3eaeed5900 100644
--- a/test/compile_test/instantiate.hpp
+++ b/test/compile_test/instantiate.hpp
@@ -275,6 +275,7 @@ void instantiate(RealType)
    boost::math::erfc_inv(v1);
    boost::math::beta(v1, v2);
    boost::math::beta(v1, v2, v3);
+   boost::math::lbeta(v1, v2);
    boost::math::betac(v1, v2, v3);
    boost::math::ibeta(v1, v2, v3);
    boost::math::ibetac(v1, v2, v3);
@@ -556,6 +557,7 @@ void instantiate(RealType)
    boost::math::erfc_inv(v1 * 1);
    boost::math::beta(v1 * 1, v2 + 0);
    boost::math::beta(v1 * 1, v2 + 0, v3 / 1);
+   boost::math::lbeta(v1 * 1, v2 + 0);
    boost::math::betac(v1 * 1, v2 + 0, v3 / 1);
    boost::math::ibeta(v1 * 1, v2 + 0, v3 / 1);
    boost::math::ibetac(v1 * 1, v2 + 0, v3 / 1);
@@ -809,6 +811,7 @@ void instantiate(RealType)
    boost::math::erfc_inv(v1, pol);
    boost::math::beta(v1, v2, pol);
    boost::math::beta(v1, v2, v3, pol);
+   boost::math::lbeta(v1, v2, pol);
    boost::math::betac(v1, v2, v3, pol);
    boost::math::ibeta(v1, v2, v3, pol);
    boost::math::ibetac(v1, v2, v3, pol);
@@ -1088,6 +1091,7 @@ void instantiate(RealType)
    test::erfc_inv(v1);
    test::beta(v1, v2);
    test::beta(v1, v2, v3);
+   test::lbeta(v1, v2);
    test::betac(v1, v2, v3);
    test::ibeta(v1, v2, v3);
    test::ibetac(v1, v2, v3);
@@ -1375,6 +1379,7 @@ void instantiate_mixed(RealType)
    }
    boost::math::erfc_inv(i);
    boost::math::beta(i, s);
+   boost::math::lbeta(i, s);
    boost::math::beta(fr, lr);
    boost::math::beta(i, s, l);
    boost::math::beta(fr, dr, lr);
@@ -1592,6 +1597,7 @@ void instantiate_mixed(RealType)
    }
    boost::math::erfc_inv(i, pol);
    boost::math::beta(i, s, pol);
+   boost::math::lbeta(i, s, pol);
    boost::math::beta(fr, lr, pol);
    boost::math::beta(i, s, l, pol);
    boost::math::beta(fr, dr, lr, pol);
@@ -1808,6 +1814,8 @@ void instantiate_mixed(RealType)
    test::erfc_inv(i);
    test::beta(i, s);
    test::beta(fr, lr);
+   test::lbeta(i, s);
+   test::lbeta(fr, lr);
    test::beta(i, s, l);
    test::beta(fr, dr, lr);
    test::betac(l, i, s);
diff --git a/test/compile_test/sf_beta_incl_test.cpp b/test/compile_test/sf_beta_incl_test.cpp
index e20431c67a..53a0804de4 100644
--- a/test/compile_test/sf_beta_incl_test.cpp
+++ b/test/compile_test/sf_beta_incl_test.cpp
@@ -21,6 +21,12 @@ void compile_and_link_test()
    check_result<long double>(boost::math::beta<long double>(l, l));
 #endif
 
+   check_result<float>(boost::math::lbeta<float, float>(f, f));
+   check_result<double>(boost::math::lbeta<double>(d, d));
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+   check_result<long double>(boost::math::lbeta<long double>(l, l));
+#endif
+
    check_result<float>(boost::math::ibeta<float>(f, f, f));
    check_result<double>(boost::math::ibeta<double>(d, d, d));
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
diff --git a/test/cuda_jamfile b/test/cuda_jamfile
index 3f717bf36e..3b835c921f 100644
--- a/test/cuda_jamfile
+++ b/test/cuda_jamfile
@@ -261,6 +261,8 @@ run test_airy_bi_prime_float.cu ;
 
 run test_beta_double.cu ;
 run test_beta_float.cu ;
+run test_lbeta_double.cu ;
+run test_lbeta_float.cu ;
 run test_betac_double.cu ;
 run test_betac_float.cu ;
 run test_ibeta_double.cu ;
diff --git a/test/test_lbeta.cpp b/test/test_lbeta.cpp
new file mode 100644
index 0000000000..28e40fd13f
--- /dev/null
+++ b/test/test_lbeta.cpp
@@ -0,0 +1,146 @@
+// Copyright Nick Thompson, 2026
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include "math_unit_test.hpp"
+#include <cmath>
+#include <stdexcept>
+#include <boost/math/concepts/real_concept.hpp>
+#include <boost/math/special_functions/beta.hpp>
+#include <boost/math/special_functions/trunc.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+
+using boost::math::lbeta;
+using boost::multiprecision::cpp_bin_float_100;
+
+// a = am * 2^ae and b = bm * 2^be are exact in every type whose exponent range holds them.
+// The values of log(beta(a, b)) were computed as lgamma(a) + lgamma(b) - lgamma(a + b) with
+// 800-digit cpp_bin_float, enough to absorb the cancellation when an argument is 2^1000.
+struct lbeta_case
+{
+   double am;
+   int ae;
+   double bm;
+   int be;
+   const char* value;
+};
+
+const lbeta_case cases[] = {
+      { 1, -60, 1, -50, "4.15898069195697740239298881123986100039915095416801455108115264409117033485238337627682905027475554197949082691e+01" },
+      { 1, -100, 3, 0, "6.93147180559945309417232121458164735161921819183174078192200771589333072536381647225831806423195466923922439014e+01" },
+      { 0.5, 0, 0.5, 0, "1.14472988584940017414342735135305871164729481291531157151362307147213776988482607978362327027548970770200981223e+00" },
+      { 0.25, 0, 0.75, 0, "1.49130347612937282885204341208214699568504488009543919857396307621883458086967343758655493377369905147301055274e+00" },
+      { 1.5, 0, 2.5, 0, "-1.62785883639038106352550113447964756065470572452570944496909696650143671799395278263983003771018504246599611185e+00" },
+      { 2, 0, 0.75, 0, "-2.71933715483641758831669494532999161982574749963589623711364445601499668091899372237820945415928865613025799190e-01" },
+      { 2.5, 0, 2.5, 0, "-2.60868808940210730038195226193165156023371556978372575559644266134412329068442796258380426388570901630403052585e+00" },
+      { 2.5, 0, 7.25, 0, "-4.90533661883930395460834936695643873217655041103926652565285221548103484033393583597181887094577679099681463945e+00" },
+      { 7.25, 0, 10, 0, "-1.15206093828576774238792925676023403474351969826024007014606103008266556359021233866230652434846476474713883000e+01" },
+      { 33, 0, 70, 0, "-6.52309611058463901652530664244886314226939493764776246951104105272718684570163848327941904447963477233497622302e+01" },
+      { 1000, 0, 2.5, 0, "-1.69865790781530187420413395240576683950617459492281475607293952992950240177344530711756692292901949091581441171e+01" },
+      { 1, 20, 2.5, 0, "-3.43726779456625527055870837198099350327971371888199561621811864168868157894165670678814799899328577816876559438e+01" },
+      { 1, 40, 0.125, 0, "-1.44631754524588046377436687345621730310702661961084719229448263705785310178427063985532264153346167268065642354e+00" },
+      { 1, 60, 3, 0, "-1.24073345320230210388286634954978816325251619811257890584231209580161422603673162680077814295275683492847917281e+02" },
+      { 7.25, 100, 1, -3, "-6.89254658305384294843871197074079636920399542538686737751765487854131887259626569180127928805023089821623285643e+00" },
+      { 1, 20, 1, 20, "-1.45364066196521333105311408462318440517413066016099383479834829660518401041863896220590715872159106794242688144e+06" },
+      { 1, 60, 1.5, 60, "-1.93982405936568553724082316071548785790406066925570806020585049192655840461166500901608060631298307929683783052e+18" },
+      { 1, 100, 1, 90, "-9.81929078485523071471385110686110593997035850585013134716105089872660111623496248017922967392287894559854586350e+27" },
+      { 2.5, 200, 2.5, 0, "-3.48579634239185123211942384088825609792055182534338388500569779067988144304556260885743443173501924555859593246e+02" },
+      { 1, -1000, 1, 1000, "6.93147180559945309417232121458176568075500134360255254120680009493393621969694715605863326996418687542001481021e+02" },
+      { 1, 1000, 1, 1000, "-1.48542634003375029367139526623878447128207418381951509618626907063551595543031141778517375655554942542595596811e+301" },
+      { 1, -1000, 1, -1000, "6.93840327740505254726649353579634744643575634494615509374800689502887015591664410321469190323415106229543482502e+02" },
+      { 1.5, 500, 7.25, 3, "-1.99483893641016901555392148959348216596399017908133451106724466122953422408704368366657864286454193673764307194e+04" },
+      { 1, 996, 1, -53, "3.67368005696770246880527428609492912902170429819419624345920988651297592807040422075224154620060976686905127353e+01" }
+};
+
+template <class Real>
+Real from_reference(const cpp_bin_float_100& x)
+{
+   return static_cast<Real>(x);
+}
+
+#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
+template <>
+boost::math::concepts::real_concept from_reference<boost::math::concepts::real_concept>(const cpp_bin_float_100& x)
+{
+   return boost::math::concepts::real_concept(static_cast<long double>(x));
+}
+#endif
+
+// Where |log(beta)| < 1, it passes through zero and relative error is meaningless, so check the absolute error.
+template <class Real>
+void test_spots(std::size_t ulps, Real absolute_tolerance)
+{
+   using std::abs;
+   using std::ldexp;
+   const int max_exponent = boost::math::itrunc(boost::math::tools::log_max_value<Real>() / boost::math::constants::ln_two<Real>()) - 2;
+   for (const lbeta_case& c : cases)
+   {
+      if ((abs(c.ae) > max_exponent) || (abs(c.be) > max_exponent))
+         continue;
+      Real a = ldexp(Real(c.am), c.ae);
+      Real b = ldexp(Real(c.bm), c.be);
+      Real expected = from_reference<Real>(cpp_bin_float_100(c.value));
+      Real computed = lbeta(a, b);
+      if (abs(expected) >= 1)
+         CHECK_ULP_CLOSE(expected, computed, ulps);
+      else
+         CHECK_LE(Real(abs(expected - computed)), absolute_tolerance);
+      CHECK_EQUAL(computed, Real(lbeta(b, a)));
+   }
+   // With one argument equal to 1, log(beta) is -log of the other. Compare with a high-precision
+   // reference rather than Real's own log, which real_concept evaluates only in double precision:
+   CHECK_ULP_CLOSE(from_reference<Real>(cpp_bin_float_100(-log(cpp_bin_float_100(3)))), Real(lbeta(Real(1), Real(3))), 1);
+   CHECK_ULP_CLOSE(from_reference<Real>(cpp_bin_float_100(-log(cpp_bin_float_100(0.25)))), Real(lbeta(Real(0.25), Real(1))), 1);
+#ifndef BOOST_NO_EXCEPTIONS
+   CHECK_THROW(lbeta(Real(0), Real(1)), std::domain_error);
+   CHECK_THROW(lbeta(Real(1), Real(-2)), std::domain_error);
+#endif
+   // NaN arguments are domain errors; ignoring the error gives NaN:
+   if (std::numeric_limits<Real>::has_quiet_NaN)
+   {
+      using ignore_domain = boost::math::policies::policy<boost::math::policies::domain_error<boost::math::policies::ignore_error>>;
+      const Real nan = std::numeric_limits<Real>::quiet_NaN();
+      CHECK_TRUE((boost::math::isnan)(Real(lbeta(nan, Real(2), ignore_domain()))));
+      CHECK_TRUE((boost::math::isnan)(Real(lbeta(Real(2), nan, ignore_domain()))));
+      CHECK_TRUE((boost::math::isnan)(Real(lbeta(nan, Real(0.5), ignore_domain()))));
+      CHECK_TRUE((boost::math::isnan)(Real(lbeta(Real(0.5), nan, ignore_domain()))));
+      CHECK_TRUE((boost::math::isnan)(Real(lbeta(nan, Real(50), ignore_domain()))));
+      CHECK_TRUE((boost::math::isnan)(Real(lbeta(nan, nan, ignore_domain()))));
+#ifndef BOOST_NO_EXCEPTIONS
+      CHECK_THROW(lbeta(nan, Real(2)), std::domain_error);
+      CHECK_THROW(lbeta(Real(2), nan), std::domain_error);
+#endif
+   }
+   // log(beta) tends to -infinity with either argument, reported as an overflow:
+   if (std::numeric_limits<Real>::has_infinity)
+   {
+      using ignore_overflow = boost::math::policies::policy<boost::math::policies::overflow_error<boost::math::policies::ignore_error>>;
+      const Real infinity = std::numeric_limits<Real>::infinity();
+      CHECK_EQUAL(Real(lbeta(infinity, Real(2.5), ignore_overflow())), Real(-infinity));
+      CHECK_EQUAL(Real(lbeta(Real(1e-30), infinity, ignore_overflow())), Real(-infinity));
+      CHECK_EQUAL(Real(lbeta(infinity, infinity, ignore_overflow())), Real(-infinity));
+#ifndef BOOST_NO_EXCEPTIONS
+      CHECK_THROW(lbeta(infinity, Real(2.5)), std::overflow_error);
+#endif
+   }
+}
+
+int main()
+{
+   test_spots<float>(2, 2 * std::numeric_limits<float>::epsilon());
+   test_spots<double>(4, 2 * std::numeric_limits<double>::epsilon());
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+   test_spots<long double>(8, 4 * std::numeric_limits<long double>::epsilon());
+#endif
+#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
+   test_spots<boost::math::concepts::real_concept>(8, boost::math::concepts::real_concept(16 * std::numeric_limits<long double>::epsilon()));
+#endif
+   using boost::multiprecision::cpp_bin_float_50;
+   test_spots<cpp_bin_float_50>(8, 8 * std::numeric_limits<cpp_bin_float_50>::epsilon());
+   // No Lanczos approximation at this precision, so this exercises the generic code path, whose
+   // errors come from beta, lgamma and tgamma_delta_ratio at this precision:
+   test_spots<cpp_bin_float_100>(64, 64 * std::numeric_limits<cpp_bin_float_100>::epsilon());
+   return boost::math::test::report_errors();
+}
diff --git a/test/test_lbeta_double.cu b/test/test_lbeta_double.cu
new file mode 100644
index 0000000000..c6f84fb717
--- /dev/null
+++ b/test/test_lbeta_double.cu
@@ -0,0 +1,139 @@
+//  Copyright John Maddock 2016.
+//  Copyright Nick Thompson 2026.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
+#define BOOST_MATH_PROMOTE_DOUBLE_POLICY false
+
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <cmath>
+#include <limits>
+#include <boost/math/special_functions/beta.hpp>
+#include <boost/math/special_functions/relative_difference.hpp>
+#include <boost/array.hpp>
+#include "cuda_managed_ptr.hpp"
+#include "stopwatch.hpp"
+
+// For the CUDA runtime routines (prefixed with "cuda_")
+#include <cuda_runtime.h>
+
+typedef double float_type;
+
+/**
+ * CUDA Kernel Device code
+ *
+ */
+__global__ void cuda_test(const float_type *in1, const float_type * in2, float_type *out, int numElements)
+{
+    using std::cos;
+    int i = blockDim.x * blockIdx.x + threadIdx.x;
+
+    if (i < numElements)
+    {
+        out[i] = boost::math::lbeta(in1[i], in2[i]);
+    }
+}
+
+template <class T> struct table_type { typedef T type; };
+typedef float_type T;
+#define SC_(x) static_cast<T>(x)
+
+#include "beta_med_data.ipp"
+#include "beta_small_data.ipp"
+
+/**
+ * Host main routine
+ */
+int main(void)
+{
+  try{
+    // Consolidate the test data:
+    std::vector<float_type> v1, v2;
+
+    for(unsigned i = 0; i < beta_med_data.size(); ++i)
+    {
+       v1.push_back(beta_med_data[i][0]);
+       v2.push_back(beta_med_data[i][1]);
+    }
+    for(unsigned i = 0; i < beta_small_data.size(); ++i)
+    {
+       v1.push_back(beta_small_data[i][0]);
+       v2.push_back(beta_small_data[i][1]);
+    }
+
+    // Error code to check return values for CUDA calls
+    cudaError_t err = cudaSuccess;
+
+    // Print the vector length to be used, and compute its size
+    int numElements = 50000;
+    std::cout << "[Vector operation on " << numElements << " elements]" << std::endl;
+
+    // Allocate the managed input vector A
+    cuda_managed_ptr<float_type> input_vector1(numElements);
+    cuda_managed_ptr<float_type> input_vector2(numElements);
+
+    // Allocate the managed output vector C
+    cuda_managed_ptr<float_type> output_vector(numElements);
+
+    // Initialize the input vectors
+    for (int i = 0; i < numElements; ++i)
+    {
+        int table_id = i % v1.size();
+        input_vector1[i] = v1[table_id];
+        input_vector2[i] = v2[table_id];
+    }
+
+    // Launch the Vector Add CUDA Kernel
+    int threadsPerBlock = 1024;
+    int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock;
+    std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl;
+
+    watch w;
+    cuda_test<<<blocksPerGrid, threadsPerBlock>>>(input_vector1.get(), input_vector2.get(), output_vector.get(), numElements);
+    cudaDeviceSynchronize();
+    std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl;
+
+    err = cudaGetLastError();
+    if (err != cudaSuccess)
+    {
+        std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    // Verify that the result vector is correct
+    std::vector<float_type> results;
+    results.reserve(numElements);
+    w.reset();
+    for(int i = 0; i < numElements; ++i)
+       results.push_back(boost::math::lbeta(input_vector1[i], input_vector2[i]));
+    double t = w.elapsed();
+    // check the results
+    for(int i = 0; i < numElements; ++i)
+    {
+        // log(beta) passes through zero, where relative error is meaningless, so compare absolutely there:
+        float_type error = std::fabs(results[i]) < 1
+            ? std::fabs(output_vector[i] - results[i]) / std::numeric_limits<float_type>::epsilon()
+            : boost::math::epsilon_difference(output_vector[i], results[i]);
+        if (error > 300)
+        {
+            std::cerr << "Result verification failed at element " << i << "!" << std::endl;
+            std::cerr << "Error rate was: " << error << "eps" << std::endl;
+            return EXIT_FAILURE;
+        }
+    }
+
+    std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl;
+    std::cout << "Done\n";
+  }
+  catch(const std::exception& e)
+  {
+    std::cerr << "Stopped with exception: " << e.what() << std::endl;
+  }
+  return 0;
+}
+
+
diff --git a/test/test_lbeta_float.cu b/test/test_lbeta_float.cu
new file mode 100644
index 0000000000..6fc05fa05c
--- /dev/null
+++ b/test/test_lbeta_float.cu
@@ -0,0 +1,137 @@
+//  Copyright John Maddock 2016.
+//  Copyright Nick Thompson 2026.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
+
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <cmath>
+#include <limits>
+#include <boost/math/special_functions/beta.hpp>
+#include <boost/math/special_functions/relative_difference.hpp>
+#include <boost/array.hpp>
+#include "cuda_managed_ptr.hpp"
+#include "stopwatch.hpp"
+
+// For the CUDA runtime routines (prefixed with "cuda_")
+#include <cuda_runtime.h>
+
+typedef float float_type;
+
+/**
+ * CUDA Kernel Device code
+ *
+ */
+__global__ void cuda_test(const float_type *in1, const float_type * in2, float_type *out, int numElements)
+{
+    using std::cos;
+    int i = blockDim.x * blockIdx.x + threadIdx.x;
+
+    if (i < numElements)
+    {
+        out[i] = boost::math::lbeta(in1[i], in2[i]);
+    }
+}
+
+template <class T> struct table_type { typedef T type; };
+typedef float_type T;
+#define SC_(x) static_cast<T>(x)
+
+#include "beta_med_data.ipp"
+#include "beta_small_data.ipp"
+
+/**
+ * Host main routine
+ */
+int main(void)
+{
+  try{
+    // Consolidate the test data:
+    std::vector<float_type> v1, v2;
+
+    for(unsigned i = 0; i < beta_med_data.size(); ++i)
+    {
+       v1.push_back(beta_med_data[i][0]);
+       v2.push_back(beta_med_data[i][1]);
+    }
+    for(unsigned i = 0; i < beta_small_data.size(); ++i)
+    {
+       v1.push_back(beta_small_data[i][0]);
+       v2.push_back(beta_small_data[i][1]);
+    }
+
+    // Error code to check return values for CUDA calls
+    cudaError_t err = cudaSuccess;
+
+    // Print the vector length to be used, and compute its size
+    int numElements = 50000;
+    std::cout << "[Vector operation on " << numElements << " elements]" << std::endl;
+
+    // Allocate the managed input vector A
+    cuda_managed_ptr<float_type> input_vector1(numElements);
+    cuda_managed_ptr<float_type> input_vector2(numElements);
+
+    // Allocate the managed output vector C
+    cuda_managed_ptr<float_type> output_vector(numElements);
+
+    // Initialize the input vectors
+    for (int i = 0; i < numElements; ++i)
+    {
+        int table_id = i % v1.size();
+        input_vector1[i] = v1[table_id];
+        input_vector2[i] = v2[table_id];
+    }
+
+    // Launch the Vector Add CUDA Kernel
+    int threadsPerBlock = 1024;
+    int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock;
+    std::cout << "CUDA kernel launch with " << blocksPerGrid << " blocks of " << threadsPerBlock << " threads" << std::endl;
+
+    watch w;
+    cuda_test<<<blocksPerGrid, threadsPerBlock>>>(input_vector1.get(), input_vector2.get(), output_vector.get(), numElements);
+    cudaDeviceSynchronize();
+    std::cout << "CUDA kernal done in " << w.elapsed() << "s" << std::endl;
+
+    err = cudaGetLastError();
+    if (err != cudaSuccess)
+    {
+        std::cerr << "Failed to launch vectorAdd kernel (error code " << cudaGetErrorString(err) << ")!" << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    // Verify that the result vector is correct
+    std::vector<float_type> results;
+    results.reserve(numElements);
+    w.reset();
+    for(int i = 0; i < numElements; ++i)
+       results.push_back(boost::math::lbeta(input_vector1[i], input_vector2[i]));
+    double t = w.elapsed();
+    // check the results
+    for(int i = 0; i < numElements; ++i)
+    {
+        // log(beta) passes through zero, where relative error is meaningless, so compare absolutely there:
+        float_type error = std::fabs(results[i]) < 1
+            ? std::fabs(output_vector[i] - results[i]) / std::numeric_limits<float_type>::epsilon()
+            : boost::math::epsilon_difference(output_vector[i], results[i]);
+        if (error > 300)
+        {
+            std::cerr << "Result verification failed at element " << i << "!" << std::endl;
+            std::cerr << "Error rate was: " << error << "eps" << std::endl;
+            return EXIT_FAILURE;
+        }
+    }
+
+    std::cout << "Test PASSED with calculation time: " << t << "s" << std::endl;
+    std::cout << "Done\n";
+  }
+  catch(const std::exception& e)
+  {
+    std::cerr << "Stopped with exception: " << e.what() << std::endl;
+  }
+  return 0;
+}
+