diff --git a/src/Numerics.Tests/SpecialFunctionsTests/SpecialFunctionsTests.cs b/src/Numerics.Tests/SpecialFunctionsTests/SpecialFunctionsTests.cs
index b1023717..51107c88 100644
--- a/src/Numerics.Tests/SpecialFunctionsTests/SpecialFunctionsTests.cs
+++ b/src/Numerics.Tests/SpecialFunctionsTests/SpecialFunctionsTests.cs
@@ -311,5 +311,26 @@ namespace MathNet.Numerics.UnitTests.SpecialFunctionsTests
{
AssertHelpers.AlmostEqualRelative(x, SpecialFunctions.Logistic(p), 14);
}
+
+ ///
+ /// Log1p function.
+ ///
+ /// Expected value.
+ /// Input X value.
+ [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);
+ }
}
}
diff --git a/src/Numerics/SpecialFunctions/Log1p.cs b/src/Numerics/SpecialFunctions/Log1p.cs
new file mode 100644
index 00000000..7ed733e2
--- /dev/null
+++ b/src/Numerics/SpecialFunctions/Log1p.cs
@@ -0,0 +1,46 @@
+using System;
+using System.Collections.Generic;
+using System.Linq;
+
+namespace MathNet.Numerics
+{
+ ///
+ /// This partial implementation of the SpecialFunctions class contains all methods related to the log1p function.
+ ///
+ public static partial class SpecialFunctions
+ {
+ ///
+ /// Computes ln(1+x) with good relative precision when |x| is small
+ ///
+ /// The parameter for which to compute the log1p function. Range: x > 0.
+ 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;
+ }
+ }
+}