• d-bugmail puremagic.com (77/77) Feb 27 2023 https://issues.dlang.org/show_bug.cgi?id=23750

d-bugmail puremagic.com writes:

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

--