"Francois Grieu" <fg****@francenet.fr> wrote in message
news:fg**************************@individual.net.. .
Many C platforms lack the C99 log1p(), which computes
log(1+x) accurately even when x is near 0.
I'm looking for a portable alternative.
Somewhere I found code which can be re-expressed as:
double rep_log1p(double x)
{
return ((1+x)-1) ? log(1+x)*(x/((1+x)-1)) : x;
}
but I do not understand the math, and I am worried of
a "dead spot" somewhere.
I ran our log1p tests against this and got:
[2046] log1p(+inf): got = -nan, want = +inf -- UNEXPECTED
log1p: Summary for 1-argument, 53-bit double function:
log1p: Testing 2040 finite function values, result was
log1p: smaller 426 times
log1p: equal 1161 times
log1p: larger 453 times.
log1p: Testing 7 special function values gave 1 mismatches.
log1p: RMS relative bit error was 0.89.
log1p: Max relative bit error was 11.86 for [1001] log1p(+0.000110055)
log1p: Max relative ulp error was 3359 ulp for [1001] log1p(+0.000110055)
64+ 32+ 16+ 8+ 4+ 2+ 1 0 1 2+ 4+ 8+ 16+ 32+ 64+
5 0 10 21 46 139 234 1161 245 106 37 16 9 6 5
I have derived
double my_log1p(double x)
{
return log(1+x)-(((1+x)-1)-x)/(1+x);
}
[the idea here is that if we define
y = 1+x after rounding
e = 1+x-y in the mathematical sense
the derivative of the function log(1+y) is 1/(1+y)
thus log(1+y+e) is close to log(1+y) + e/(1+y) ]
Both functions seem to work, though(my_log1p) seems
somewhat more accurae, and faster.
For this version, our tests gave:
[2043] log1p(-1): got = -nan, want = -inf -- UNEXPECTED
[2046] log1p(+inf): got = -nan, want = +inf -- UNEXPECTED
log1p: Summary for 1-argument, 53-bit double function:
log1p: Testing 2040 finite function values, result was
log1p: smaller 1 times
log1p: equal 2033 times
log1p: larger 6 times.
log1p: Testing 7 special function values gave 2 mismatches.
log1p: RMS relative bit error was -9.00.
log1p: Max relative bit error was 1.00 for [2001] log1p(-2.22045e-16)
log1p: Max relative ulp error was 2 ulp for [2001] log1p(-2.22045e-16)
64+ 32+ 16+ 8+ 4+ 2+ 1 0 1 2+ 4+ 8+ 16+ 32+ 64+
0 0 0 0 0 0 2 2033 4 1 0 0 0 0 0
Any comment on the accuracy or portability ?
Our 53-bit version of log1p gave:
log1p: Summary for 1-argument, 53-bit double function:
log1p: Testing 2040 finite function values, result was
log1p: smaller 15 times
log1p: equal 1998 times
log1p: larger 27 times.
log1p: Testing 7 special function values gave 0 mismatches.
log1p: RMS relative bit error was -7.86.
log1p: Max relative bit error was 1.00 for [2001] log1p(-2.22045e-16)
log1p: Max relative ulp error was 2 ulp for [2001] log1p(-2.22045e-16)
64+ 32+ 16+ 8+ 4+ 2+ 1 0 1 2+ 4+ 8+ 16+ 32+ 64+
0 0 0 0 0 0 16 1998 25 1 0 0 0 0 0
So your function is actually slightly more accurate than ours,
not quite as robust, and probably not quite as fast. (We use
a polynomial approximation over the difficult region.)
BTW, for all tests I used mingw with our log, which is pretty good.
In brief:
log: RMS relative bit error was -9.79.
log: Max relative bit error was 0.90 for [1706] log(+1.70613)
log: Max relative ulp error was 1 ulp for [1048] log(+1.0471)
64+ 32+ 16+ 8+ 4+ 2+ 1 0 1 2+ 4+ 8+ 16+ 32+ 64+
0 0 0 0 0 0 0 2048 2 0 0 0 0 0 0
HTH,
P.J. Plauger
Dinkumware, Ltd.
http://www.dinkumware.com