# HG changeset patch
# User Rik <rik@octave.org>
# Date 1444607717 25200
# Node ID c3c052b9192a73afa67dbfe50a60c8e2f9795393
# Parent  f9c991dc5c1a6d93085719bff5cd8ba3f9a4c315
Improve performance and error reporting of betainv, gaminv (bug #34363).

* betainv.m: Replace for loop with do-until loop.  Shorten max loop cycles to
40, rather than 10,000.  Issue warning if algorithm fails.

* gaminv.m: Replace for loop with do-until loop.  Shorten max loop cycles to
40, rather than 100.  Issue warning if algorithm fails.

diff --git a/scripts/statistics/distributions/betainv.m b/scripts/statistics/distributions/betainv.m
--- a/scripts/statistics/distributions/betainv.m
+++ b/scripts/statistics/distributions/betainv.m
@@ -81,8 +81,10 @@
       y(l) = 1 - sqrt (myeps) * ones (length (l), 1);
     endif
 
-    y_old = y;
-    for i = 1 : 10000
+    y_new = y;
+    loopcnt = 0;
+    do
+      y_old = y_new;
       h     = (betacdf (y_old, a, b) - x) ./ betapdf (y_old, a, b);
       y_new = y_old - h;
       ind   = find (y_new <= myeps);
@@ -94,11 +96,11 @@
         y_new(ind) = 1 - (1 - y_old(ind)) / 10;
       endif
       h = y_old - y_new;
-      if (max (abs (h)) < sqrt (myeps))
-        break;
-      endif
-      y_old = y_new;
-    endfor
+    until (max (abs (h)) < sqrt (myeps) || ++loopcnt == 40)
+
+    if (loopcnt == 40)
+      warning ("betainv: calculation failed to converge for some values");
+    endif
 
     inv(k) = y_new;
   endif
diff --git a/scripts/statistics/distributions/gaminv.m b/scripts/statistics/distributions/gaminv.m
--- a/scripts/statistics/distributions/gaminv.m
+++ b/scripts/statistics/distributions/gaminv.m
@@ -79,8 +79,10 @@
       y(l) = sqrt (myeps) * ones (length (l), 1);
     endif
 
-    y_old = y;
-    for i = 1 : 100
+    y_new = y;
+    loopcnt = 0;
+    do
+      y_old = y_new;
       h     = (gamcdf (y_old, a, b) - x) ./ gampdf (y_old, a, b);
       y_new = y_old - h;
       ind   = find (y_new <= myeps);
@@ -88,13 +90,14 @@
         y_new(ind) = y_old(ind) / 10;
         h = y_old - y_new;
       endif
-      if (max (abs (h)) < sqrt (myeps))
-        break;
-      endif
-      y_old = y_new;
-    endfor
+    until (max (abs (h)) < sqrt (myeps) || ++loopcnt == 40)
+
+    if (loopcnt == 40)
+      warning ("gaminv: calculation failed to converge for some values");
+    endif
 
     inv(k) = y_new;
+
   endif
 
 endfunction