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..228ccf7f
--- /dev/null
+++ b/src/Numerics/SpecialFunctions/ExponentialIntegral.cs
@@ -0,0 +1,141 @@
+//
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2014 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+//
+
+//
+// Ashley Messer
+//
+
+// ReSharper disable CheckNamespace
+namespace MathNet.Numerics
+// ReSharper restore CheckNamespace
+{
+ using System;
+
+ 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..6c37940b
--- /dev/null
+++ b/src/UnitTests/SpecialFunctionsTests/ExponentialIntegralTests.cs
@@ -0,0 +1,53 @@
+
+namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
+{
+ using System;
+ using NUnit.Framework;
+
+ ///
+ /// Exponential Integral 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 1b2a7f83..1836e16d 100644
--- a/src/UnitTests/UnitTests.csproj
+++ b/src/UnitTests/UnitTests.csproj
@@ -368,6 +368,7 @@
+