# HG changeset patch
# User Rik <rik@octave.org>
# Date 1650407431 25200
# Node ID 74f8854958d204299ab88fe023b2471a0fe041a7
# Parent  e7032ef9a9c37d5c45213b03848a94a637f5c758# Parent  fe898ae23e0e6ff113d9d51e60b2bd974531c7f5
maint: merge stable to default.

diff -r e7032ef9a9c3 -r 74f8854958d2 scripts/specfun/betainc.m
--- a/scripts/specfun/betainc.m	Mon Apr 18 20:32:33 2022 -0400
+++ b/scripts/specfun/betainc.m	Tue Apr 19 15:30:31 2022 -0700
@@ -139,12 +139,16 @@
 
   if (strcmpi (tail, "lower"))
     I(a_b_one) = x(a_b_one);
-    I(a_one) = 1 - (1 - x(a_one)) .^ b(a_one);
+    ## See bug #62329.
+    ## equivalent to "1 - (1 - x(a_one)) .^ b(a_one)", but less roundoff error
+    I(a_one) = - expm1 (log1p (- x(a_one)) .* b(a_one));
     I(b_one) = x(b_one) .^ a(b_one);
   elseif (strcmpi (tail, "upper"))
     I(a_b_one) = 1 - x(a_b_one);
-    I(a_one) = (1 - x(a_one)) .^ b(a_one);
-    I(b_one) = 1 - x(b_one) .^ a(b_one);
+    ## equivalent to "(1 - x(a_one)) .^ b(a_one)", but less roundoff error
+    I(a_one) = exp (log1p (- x(a_one)) .* b(a_one));
+    ## equivalent to "1 - x(b_one) .^ a(b_one)", but less roundoff error
+    I(b_one) = - expm1 (log (x(b_one)) .* a(b_one));
   endif
 
   ## Non-Trivial cases
@@ -245,11 +249,19 @@
 %! [a,b] = ndgrid (linspace (1e-4, 100, 20), linspace (1e-4, 100, 20));
 %! assert (betainc (0, a, b), zeros (20));
 %! assert (betainc (1, a, b), ones (20));
+%! assert (betainc (0, a, b, "upper"), ones (20));
+%! assert (betainc (1, a, b, "upper"), zeros (20));
 
 %!test <*34405>
 %! assert (betainc (NaN, 1, 2), NaN);
 %! assert (betainc (0.5, 1, Inf), 1);
 
+%!test <*62329>
+%! assert (betainc (2e-20, 1, 0.5), 1e-20, -1e-15)
+%! assert (betainc (2e-5, 1, 0.5), 2e-5 / (1 + sqrt (1 - 2e-5)), -1e-15)
+%! assert (betainc (0.99, 1, 0.5, "upper"), 0.1, -1e-15)
+%! assert (betainc (0.99, 0.5, 1, "upper"), - expm1 (log (0.99)/2), -1e-15)
+
 ## Test input validation
 %!error <Invalid call> betainc ()
 %!error <Invalid call> betainc (1)