diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 42adebf0..c741b621 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -199,6 +199,7 @@
+
diff --git a/src/Numerics/SpecialFunctions/ExponentialIntegral.cs b/src/Numerics/SpecialFunctions/ExponentialIntegral.cs
new file mode 100644
index 00000000..90b388c7
--- /dev/null
+++ b/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
+ {
+ ///
+ /// Computes the Exponential Integral function.
+ ///
+ /// The argument of the Exponential Integral function.
+ /// The value of the Exponential Integral function.
+ ///
+ /// 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
+ ///
+ ///
+ /// for x > 1 uses continued fraction approach that is often used to compute incomplete gamma.
+ /// for 0 < x <= 1 uses taylor series expansion
+ ///
+ /// Our unit tests suggest that the accuracy of the Exponential Integral function is correct up to 13 floating point digits.
+ ///
+ 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));
+ }
+ }
+ }
+}
diff --git a/src/UnitTests/SpecialFunctionsTests/ExponentialIntegralTests.cs b/src/UnitTests/SpecialFunctionsTests/ExponentialIntegralTests.cs
new file mode 100644
index 00000000..82097ff0
--- /dev/null
+++ b/src/UnitTests/SpecialFunctionsTests/ExponentialIntegralTests.cs
@@ -0,0 +1,53 @@
+
+namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
+{
+ using System;
+ using NUnit.Framework;
+
+ ///
+ /// Factorial tests.
+ ///
+ [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);
+ }
+
+ }
+}
diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj
index e2ca1269..bda8e5d5 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -375,6 +375,7 @@
+