Merge pull request mathnet#887 from delicioustuna/specialfunction-log1p

add log1p function

Author: cdrnet

src/Numerics.Tests/SpecialFunctionsTests/SpecialFunctionsTests.cs

@@ -311,5 +311,26 @@ public void Logistic(double p, double x)
         {
             AssertHelpers.AlmostEqualRelative(x, SpecialFunctions.Logistic(p), 14);
         }
+
+        /// <summary>
+        /// Log1p function.
+        /// </summary>
+        /// <param name="y">Expected value.</param>
+        /// <param name="x">Input X value.</param>
+        [TestCase(-1.0, double.NegativeInfinity)]
+        [TestCase(-0.01, -0.010050335853501442)]
+        [TestCase(-0.0001, -0.00010000500033335834)]
+        [TestCase(-1.0e-6, -1.0000005000003334e-6)]
+        [TestCase(-1.0e-8, -1.0000000050000001e-8)]
+        [TestCase(0.0, 0.0)]
+        [TestCase(1.0e-8, 9.999999950000001e-9)]
+        [TestCase(1.0e-6, 9.999995000003334e-7)]
+        [TestCase(0.0001, 9.999500033330834e-5)]
+        [TestCase(0.01, 0.009950330853168083)]
+        [TestCase(1.0, 0.6931471805599453)]
+        public void Log1p(double x, double y)
+        {
+            AssertHelpers.AlmostEqualRelative(y, SpecialFunctions.Log1p(x), 14);
+        }
     }
 }

src/Numerics/SpecialFunctions/Log1p.cs

@@ -0,0 +1,46 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+
+namespace MathNet.Numerics
+{
+    /// <summary>
+    /// This partial implementation of the SpecialFunctions class contains all methods related to the log1p function.
+    /// </summary>
+    public static partial class SpecialFunctions
+    {
+        /// <summary>
+        /// Computes ln(1+x) with good relative precision when |x| is small
+        /// </summary>
+        /// <param name="x">The parameter for which to compute the log1p function. Range: x > 0.</param>
+        public static double Log1p(double x)
+        {
+            double y0 = Math.Log(1.0 + x);
+
+            if ((-0.2928 < x) && (x < 0.4142))
+            {
+                double y = y0;
+
+                if (y == 0.0)
+                {
+                    y = 1.0;
+                }
+                else if ((y < -0.69) || (y > 0.4))
+                {
+                    y = (Math.Exp(y) - 1.0) / y;
+                }
+                else
+                {
+                    double t = y / 2.0;
+                    y = Math.Exp(t) * Math.Sinh(t) / t;
+                }
+
+                double s = y0 * y;
+                double r = (s - x) / (s + 1.0);
+                y0 = y0 - r * (6 - r) / (6 - 4 * r);
+            }
+
+            return y0;
+        }
+    }
+}