Bug 23750 – log1p for floats/doubles not actually providing extra accuracy

As discussed in issue 23677 [1], overloads for log1p were recently added for floats and doubles. However, these versions don't seem to actually be providing the additional accuracy that the CEPHES implementation that it is based on provides. These functions forward to a logImpl function that seems to have some the CEPHES code, but doesn't actually produce the more accurate result that CEPHES provides. I adapted the code below based on a CEPHES mirror [2] for doubles. This implementation of log1p provides a much more accurate result than the implementation in phobos (not shown, but comparing to the result before these overloads were added). Moreover, it doesn't experience the same failure when the value of x gets smaller. I compiled with DMD 2.102. ```d enum double SQRTH = 0.70710678118654752440; enum double SQRT2 = 1.41421356237309504880; double log1p(double x) { import std.math.exponential: log; import std.math.algebraic: poly; static immutable double[] LP = [ 2.0039553499201281259648E1, 5.7112963590585538103336E1, 6.0949667980987787057556E1, 2.9911919328553073277375E1, 6.5787325942061044846969E0, 4.9854102823193375972212E-1, 4.5270000862445199635215E-5, ]; static immutable double[] LQ = [ 6.0118660497603843919306E1, 2.1642788614495947685003E2, 3.0909872225312059774938E2, 2.2176239823732856465394E2, 8.3047565967967209469434E1, 1.5062909083469192043167E1, 1.0000000000000000000000E0, ]; double z = 1.0 + x; if ((z < SQRTH) || (z > SQRT2)) return log(z); z = x * x; z = -0.5 * z + x * (z * poly(x, LP) / poly(x, LQ)); return x + z; } void main() { import std.stdio; static import std.math; double x = 1e-15; double y = std.math.log(1 + x); double z1 = log1p(x); double z2 = std.math.log1p(x); double z3 = log1p(x / 10); double z4 = std.math.log1p(x / 10); writefln("the value of x is %.20e", x);//prints around 1.0e-15 writefln("the value of y is %.20e", y); //prints around 1.11e-15 writefln("the value of z1 is %.20e", z1); //prints around 9.99e-16 writefln("the value of z2 is %.20e", z2); //prints around 1.11e-15 writefln("the value of z3 is %.20e", z3); //prints around 1.0e-15 writefln("the value of z4 is %.20e", z4); //prints 0 } ``` [1] https://issues.dlang.org/show_bug.cgi?id=23677 [2] https://github.com/jeremybarnes/cephes/blob/master/cprob/unity.c

@ibuclaw created dlang/phobos pull request #8701 "fix Issue 23750 - log1p for floats/doubles not actually providing extra accuracy" fixing this issue: - fix Issue 23750 - log1p for floats/doubles not actually providing extra accuracy https://github.com/dlang/phobos/pull/8701

dlang/phobos pull request #8701 "fix Issue 23750 - log1p for floats/doubles not actually providing extra accuracy" was merged into stable: - 0a230b0d56120a060a2953eed46a2498b89d26fd by Iain Buclaw: fix Issue 23750 - log1p for floats/doubles not actually providing extra accuracy https://github.com/dlang/phobos/pull/8701

dlang/phobos pull request #8704 "merge stable" was merged into master: - e74d997522f45fd080032112920b88a9f74db32f by Iain Buclaw: fix Issue 23750 - log1p for floats/doubles not actually providing extra accuracy https://github.com/dlang/phobos/pull/8704