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: nob...@puremagic.com
Reporter: john.michael.h...@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
--
