implemented exponential integral for issue #228 included tests for small x, large x, and n=1, 2, and special cases

src/Numerics/Numerics.csproj

@@ -199,6 +199,7 @@
     <Compile Include="RootFinding\Bisection.cs" />
     <Compile Include="SpecialFunctions\Evaluate.cs" />
     <Compile Include="ExcelFunctions.cs" />
+    <Compile Include="SpecialFunctions\ExponentialIntegral.cs" />
     <Compile Include="SpecialFunctions\ModifiedStruve.cs" />
     <Compile Include="SpecialFunctions\ModifiedBessel.cs" />
     <Compile Include="SpecialFunctions\Logistic.cs" />

src/Numerics/SpecialFunctions/ExponentialIntegral.cs

@@ -0,0 +1,108 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+using System.Text;
+
+namespace MathNet.Numerics
+{
+    public static partial class SpecialFunctions
+    {
+        /// <summary>
+        /// Computes the Exponential Integral function.
+        /// </summary>
+        /// <param name="x">The argument of the Exponential Integral function.</param>
+        /// <returns>The value of the Exponential Integral function.</returns>
+        /// <remarks>
+        /// <para>This implementation of the computation of the Exponential Integral function follows the derivation in
+        ///     "Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55", Abramowitz, M., and Stegun, I.A. 1964,  reprinted 1968 by
+        ///     Dover Publications, New York), Chapters 6, 7, and 26.
+        ///     AND
+        ///     "Advanced mathematical methods for scientists and engineers", Bender, Carl M.; Steven A. Orszag (1978). page 253
+        /// </para>
+        /// <para>
+        ///     for x > 1  uses continued fraction approach that is often used to compute incomplete gamma.
+        ///     for 0 < x <= 1 uses taylor series expansion
+        /// </para>
+        /// <para>Our unit tests suggest that the accuracy of the Exponential Integral function is correct up to 13 floating point digits.</para>
+        /// </remarks>
+        public static double ExponentialIntegral(double x, int n)
+        {
+            //parameter validation
+            if (n < 0 || x < 0.0 ) {
+                throw new ArgumentOutOfRangeException(string.Format("x and n must be positive: x={0}, n={1}", x, n));                
+            }
+
+            const double epsilon = 0.00000000000000001;
+            int maxIterations = 100;
+            int i, ii;
+            double ndbl = (double)n;
+            double result = double.NaN;
+            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;
+            double psi;
+            double a, b, c, d, h; //variables for continued fraction
+
+            //special cases
+            if (n == 0)
+            {
+                result = Math.Exp( -1.0d * x ) / x;
+                return result;
+            }
+            else if (x == 0.0d)
+            {
+                result = 1.0d / (ndbl - 1.0d);
+                return result;
+            }
+            //general cases
+            //continued fraction for large x
+            if (x > 1.0d)                
+            { 
+                b = x + ((double)n);
+                c = 1.0d / nearDoubleMin;
+                d = 1.0d / b;
+                h = d;
+                for (i = 1; i <= maxIterations; i++)
+                {
+                    a = -1.0d * ((double)i) * ((ndbl - 1.0d) + (double)i);
+                    b += 2.0d;
+                    d = 1.0d / (a * d + b); 
+                    c = b + a / c;
+                    del = c * d;
+                    h = h * del;
+                    if (Math.Abs(del - 1.0d) < epsilon)
+                    {
+                        result = h * Math.Exp( -x );
+                        return result;
+                    }
+                }
+                throw new ArithmeticException(string.Format("continued fraction failed to converge for x={0}, n={1})", x, n));
+            }
+            //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.
+                for (i = 1; i <= maxIterations; i++)
+                {
+                    factorial *= (-1.0d * x / ((double)i));
+                    if (i != (ndbl - 1.0d)) { del = -factorial / (i - (ndbl - 1.0d)); }
+                    else
+                    {
+                        psi = -1.0d * Constants.EulerMascheroni; 
+                        for (ii = 1; ii <= (ndbl - 1.0d); ii++)
+                        {
+                            psi += (1.0d / ((double)ii));
+                        }
+                        del = factorial * (-1.0d * Math.Log(x) + psi);
+                    }
+                    result += del;
+                    if (Math.Abs(del) < Math.Abs(result) * epsilon)
+                    {
+                        return result;
+                    }
+                }
+                throw new ArithmeticException(string.Format("series failed to converge for x={0}, n={1})", x, n));
+            }
+        }
+    }
+}

src/UnitTests/SpecialFunctionsTests/ExponentialIntegralTests.cs

@@ -0,0 +1,53 @@
+
+namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
+{
+    using System;
+    using NUnit.Framework;
+
+    /// <summary>
+    /// Factorial tests.
+    /// </summary>
+    [TestFixture, Category("Distributions")]
+    public class ExponentialIntegralTests
+    {
+
+        [TestCase(0.001d, 6.33153936413614904)]
+        [TestCase(0.1d, 1.82292395841939059)]
+        [TestCase(1.0d, 0.219383934395520286d)]
+        [TestCase(2.0d, 0.0489005107080611248d)]
+        [TestCase(2.5d, 0.0249149178702697399)]
+        [TestCase(10.0d, 4.15696892968532464e-06)]
+        public void ExponentialIntegral_Matches_MATLAB_and_R_expint_E1(double x, double result)
+        {
+            double actual = SpecialFunctions.ExponentialIntegral( x, 1 );
+            double delta = Math.Abs( result - actual );
+            AssertHelpers.AlmostEqualRelative( result, actual, 13 );
+        }
+
+        [TestCase(0.001d, 2, 0.992668960469238915)]
+        [TestCase(0.1d, 2, 0.722545022194020392)]
+        [TestCase(1.0d, 2, 0.148495506775922048)]
+        [TestCase(2.0d, 2, 0.0375342618204904527)]
+        [TestCase(10.0d, 2, 3.830240465631608e-06)]
+        public void ExponentialIntegral_Matches_R_expint_En(double x, int n, double result)
+        {
+            double actual = SpecialFunctions.ExponentialIntegral(x, n);
+            double delta = Math.Abs(result - actual);
+            AssertHelpers.AlmostEqualRelative(result, actual, 13);
+        }
+
+
+        [TestCase(0.001d, 0, 999.000499833375)]
+        [TestCase(0.1d, 0, 9.048374180359595)]
+        [TestCase(1.0d, 0, 0.3678794411714423)]
+        [TestCase(2.0d, 0, 0.06766764161830635)]
+        [TestCase(10.0d, 0, 4.539992976248485e-06)]
+        public void ExponentialIntegral_SpecialCase_EXP_Matches_from_R_expint_En(double x, int n, double result)
+        {
+            double actual = SpecialFunctions.ExponentialIntegral(x, n);
+            double delta = Math.Abs(result - actual);
+            AssertHelpers.AlmostEqualRelative(result, actual, 13);
+        }
+
+    }
+}

src/UnitTests/UnitTests.csproj

@@ -375,6 +375,7 @@
     <Compile Include="RootFindingTests\NewtonRaphsonTest.cs" />
     <Compile Include="SortingTests.cs" />
     <Compile Include="ExcelTests.cs" />
+    <Compile Include="SpecialFunctionsTests\ExponentialIntegralTests.cs" />
     <Compile Include="SpecialFunctionsTests\ModifiedStruveTests.cs" />
     <Compile Include="SpecialFunctionsTests\ModifiedBesselTests.cs" />
     <Compile Include="SpecialFunctionsTests\ErfTests.cs" />