Merge pull request #267 from glyphard/gh_#228_exponential_integral

Gh #228 exponential integral

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,141 @@
+// <copyright file="ExponentialIntegral.cs" company="Math.NET">
+// 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.
+// </copyright>
+
+// <contribution>
+//    Ashley Messer
+// </contribution>
+
+// ReSharper disable CheckNamespace
+namespace MathNet.Numerics
+// ReSharper restore CheckNamespace
+{
+    using System;
+
+    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>
+    /// Exponential Integral 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

@@ -368,6 +368,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" />