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 (C) 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 possible "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(y) is 1/y
thus log(y+e) is close to log(y) + e/y ]


Both functions seem to work, though my_log1p() seems
somewhat more accurate, and faster.


Any comment on the accuracy or portability ?

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

If no, is there a way to insure this ?


François Grieu

 

[Pierre Asselin]

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 (C) 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 possible "dead spot" somewhere.

This is discussed in David Goldberg, "What Every Computer Scientist
Should Know About Floating-Point Arithmetic", ACM Computing Surveys,
vol. 23 No. 1 pp 5-48 (1 Marh 1991). From the bibliography therein,
the trick comes from the HP-15C calculator.

The idea is, first, that log(1+x) == x*(log(1+x)/x) and, second,
that the function F(u) = log(1+u)/u varies slowly enough that
evaluating it at u == ((1+x)-1) instead of u == x makes no difference.
The log function has to have good relative error, so that the loss
of accuracy in log(1+x) comes from the (1+x) and not the logarithm.
You just need special handling when u= (1+x)-1 == 0, where F(u)= 1.
You do not need a special case when x is large, as the ratio
x/((1+x)-1) approaches unity and drops out of the picture.

You may want special handling for x<0, e.g. use
log(1+x) = -log(1 + (-x/(1+x)) ) .

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(y) is 1/y
thus log(y+e) is close to log(y) + e/y ]

Both functions seem to work, though my_log1p() seems
somewhat more accurate, and faster.

Your formula seems to have all the right properties for x >= 0.

Any comment on the accuracy or portability ?

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

I was wondering about that. I don't have the standard in front of me.

If no, is there a way to insure this ?

IIRC, assign to explicit temporaries.
u= 1+x; u-= 1; return (u ? x*(log(1+u)/u) : x);

 

[Tim Prince]

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 (C) 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 possible "dead spot" somewhere.

This is discussed in David Goldberg, "What Every Computer Scientist
Should Know About Floating-Point Arithmetic", ACM Computing Surveys,
vol. 23 No. 1 pp 5-48 (1 Marh 1991). From the bibliography therein,
the trick comes from the HP-15C calculator.

The idea is, first, that log(1+x) == x*(log(1+x)/x) and, second,
that the function F(u) = log(1+u)/u varies slowly enough that
evaluating it at u == ((1+x)-1) instead of u == x makes no difference.
The log function has to have good relative error, so that the loss
of accuracy in log(1+x) comes from the (1+x) and not the logarithm.
You just need special handling when u= (1+x)-1 == 0, where F(u)= 1.
You do not need a special case when x is large, as the ratio
x/((1+x)-1) approaches unity and drops out of the picture.

You may want special handling for x<0, e.g. use
log(1+x) = -log(1 + (-x/(1+x)) ) .

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(y) is 1/y
thus log(y+e) is close to log(y) + e/y ]

Both functions seem to work, though my_log1p() seems
somewhat more accurate, and faster.

Your formula seems to have all the right properties for x >= 0.

Any comment on the accuracy or portability ?

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

I was wondering about that. I don't have the standard in front of me.

If no, is there a way to insure this ?

IIRC, assign to explicit temporaries.
u= 1+x; u-= 1; return (u ? x*(log(1+u)/u) : x);

François Grieu

Evidently, the two of you agree that "portable" means not depending on
adherence to C99, and one of you considers that it means in accordance with
the state of C prior to C89. I thought I read that some French people held
out for C89 to require observance of parentheses. I'm not certain that
compilers which require special options to observe parentheses may not have
similar problems with temporaries.