Portable log1p()

[Francois Grieu]

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 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.


Any comment on the accuracy or portability ?

In particular, is it guaranteed that a conformant
implementation will NOT optimize things, defeating
the attempt to get accurate results ?

If no, is there a way to insure this ?


François Grieu

 

[P.J. Plauger]

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

 

[Francois Grieu]

ran our log1p tests

Testing agaisnt professional, name brand validation
suite; now THAT'S the power of usenet! Much appreciated!

double my_log1p(double x)
{
return log(1+x)-(((1+x)-1)-x)/(1+x);
}

[2043] log1p(-1): got = -nan, want = -inf -- UNEXPECTED
[2046] log1p(+inf): got = -nan, want = +inf -- UNEXPECTED

Indeed, this should not be.

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.)

Speed is very dependent on if log() is in hardware or not.
When it is, my_log1p() is hard to beat in the 0.002 to 0.5
and and -0.25 to -0.002 ranges, I guess.
When log() is in SW, polynomial approximation seems the way
to go.


Meanwhile, a private email pointed out that the above math
is close to the one used by GSL (GNU Scientific Library),

#include <GNU_GPL_v2.h>
double gsl_log1p (const double x)
{
volatile double y;
y = 1 + x;
return log(y)-((y-1)-x)/y ; /* cancels errors with IEEE arithmetic */
}

which however has an interesting extra twist: the use of
"volatile" prevents - I believe/hope somewhat reliably -
the compiler from "simplifying" ((1+x)-1)-x, or performing
similar damage.

I'm still trying to find what the C standards says about how
an implementation is allowed to reshuffle expressions like
log(1+x)-(((1+x)-1)-x)/(1+x)

In also worried of implementations where computations
are carried on e.g 80 bits and rounded to 64 bits when stored.


François Grieu

 

[Kevin D. Quitt]

I'm still trying to find what the C standards says about how
an implementation is allowed to reshuffle expressions like
log(1+x)-(((1+x)-1)-x)/(1+x)

If nothing else, the "as-if" rule applies: it can do anything it wants as
long as the results are identical - as if it made no change.