diff --git a/src/Numerics/SpecialFunctions/Beta.cs b/src/Numerics/SpecialFunctions/Beta.cs index ae99548d..9e3064a0 100644 --- a/src/Numerics/SpecialFunctions/Beta.cs +++ b/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 /// The lower incomplete (unregularized) beta function. 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); } /// @@ -115,15 +115,14 @@ namespace MathNet.Numerics } var bt = (x == 0.0 || x == 1.0) - ? 0.0 - : Math.Exp(GammaLn(a + b) - GammaLn(a) - GammaLn(b) + (a * Math.Log(x)) + (b * Math.Log(1.0 - x))); + ? 0.0 + : 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)); - d = 1.0 + (aa * d); + var aa = m*(b - m)*x/((qam + m2)*(a + m2)); + 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; - h *= d * c; - aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2)); - d = 1.0 + (aa * d); + d = 1.0/d; + h *= d*c; + aa = -(a + m)*(qab + m)*x/((a + m2)*(qap + m2)); + 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; - var del = d * c; + d = 1.0/d; + var del = d*c; h *= del; if (Math.Abs(del - 1.0) <= eps) { - if (symmetryTransformation) - { - return 1.0 - (bt * h / a); - } - - return bt * h / a; + return symmetryTransformation ? 1.0 - (bt*h/a) : bt*h/a; } } - throw new ArgumentException(Resources.ArgumentTooLargeForIterationLimit); + return symmetryTransformation ? 1.0 - (bt*h/a) : bt*h/a; } } }