https://issues.dlang.org/show_bug.cgi?id=23750

          Issue ID: 23750
           Summary: log1p for floats/doubles not actually providing extra
                    accuracy
           Product: D
           Version: D2
          Hardware: x86
                OS: Windows
            Status: NEW
          Severity: normal
          Priority: P1
         Component: phobos
          Assignee: nobody@puremagic.com
          Reporter: john.michael.hall@gmail.com

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

--