From 53a464c86c3568248b20387ee0bd8857f558041a Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Sun, 28 Dec 2014 23:00:33 +0100 Subject: [PATCH] Special Functions: document the n-parameter of the exponential integral --- .../SpecialFunctions/ExponentialIntegral.cs | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/Numerics/SpecialFunctions/ExponentialIntegral.cs b/src/Numerics/SpecialFunctions/ExponentialIntegral.cs index 66a159e0..5c4d389d 100644 --- a/src/Numerics/SpecialFunctions/ExponentialIntegral.cs +++ b/src/Numerics/SpecialFunctions/ExponentialIntegral.cs @@ -41,9 +41,10 @@ namespace MathNet.Numerics public static partial class SpecialFunctions { /// - /// Computes the Exponential Integral function. + /// Computes the generalized Exponential Integral function (En). /// /// The argument of the Exponential Integral function. + /// Integer power of the denominator term. Generalization index. /// The value of the Exponential Integral function. /// /// This implementation of the computation of the Exponential Integral function follows the derivation in @@ -70,7 +71,7 @@ namespace MathNet.Numerics int maxIterations = 100; int i, ii; double ndbl = (double)n; - double result = double.NaN; + double result; double nearDoubleMin = 1e-100; //needs a very small value that is not quite as small as the lowest value double can take double factorial = 1.0d; double del; @@ -80,13 +81,11 @@ namespace MathNet.Numerics //special cases if (n == 0) { - result = Math.Exp(-1.0d*x)/x; - return result; + return Math.Exp(-1.0d*x)/x; } else if (x == 0.0d) { - result = 1.0d/(ndbl - 1.0d); - return result; + return 1.0d/(ndbl - 1.0d); } //general cases //continued fraction for large x @@ -106,13 +105,12 @@ namespace MathNet.Numerics h = h*del; if (Math.Abs(del - 1.0d) < epsilon) { - result = h*Math.Exp(-x); - return result; + return h*Math.Exp(-x); } } throw new ArithmeticException(string.Format("continued fraction failed to converge for x={0}, n={1})", x, n)); } - //series computation for small x + //series computation for small x else { result = ((ndbl - 1.0d) != 0 ? 1.0/(ndbl - 1.0d) : (-1.0d*Math.Log(x) - Constants.EulerMascheroni)); //Set first term.