# HG changeset patch
# User Nir Krakauer <nkrakauer@ccny.cuny.edu>
# Date 1650406167 25200
# Node ID fe898ae23e0e6ff113d9d51e60b2bd974531c7f5
# Parent  9c7822f4fc9c40405062134f7295930489ac10ee
betainc.m: Use sophisticated technique for calculating exponents to avoid innacuracies (bug #62329)

* betainc.m: Use expm1() and log1p() to replace exponentiation operator ".^" in
special case code when a or b equals 1.  Add BIST tests for bug #62329.

diff -r 9c7822f4fc9c -r fe898ae23e0e scripts/specfun/betainc.m
--- a/scripts/specfun/betainc.m	Sat Apr 16 12:43:24 2022 -0700
+++ b/scripts/specfun/betainc.m	Tue Apr 19 15:09:27 2022 -0700
@@ -139,12 +139,16 @@
 
   if (strcmpi (tail, "lower"))
     y(a_b_one) = x(a_b_one);
-    y(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
+    y(a_one) = - expm1 (log1p (- x(a_one)) .* b(a_one));
     y(b_one) = x(b_one) .^ a(b_one);
   elseif (strcmpi (tail, "upper"))
     y(a_b_one) = 1 - x(a_b_one);
-    y(a_one) = (1 - x(a_one)) .^ b(a_one);
-    y(b_one) = 1 - x(b_one) .^ a(b_one);
+    ## equivalent to "(1 - x(a_one)) .^ b(a_one)", but less roundoff error
+    y(a_one) = exp (log1p (- x(a_one)) .* b(a_one));
+    ## equivalent to "1 - x(b_one) .^ a(b_one)", but less roundoff error
+    y(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)