Special Functions: BetaRegularized always returns result (after reaching eps or 140 iterations)

src/Numerics/SpecialFunctions/Beta.cs

@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2011 Math.NET
+// Copyright (c) 2009-2014 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -86,7 +86,7 @@ namespace MathNet.Numerics
         /// <returns>The lower incomplete (unregularized) beta function.</returns>
         public static double BetaIncomplete(double a, double b, double x)
         {
-            return BetaRegularized(a, b, x) * Beta(a, b);
+            return BetaRegularized(a, b, x)*Beta(a, b);
         }
 
         /// <summary>
@@ -115,15 +115,14 @@ namespace MathNet.Numerics
             }
 
             var bt = (x == 0.0 || x == 1.0)
-                         ? 0.0
+                ? 0.0
-                         : Math.Exp(GammaLn(a + b) - GammaLn(a) - GammaLn(b) + (a * Math.Log(x)) + (b * Math.Log(1.0 - x)));
+                : Math.Exp(GammaLn(a + b) - GammaLn(a) - GammaLn(b) + (a*Math.Log(x)) + (b*Math.Log(1.0 - x)));
 
-            var symmetryTransformation = x >= (a + 1.0) / (a + b + 2.0);
+            var symmetryTransformation = x >= (a + 1.0)/(a + b + 2.0);
 
             /* Continued fraction representation */
-            const int maxIterations = 120;
             var eps = Precision.DoublePrecision;
-            var fpmin = 0.0.Increment() / eps;
+            var fpmin = 0.0.Increment()/eps;
 
             if (symmetryTransformation)
             {
@@ -137,65 +136,60 @@ namespace MathNet.Numerics
             var qap = a + 1.0;
             var qam = a - 1.0;
             var c = 1.0;
-            var d = 1.0 - (qab * x / qap);
+            var d = 1.0 - (qab*x/qap);
 
             if (Math.Abs(d) < fpmin)
             {
                 d = fpmin;
             }
 
-            d = 1.0 / d;
+            d = 1.0/d;
             var h = d;
 
-            for (int m = 1, m2 = 2; m <= maxIterations; m++, m2 += 2)
+            for (int m = 1, m2 = 2; m <= 140; m++, m2 += 2)
             {
-                var aa = m * (b - m) * x / ((qam + m2) * (a + m2));
+                var aa = m*(b - m)*x/((qam + m2)*(a + m2));
-                d = 1.0 + (aa * d);
+                d = 1.0 + (aa*d);
 
                 if (Math.Abs(d) < fpmin)
                 {
                     d = fpmin;
                 }
 
-                c = 1.0 + (aa / c);
+                c = 1.0 + (aa/c);
                 if (Math.Abs(c) < fpmin)
                 {
                     c = fpmin;
                 }
 
-                d = 1.0 / d;
+                d = 1.0/d;
-                h *= d * c;
+                h *= d*c;
-                aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
+                aa = -(a + m)*(qab + m)*x/((a + m2)*(qap + m2));
-                d = 1.0 + (aa * d);
+                d = 1.0 + (aa*d);
 
                 if (Math.Abs(d) < fpmin)
                 {
                     d = fpmin;
                 }
 
-                c = 1.0 + (aa / c);
+                c = 1.0 + (aa/c);
 
                 if (Math.Abs(c) < fpmin)
                 {
                     c = fpmin;
                 }
 
-                d = 1.0 / d;
+                d = 1.0/d;
-                var del = d * c;
+                var del = d*c;
                 h *= del;
 
                 if (Math.Abs(del - 1.0) <= eps)
                 {
-                    if (symmetryTransformation)
+                    return symmetryTransformation ? 1.0 - (bt*h/a) : bt*h/a;
-                    {
-                        return 1.0 - (bt * h / a);
-                    }
-
-                    return bt * h / a;
                 }
             }
 
-            throw new ArgumentException(Resources.ArgumentTooLargeForIterationLimit);
+            return symmetryTransformation ? 1.0 - (bt*h/a) : bt*h/a;
         }
     }
 }